psol_rbgpre Subroutine

   subroutine psol_rbgpre(b, bd, ipbd)
   !! This routine solves a linear system \(A_R x = b\), using the LU factors of the diagonal 
   !! blocks computed in [[jac_rbgpre]] and mesh parameters passed as module variables. The
   !! solution is carried out by [[dgesl]] from LINPACK.
      
      use linpack, only: dgesl

      real(rk), intent(inout) :: b(*)
         !! Right-hand side vector on entry and solution vector on return.
      real(rk), intent(in) :: bd(*)
         !! LU factors of the diagonal blocks. `bd` corresponds to the segment `rwp` of `rwork`.
      integer, intent(in) :: ipbd(*)
         !! Pivots for the LU factorizations. `ipbd` corresponds to the segment `iwp` of `iwork`.

      integer :: ib, ibd, ig0, igm1, igx, igy, iip, jx, jy

      ib = 1
      do jy = 1, meshy
         igy = jigy(jy)
         ig0 = (igy - 1)*ngx
         do jx = 1, meshx
            igx = jigx(jx)
            igm1 = igx - 1 + ig0
            ibd = 1 + igm1*mpsq
            iip = 1 + igm1*mp
            call dgesl(bd(ibd), mp, mp, ipbd(iip), b(ib), 0)
            ib = ib + mp
         end do
      end do

   end subroutine psol_rbgpre