Aims: Linear magnetohydrostatic (MHS) models of solar magnetic fields balance plasma pressure gradients, gravity and Lorentz forces where the current density is composed of a linear force-free component and a cross-field component that depends on gravitational stratification. In this paper, we investigate an efficient numerical procedure for calculating such equilibria. Methods: The MHS equations are reduced to two scalar elliptic equations – one on the lower boundary and the other within the interior of the computational domain. The normal component of the magnetic field is prescribed on the lower boundary and a multigrid method is applied on both this boundary and within the domain to find the poloidal scalar potential. Once solved to a desired accuracy, the magnetic field, plasma pressure and density are found using a finite difference method. Results: We investigate the effects of the cross-field currents on the linear MHS equilibria. Force-free and non-force-free examples are given to demonstrate the numerical scheme and an analysis of speed-up due to parallelization on a graphics processing unit (GPU) is presented. It is shown that speed-ups of ×30 are readily achievable.