  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 x30
are readily achievable.
