jac_ilupre Subroutine

   subroutine jac_ilupre( &
      res, ires, neq, t, y, ydot, rewt, savr, wk, h, cj, rwp, iwp, ierr, rpar, ipar)
   !! This subroutine uses finite-differences to calculate the Jacobian matrix in sparse format,
   !! and then performs an incomplete LU decomposition using either [[ilut]] or [[ilutp]] from
   !! SPARSKIT.

      use dsparskit, only: amudia, dvperm, prtmt, roscal

      procedure(res_t) :: res
         !! User-defined residuals routine.
      integer, intent(out) :: ires
         !! Error flag set by `res`.
      integer, intent(in) :: neq
         !! Problem size.
      real(rk), intent(in) :: t
         !! Current independent variable.
      real(rk), intent(inout) :: y(neq)
         !! Current dependent variables.
      real(rk), intent(inout) :: ydot(neq)
         !! Current derivatives of dependent variables.
      real(rk), intent(in) :: rewt(neq)
         !! Reciprocal error weights for scaling `y` and `ydot`.
      real(rk), intent(inout) :: savr(neq)
         !! Current residual evaluated at `(t, y, ydot)`.
      real(rk), intent(inout) :: wk(neq)
         !! Real work space available to this subroutine.
      real(rk), intent(in) :: h
         !! Current step size.
      real(rk), intent(in) :: cj
         !! Scalar used in forming the system Jacobian.
      real(rk), intent(inout) :: rwp(*)
         !! Matrix elements of ILU.
      integer, intent(inout) :: iwp(*)
         !! Array indices for elements of ILU.
      integer, intent(inout) :: ierr
         !! Error flag (0 means success, else failure).
      real(rk), intent(inout) :: rpar(*)
         !! Real array used for communication between the calling program and user routines.
      integer, intent(inout) :: ipar(*)
         !! Integer array used for communication between the calling program and user routines.

      external ::  xerrwd !@todo: replace by module

      character(len=8), parameter :: PMETH(4) = [character(len=8) :: 'ILUT', 'ILUTP', 'ILU0', 'MILU0']
      real(rk) :: tolilut, permtol, sqrtn
      integer :: i, lbw, ubw, lenplumx, ljac, ljaci, ljacj, liperm, lrownrms, lrwk1, liwk1, &
                 ifmt, lenpfac, lenplufac, lfililut, ipremeth, neqp1, nnzmx, lplu, lju, ljlu, &
                 lperm, lqperm, llevels, lmask
      integer :: isrnorm  ! =1 causes row normalization of JAC.
      integer :: normtype ! =0,1,2 for max-norm, 1-norm, or 2-norm row scaling
      integer :: ireorder ! =1 causes row and column reordering of JAC.
      integer :: jacout   ! =1 causes the Jacobian matrix and SAVR to be written to a file and
                          ! then exit with IRES = -2 to signal a stop to DDASPK.
      integer :: iunit    ! logical unit number to use when JACOUT == 1
      integer :: jscalcol ! =1 causes the columns of the Jacobian matrix to be scaled by REWT-inverse
      integer :: nre      ! number of RES calls needed to evaluate Jacobian. NRE is returned in IPAR(30).
      character(len=8) :: premeth
      character(len=72) :: title
      character(len=80) :: msg

      ! Zero out NRE counter (number of RES calls needed to evaluate Jacobian)
      nre = 0

      ! Load values from IPAR and RPAR.
      lbw = ipar(1)
      ubw = ipar(2)
      lenpfac = ipar(3)
      lenplufac = ipar(4)
      ipremeth = ipar(5)
      lfililut = ipar(6)
      ireorder = ipar(7)
      isrnorm = ipar(8)
      normtype = ipar(9)
      jacout = ipar(10)
      jscalcol = ipar(11)
      tolilut = rpar(1)
      permtol = rpar(2)
      premeth = pmeth(ipremeth)

      ! Set pointers into the WP and IWP arrays.
      neqp1 = neq + 1
      nnzmx = lenpfac*neq
      lenplumx = nnzmx + lenplufac*neq

      ! Set up pointers into WP
      ljac = 1
      lrownrms = nnzmx + ljac
      if (isrnorm == 1) then
         lplu = lrownrms + neq
      else
         lplu = lrownrms
      end if
      lrwk1 = lplu + lenplumx

      ! Set up pointers into IWP
      ljaci = 1
      ljacj = ljaci + neqp1
      lju = ljacj + nnzmx
      ljlu = lju + lenplumx

      ! Calculate Jacobian matrix.
      ierr = 0
      call jcalc(neq, t, y, ydot, savr, lbw, ubw, wk, rewt, res, h, cj, nnzmx, rwp(ljac), &
                 iwp(ljacj), iwp(ljaci), rwp(lplu), iwp(ljlu), iwp(lju), ipar, rpar, ires, &
                 nre, ierr)
      if (ires < 0) return
      if (ierr /= 0) return

      ! Save NRE value for user output.
      ipar(30) = ipar(30) + nre

      ! Modify pointers into IWP
      ljlu = lju + neqp1
      if (ireorder /= 0) then
         lperm = ljlu + lenplumx
         lqperm = lperm + neq
         liwk1 = lqperm + neq
         llevels = ljlu + nnzmx          ! assumes that LENPLUFAC >= 2.
         lmask = llevels + neq
      else
         lperm = 0
         lqperm = 0
         llevels = 0
         lmask = 0
         liwk1 = ljlu + lenplumx
      end if

      if (premeth == 'ILUTP') then
         liperm = liwk1 + 2*neq
      else
         liperm = liwk1
      end if

      ! Multiply Jacobian columns by inverse of scaling vector REWT.
      ! In PSOLILU, the WGHT array equals REWT/SQRT(NEQ), so we must be consistent here.
      if (jscalcol == 1) then
         sqrtn = sqrt(real(neq))
         do i = 1, neq
            wk(i) = sqrtn/rewt(i)
         end do
         call amudia(neq, 0, rwp(ljac), iwp(ljacj), iwp(ljaci), wk, rwp(ljac), iwp(ljacj), iwp(ljaci))
      end if

      ! Normalize Jacobian rows, if desired.
      if (isrnorm == 1) then
         call roscal(neq, 0, normtype, rwp(ljac), iwp(ljacj), iwp(ljaci), &
                     rwp(lrownrms), rwp(ljac), iwp(ljacj), iwp(ljaci), ierr)
         if (ierr /= 0) return
      end if

      ! Reorder Jacobian rows and columns, if desired.
      if (ireorder == 1) then
         call jreord(neq, nnzmx, &
                     rwp(ljac), iwp(ljacj), iwp(ljaci), &
                     rwp(lplu), iwp(ljlu), iwp(lju), &
                     iwp(lperm), iwp(lqperm), iwp(llevels), iwp(lmask))
      end if

      ! Write matrix JAC and scaled RES to file if JACOUT == 1.
      if (jacout == 1) then
         iunit = ipar(29)
         if (isrnorm == 1) then
            do i = 1, neq
               savr(i) = savr(i)*rwp(lrownrms + i - 1)
            end do
         end if
         if (ireorder /= 0) call dvperm(neq, savr, iwp(lperm))
         title = 'DASPK Test Matrix '
         ifmt = 15
         call prtmt(neq, neq, rwp(ljac), iwp(ljacj), iwp(ljaci), savr, &
                    'NN', title, 'SPARSKIT', 'RUA', ifmt, 3, iunit)
         msg = 'DJACILU -- Jacobian Matrix written to file.'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         ierr = 1
         ires = -2
         return
      end if

      ! Compute ILU decomposition.
      call jilu(neq, nnzmx, rwp(ljac), iwp(ljacj), iwp(ljaci), iwp(lju), rwp(lplu), iwp(ljlu), &
                rwp(lrwk1), iwp(liwk1), lenplumx, tolilut, &
                lfililut, permtol, premeth, iwp(liperm), ierr)
      if ((ierr == -2) .or. (ierr == -3)) then
         ires = -2 ! Stop since more storage needed.
      end if

      ! Save pointers for use in DPSOLILU into IPAR array.
      ipar(21) = lplu
      ipar(22) = lju
      ipar(23) = ljlu
      ipar(24) = lrownrms
      ipar(25) = lperm
      ipar(26) = lqperm

   end subroutine jac_ilupre