psol_ilupre Subroutine

   subroutine psol_ilupre( &
      neq, t, y, ydot, r0, wk, cj, wght, rwp, iwp, b, epslin, ierr, rpar, ipar)
   !! This subroutine solves the linear system \(P x = b\) for the sparse preconditioner \(P\),
   !! given a vector \(b\), using the incomplete LU decomposition produced by [[jac_ilupre]].
   !! The solution is carried out by [[lusol]] from SPARSKIT.

      use dsparskit, only: lusol, dvperm

      integer, intent(in) :: neq
         !! Problem size.
      real(rk), intent(in) :: t
         !! Current independent variable (not used).
      real(rk), intent(in) :: y(*)
         !! Current dependent variables (not used).
      real(rk), intent(in) :: ydot(*)
         !! Current derivatives of dependent variables (not used).
      real(rk), intent(in) :: r0(*)
         !! Current residual evaluated at `(t, y, ydot)` (not used).
      real(rk), intent(in) :: wk(*)
         !! Real work space available to this subroutine.
      real(rk), intent(in) :: cj
         !! Scalar used in forming the system Jacobian.
      real(rk), intent(in) :: wght(*)
         !! Error weights for computing norms.
      real(rk), intent(inout) :: rwp(*)
         !! Matrix elements of ILU.
      integer, intent(inout) :: iwp(*)
         !! Array indices for elements of ILU
      real(rk), intent(inout) :: b(*)
         !! Right-hand side vector on input; solution on output.
      real(rk), intent(in) :: epslin
         !! Tolerance for linear system (not used).
      integer, intent(inout) :: ierr
         !! Error flag (not used).
      real(rk), intent(inout) :: rpar(*)
         !! Real array used for communication between the calling program and user routines
         !! (not used).
      integer, intent(inout) :: ipar(*)
         !! Integer array used for communication between the calling program and user routines.

      integer :: i, lplu, lju, ljlu, lrownrms, lperm, lqperm, ireorder, isrnorm, ipremeth, &
                 jscalcol

      ! Load IPREMETH, IREORDER and ISRNORM values from IPAR.
      ipremeth = ipar(5)
      ireorder = ipar(7)
      isrnorm = ipar(8)
      jscalcol = ipar(11)

      ! Load pointers into RWP and IWP arrays.
      lplu = ipar(21)
      lju = ipar(22)
      ljlu = ipar(23)
      lrownrms = ipar(24)
      lperm = ipar(25)
      lqperm = ipar(26)

      ! Scale c by multiplying by row-normalization factors, if used.
      if (isrnorm == 1) then
         do i = 1, neq
            b(i) = b(i)*rwp(lrownrms + i - 1)
         end do
      end if

      ! Solve P*x=b for a preconditioner stored as a sparse matrix in compressed sparse row
      ! format. If rows and columns of P were reordered (permuted), permute b, then use inverse
      ! permutation on x.
      if (ipremeth == 1 .or. ipremeth == 2) then
         if (ireorder == 1) call dvperm(neq, b, iwp(lperm))
         call lusol(neq, b, wk, rwp(lplu), iwp(ljlu), iwp(lju))
         if (ireorder == 1) call dvperm(neq, wk, iwp(lqperm))
      end if

      ! Unscale x by dividing by column scaling vector WGHT.
      if (jscalcol == 1) then
         b(1:neq) = wk(1:neq)/wght(1:neq)
      else
         b(1:neq) = wk(1:neq)
      end if

      ierr = 0

   end subroutine psol_ilupre