daskr_ilupre.f90 Source File


This file depends on

sourcefile~~daskr_ilupre.f90~~EfferentGraph sourcefile~daskr_ilupre.f90 daskr_ilupre.f90 sourcefile~daskr_kinds.f90 daskr_kinds.F90 sourcefile~daskr_ilupre.f90->sourcefile~daskr_kinds.f90 sourcefile~daskr_new.f90 daskr_new.f90 sourcefile~daskr_ilupre.f90->sourcefile~daskr_new.f90 sourcefile~dsparskit.f dsparskit.f sourcefile~daskr_ilupre.f90->sourcefile~dsparskit.f sourcefile~daskr_new.f90->sourcefile~daskr_kinds.f90 sourcefile~blas_interfaces.f90 blas_interfaces.f90 sourcefile~daskr_new.f90->sourcefile~blas_interfaces.f90 sourcefile~linpack.f90 linpack.f90 sourcefile~daskr_new.f90->sourcefile~linpack.f90 sourcefile~blas_interfaces.f90->sourcefile~daskr_kinds.f90 sourcefile~linpack.f90->sourcefile~daskr_kinds.f90 sourcefile~linpack.f90->sourcefile~blas_interfaces.f90

Files dependent on this one

sourcefile~~daskr_ilupre.f90~~AfferentGraph sourcefile~daskr_ilupre.f90 daskr_ilupre.f90 sourcefile~example_heatilu.f90 example_heatilu.f90 sourcefile~example_heatilu.f90->sourcefile~daskr_ilupre.f90

Source Code

!----------------------------------------------------------------------------------------------
! Adapted from original Fortran code in `original/preconds/dilupre.f`
!----------------------------------------------------------------------------------------------

module daskr_ilupre
!! # Preconditioner Routines for Sparse Problems
!!
!! This module provides a general-purpose sparse incomplete LU (ILU) preconditioner for use with 
!! the [[daskr]] solver, with the Krylov linear system method. 
!!
!! When using [[daskr]] to solve a problem \(G(t,y,\dot{y}) = 0\), whose Jacobian
!! \( J = \partial G/ \partial y + c_J \partial G/ \partial \dot{y} \), where \(c_J\) is a
!! scalar, is a general sparse matrix, the routines [[jac_ilupre]] and [[psol_ilupre]] can be
!! used to generate an approximation to \(J\) as the preconditioner and to solve the resulting
!! sparse linear system, in conjunction with the Krylov method option.
!!
!! Internally, the incomplete LU factorization is achieved via one of two routines — [[ilut]]
!! or [[ilutp]] — from the SPARSKIT library. The routine [[ilut]] performs an ILU factorization
!! of a sparse matrix using a dual thresholding technique based on a drop tolerance (`tolilut`)
!! and a level of fill-in parameter (`lfililut`). The parameter `lfililut` controls the amount
!! of fill-in allowed in the factorization (limited to a maximum of `2*lfililut*neq`, but normally
!! much less). Increasing `lfililut` will generally make the ILU factorization more accurate.
!! The parameter `tolilut` also controls the accuracy of the ILU factorization via a drop tolerance
!! based on element size. Decreasing `tolilut` will increase the amount of fill-in and make for
!! a more accurate factorization. The routine [[ilutp]] is a variant of [[ilut]] that in addition
!! performs pivoting based on a tolerance ratio `permtol`.
!!
!! An important aspect of using incomplete factorization techniques is that of reordering the
!! rows and columns in the Jacobian matrix \(J\) before performing the ILU. In this package,
!! this is accomplished via the parameter `ireorder`, which when equal to 1 performs a reverse
!! Cuthill-McKee (RCM) reordering before performing the ILU factorization. Based on the limited
!! amount of testing done so far, RCM seems the best overall choice. It is possible to include
!! a different reordering technique if desired.
!!
!! ## Usage
!!  
!! To use these routines in conjunction with [[daskr]], the user's calling program should include
!! the following, in addition to setting the other [[daskr]] input parameters:
!
!! * Dimension the array `ipar` to have length at least 30, and load the following parameters
!!   into `ipar` as:
!!   
!!    | Index  | Name         | Description                                                      |
!!    |--------|--------------|------------------------------------------------------------------|
!!    | 1      | `ml`         | The lower bandwidth used in calculating \(J\).                   |
!!    | 2      | `mu`         | The upper bandwidth used in calculating \(J\).                   |
!!    | 3      | `lenpfac`    | The average number of nonzeros in a row of \(J\). The maximum of nonzeros allowed in \(J\) is `nnzmx = lenpfac*neq`. `lenpfac >= 2`. |
!!    | 4      | `lenplufac`  | The average amount of fill-in per row in the factored \(J\). The maximum number of nonzeros allowed in the factored \(J\) is `lenplumx = nnzmx + lenplufac*neq`. `lenplufac >=2`. |
!!    | 5      | `ipremeth`   | Preconditioner type flag. `=1` means [[ilut]] and `=2` means [[ilutp]]. |
!!    | 6      | `lfililut`   | Fill-in parameter for [[ilut]] and [[ilutp]]. The largest `lfililut` elements per row of the L and U factors are kept. Each row of L and U will have a maximum of `lfililut` elements in addition to their original number of nonzero elements. |
!!    | 7      | `ireorder`   | Reordering flag. `=0` means no reordering of \(J\) rows/columns before incomplete factorization. `=1` means reverse Cuthill-McKee (RCM) reordering is performed. |
!!    | 8      | `isrnorm`    | Row norm flag. `=1` means compute and use row norms as scalings in the preconditioner system \(P x = b\). `=0` means no row norm scaling. |
!!    | 9      | `normtype`   | Type of row norm scaling for `isrnorm`. `=0` means max-norm. `=1` means 1-norm. `=2` means 2-norm. |
!!    | 10     | `jacout`     | Output Jacobian flag. `=1` means write \(J\) and initial residual \(G\) to a file (logical unit in `ipar(29)`). Integration halts with `ires = -2`. `=0` means no output. Storage format is Boeing-Harwell. |
!!    | 11     | `jscalcol`   | Flag for scaling \(J\) columns by the inverses of elements in the `ewt` array. `=0` means no scaling. `=1` means perform scaling. |
!!    | 21:28  | —            | Used to hold pointer information.                                |
!!    | 29     | `jacout_unit`| Logical unit number for matrix output file. Used only if `jacout = 1`. |
!!    | 30     | `rescalls`   | On return from [[daskr]], holds the number of calls to the `res` routine used in preconditioner evaluations. | 
!!
!! * Dimension the array `rpar` to have length at least 2, and load the parameters shown in the
!!   table below into `rpar`.
!!  
!!    | Index | Name      | Description                                                          |
!!    |-------|-----------|----------------------------------------------------------------------|
!!    | 1     | `tolilut` | Drop tolerance for use by [[ilut]] and [[ilutp]]. `tolilut >= 0`. Larger values cause less fill-in. Good values range from 0.001 to 0.01.                                   |
!!    | 2     | `permtol` | Tolerance ratio used in determining column pivoting by ILUTP. `permtol >= 0`. Good values are from 0.1 to 0.01. Two columns are permuted only if `a(i,j)*permtol > a(i,i)`. |
!!  
!! * The two parameters `tolilut` and `lfililut` give the user a great deal of flexibility. One 
!!   can use `tolilut = 0` to get a strategy based on keeping the largest elements in each row 
!!   of \(L\) and \(U\). Taking `tolilut /= 0` but `lfililut = neq` will give the usual threshold
!!   strategy (however, fill-in is then unpredictable).
!!
!! * Set `info(12) = 1` to select the Krylov iterative method and `info(15) = 1` to indicate
!!   that a `jac` routine exists. Then in the call to [[daskr]], pass the procedure names
!!   `jac_ilupre` and `psol_ilupre` as the arguments `jac` and `psol`, respectively.
!!
!! * The [[daskr]] work arrays `rwork` and `iwork` must include segments `rwp` and `iwp` for use
!!   by [[jac_ilupre]] and [[psol_ilupre]]. The lengths of these depend on the problem size, 
!!   half-bandwidths, and other parameters as shown in the table below. Load these lengths in
!!   `iwork` as `iwork(27) = lrwp` and `iwork(28) = liwp` and include these values in the declared
!!   size of `rwork` and `iwork`, respectively.
!!
!!    | Variable | Length                                                                              |
!!    |----------|-------------------------------------------------------------------------------------|
!!    | `lrwp`   | `2*lenpfac*neq + lenplufac*neq + isrnorm*neq + neq`                                 |
!!    | `liwp`   | `3*neq + 1 + 3*lenpfac*neq + 2*lenplufac*neq + 2*ireorder*neq + 2*(ipremeth-1)*neq` |
!!      
!! ## Example
!!
!! The program [[example_heatilu]] demonstrates the use of this preconditioner.

   use daskr_kinds, only: rk, zero, one
   use daskr, only: res_t
   implicit none
   private

   public :: setup_ilupre, jac_ilupre, psol_ilupre

contains

   pure subroutine setup_ilupre(neq, lrwp, liwp, rpar, ipar, ierr, lrwp_min, liwp_min)
   !! Setup routine for the incomplete LU preconditioner. This routine checks the user input
   !! and calculates the minimum length needed for the preconditioner workspace arrays.

      integer, intent(in) :: neq
         !! Problem size.
      integer, intent(in) :: lrwp
         !! Current length of `rwp`.
      integer, intent(in) :: liwp
         !! Current length of `iwp`.
      real(rk), intent(in) :: rpar(*)
         !! User real workspace.
      integer, intent(in) :: ipar(*)
         !! User integer workspace.
      integer, intent(out) :: ierr
         !! Error flag (0 means success, else failure):
         !! `1 <= ierr <= 11` means there's an illegal value for `ipar(ierr)`;
         !! `ierr = 12` means `ipar(29)` is illegal;
         !! `ierr = 21` means `rpar(1)` is illegal;
         !! `ierr = 22` means `rpar(2)` is illegal;
         !! `ierr = 30` means more `wp` length is needed;
         !! `ierr = 31` means more `iwp` length is needed.
      integer, intent(out) :: lrwp_min
         !! Minimum `rwp` length needed.
      integer, intent(out) :: liwp_min
         !! Minimum `iwp` length needed.

      integer :: lbw, ubw, lenplumx, ljac, ljaci, ljacj, lrownrms, lrwk1, liwk1, 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 ierr = 1 to signal a stop to daskr.
      integer :: jscalcol ! =1 causes the columns of the Jacobian matrix to be scaled by EWT-inverse
      real(rk)  :: tolilut, permtol

      ! Load values from IPAR and RPAR. Check for illegal values.
      lbw = ipar(1)        ! LBW must be > 0
      if (lbw <= 0) then
         ierr = 1
         return
      end if

      ubw = ipar(2)        ! UBW must be > 0
      if (ubw <= 0) then
         ierr = 2
         return
      end if

      lenpfac = ipar(3)    ! LENPFAC must be >= 2
      if (lenpfac <= 1) then
         ierr = 3
         return
      end if

      lenplufac = ipar(4)  ! LENPLUFAC must be >= 2
      if (lenplufac <= 1) then
         ierr = 4
         return
      end if

      ipremeth = ipar(5)   ! IPREMETH must be == 1 or 2 currently
      if (ipremeth /= 1 .and. ipremeth /= 2) then
         ierr = 5
         return
      end if

      lfililut = ipar(6)   ! LFILILUT must be >= 0
      if (lfililut < 0) then
         ierr = 6
         return
      end if

      ireorder = ipar(7)   ! IREORDER must be 0 or 1
      if ((ireorder < 0) .or. (ireorder > 1)) then
         ierr = 7
         return
      end if

      isrnorm = ipar(8)    ! ISRNORM must be 0 or 1
      if ((isrnorm < 0) .or. (isrnorm > 1)) then
         ierr = 8
         return
      end if

      normtype = ipar(9)   ! NORMTYPE must be 0, 1, or 2
      if ((normtype < 0) .or. (normtype > 2)) then
         ierr = 9
         return
      end if

      jacout = ipar(10)    ! JACOUT must be 0 or 1
      if ((jacout < 0) .or. (jacout > 1)) then
         ierr = 10
         return
      end if

      jscalcol = ipar(11)    ! JSCALCOL must be 0 or 1
      if ((jscalcol < 0) .or. (jscalcol > 1)) then
         ierr = 11
         return
      end if

      if (jacout == 1) then  ! IPAR(29) must be > 0
         if (ipar(29) <= 0) then
            ierr = 12
            return
         end if
      end if

      tolilut = rpar(1)    ! TOLILUT must be >= 0.0
      if (tolilut < zero) then
         ierr = 21
         return
      end if

      if (ipremeth == 2) then
         permtol = rpar(2)        ! PERMTOL must be >= 0.0
         if (permtol < zero) then
            ierr = 22
            return
         end if
      end if

      ! Calculate minimum work lengths for 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
      lrwp_min = lrwk1 + neq - 1
      if (lrwp < lrwp_min) then
         ierr = 30  ! more WP length needed.
         return
      end if

      ! Set up pointers into IWP
      ljaci = 1
      ljacj = ljaci + neqp1
      lju = ljacj + nnzmx
      ljlu = lju + max(lenplumx, 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

      liwp_min = liwk1 + 2*neq - 1
      if (ipremeth == 2) liwp_min = liwp_min + 2*neq
      if (liwp < liwp_min) then
         ierr = 31 ! more IWP length needed.
         return
      end if

      ierr = 0

   end subroutine setup_ilupre

   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

   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

   subroutine jcalc( &
      neq, t, y, ydot, r0, ml, mu, r1, rewt, res, h, cj, nnzmx, jac, ja, ia, rcoo, jcoo, &
      icoo, ipar, rpar, ires, nre, ierr)
   !! This subroutine calculates the Jacobian matrix by one-sided finite-differences. Lower and
   !! upper bandwidths are used to select the elements to be computed. The Jacobian is stored
   !! in compressed sparse row format.

      use dsparskit, only: coocsr

      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) :: r0(neq)
         !! Current residual evaluated at `(t, y, ydot)`.
      integer, intent(in) :: ml
         !! Lower bandwidth.
      integer, intent(in) :: mu
         !! Upper bandwidth.
      real(rk), intent(inout) :: r1(neq)
         !! Real work space available to this subroutine.
      real(rk), intent(in) :: rewt(neq)
         !! Reciprocal error weights for scaling `y` and `ydot`.
      procedure(res_t) :: res
         !! User-defined residuals routine.
      real(rk), intent(in) ::  h
         !! Current step size.
      real(rk), intent(in) :: cj
         !! Scalar used in forming the system Jacobian.
      integer, intent(in) :: nnzmx
         !! Maximum number of nonzeros in Jacobian.
      real(rk), intent(out) :: jac(nnzmx)
         !! Nonzero Jacobian elements.
      integer, intent(out) :: ja(nnzmx)
         !! Column indices of nonzero Jacobian elements.
      integer, intent(out) :: ia(neq + 1)
         !! Pointers to beginning of each row in `jac` and `ja`.
      real(rk), intent(inout) :: rcoo(nnzmx)
         !! Nonzero Jacobian elements.
      integer, intent(out) :: jcoo(nnzmx)
         !! Column indices of nonzero Jacobian elements.
      integer, intent(out) :: icoo(nnzmx)
         !! Row indices of nonzero Jacobian elements.
      integer, intent(inout) :: ipar(*)
         !! User integer workspace.
      real(rk), intent(inout) :: rpar(*)
         !! User real workspace.
      integer, intent(out) :: ires
         !! Error flag for `res` routine.
      integer, intent(inout) :: nre
         !! Number of calls to `res`.
      integer, intent(out) :: ierr
         !! Error flag.

      external :: xerrwd

      integer :: nnz, i, i1, i2, j, jj, mba, meband, meb1, mband
      real(rk) :: jacelem, squround, del, delinv
      character(len=80) :: msg

      ! Set band parameters.
      nnz = 1
      mband = ml + mu + 1
      mba = min(mband, neq)
      meband = mband + ml
      meb1 = meband - 1

      ! Set the machine unit roundoff UROUND and SQRT(UROUND), used to set increments in the
      ! difference quotient procedure.
      squround = sqrt(epsilon(one))

      ! Initial error flags.
      ierr = 0
      ires = 0

      ! Make MBA calls to RES to approximate the Jacobian.
      ! Here, R0(1),...,R0(neq) contains the base RES value, and
      ! R1(1),...,R1(NEQ) contains the perturbed values of RES.
      do j = 1, mba

         do jj = j, neq, mband
            jac(jj) = y(jj)
            jac(jj + neq) = ydot(jj)
            del = squround*max(abs(y(jj)), abs(h*ydot(jj)), abs(one/rewt(jj)))
            del = sign(del, h*ydot(jj))
            del = (y(jj) + del) - y(jj)
            y(jj) = y(jj) + del
            ydot(jj) = ydot(jj) + cj*del
         end do

         call res(t, y, ydot, cj, r1, ires, rpar, ipar)
         if (ires < 0) return

         nre = nre + 1
         do jj = j, neq, mband
            y(jj) = jac(jj)
            ydot(jj) = jac(jj + neq)
            del = squround*max(abs(y(jj)), abs(h*ydot(jj)), abs(one/rewt(jj)))
            del = sign(del, h*ydot(jj))
            del = (y(jj) + del) - y(jj)
            delinv = one/del
            i1 = max(1, (jj - mu))
            i2 = min(neq, (jj + ml))
            do i = i1, i2
               ! Calculate possibly nonzero Jacobian elements for this variable,
               ! and store nonzero elements in coordinate format.
               jacelem = (r1(i) - r0(i))*delinv
               if (jacelem /= zero) then
                  if (nnz > nnzmx) then
                     msg = 'DJCALC -- More storage needed for Jacobian.'
                     call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
                     msg = 'DJCALC -- Increase LENPFAC.'
                     call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
                     msg = 'DJCALC -- Storage exceeded at (I,J) = (I1,I2)'
                     call xerrwd(msg, 80, 0, 0, 2, i, jj, 0, 0.0, 0.0)
                     ierr = 1
                     ires = -2
                     return
                  end if
                  rcoo(nnz) = jacelem
                  jcoo(nnz) = jj
                  icoo(nnz) = i
                  nnz = nnz + 1
               end if
            end do
         end do
      end do
      nnz = nnz - 1

      ! Convert Jacobian from coordinate to compressed sparse row format.
      call coocsr(neq, nnz, rcoo, icoo, jcoo, jac, ja, ia)

   end subroutine jcalc

   subroutine jilu( &
      neq, nnzmx, jac, ja, ia, ju, plu, jlu, rwk1, iwk1, lenplumx, tolilut, &
      lfililut, permtol, premeth, iperm, ierr)
   !! This subroutine computes the incomplete LU decomposition of the Jacobian matrix and returns
   !! it in any given storage format.

      use dsparskit, only: ilut, ilutp
      
      integer, intent(in) :: neq
         !! Problem size.
      integer, intent(in) :: nnzmx
         !! Maximum number of nonzeros in Jacobian.
      real(rk), intent(in) :: jac(nnzmx)
         !! Nonzero Jacobian elements.
      integer, intent(in) :: ja(nnzmx)
         !! Column indices of nonzero Jacobian elements.
      integer, intent(in) :: ia(neq + 1)
         !! Pointers to beginning of each row in `jac` and `ja`.
      integer, intent(out) :: ju(neq)
         !! Pointer to beginning of each row of U in matrix `plu` and `jlu`.
      real(rk), intent(out) :: plu(lenplumx)
         !! Matrix elements of ILU.
      integer, intent(out) :: jlu(lenplumx)
         !! Sizes and array indices for elements of ILU.
      real(rk), intent(inout) :: rwk1(neq)
         !! Real work space.
      integer, intent(inout) :: iwk1(2*neq)
         !! Integer work space.
      integer, intent(in) :: lenplumx
         !! Length of `plu` and `jlu` arrays.
      real(rk), intent(in) :: tolilut
         !! Tolerance for ILU decomposition.
      integer, intent(in) :: lfililut
         !! @todo: find description
      real(rk), intent(in) :: permtol
         !! Tolerance for permutation.
      character(len=8), intent(in) :: premeth
         !! Preconditioner method.
      integer, intent(inout) :: iperm(2*neq)
         !! Integer work space.
      integer, intent(out) :: ierr
         !! Error flag (0 means success, else failure).

      external :: xerrwd ! @todo: replace by module

      character(len=80) :: msg
      logical :: err

      err = .false.

      if (premeth == 'ILUT') then
         ! Use incomplete factorization routine ILUT from SparsKit.
         call ilut(neq, jac, ja, ia, lfililut, tolilut, plu, jlu, ju, lenplumx, rwk1, iwk1, ierr)
         if (ierr /= 0) then
            msg = 'DJILU -- Error return from ILUT: IERR = (I1)'
            call xerrwd(msg, 80, 0, 0, 1, ierr, 0, 0, 0.0, 0.0)
            err = .true.
         end if
      elseif (premeth == 'ILUTP') then
         ! Use incomplete factorization routine ILUTP from SparsKit.
         call ilutp(neq, jac, ja, ia, lfililut, tolilut, permtol, neq, plu, jlu, ju, lenplumx, &
                    rwk1, iwk1, iperm, ierr)
         if (ierr /= 0) then
            msg = 'DJILU -- Error return from ILUTP: IERR = (I1)'
            call xerrwd(msg, 80, 0, 0, 1, ierr, 0, 0, 0.0, 0.0)
            err = .true.
         end if
      end if

      ! Put in other options here for incomplete factorizations.
      ! @todo: refactor this
      if (err) then
         msg = 'DJILU -- IERR /= 0 means one of the following has occurred:'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '    IERR >  0   --> Zero pivot encountered at step number IERR.'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '    IERR = -1   --> Error. input matrix may be wrong.'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '                     (The elimination process has generated a'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '                     row in L or U with length > NEQ.)'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '    IERR = -2   --> Matrix L overflows.'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '    IERR = -3   --> Matrix U overflows.'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '    IERR = -4   --> Illegal value for LFILILUT.'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '    IERR = -5   --> Zero row encountered.'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '    '
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '    For IERR = -2 or -3, increase the value of LENPLUFAC or'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '    decrease the value of LFILILUT if LENPLUFAC cannot be'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
         msg = '    increased.'
         call xerrwd(msg, 80, 0, 0, 0, 0, 0, 0, 0.0, 0.0)
      end if

   end subroutine jilu

   subroutine jreord(neq, nnzmx, jac, ja, ia, awk, jwk, iwk, perm, qperm, levels, mask)
   !! This subroutine reorders the Jacobian matrix.

      use dsparskit, only: atob, bfs, rversp, dperm

      integer, intent(in) :: neq
         !! Problem size.
      integer, intent(in) :: nnzmx
         !! Maximum number of nonzeroes in Jacobian.
      real(rk), intent(in) :: jac(nnzmx)
         !! Nonzero Jacobian elements.
      integer, intent(in) :: ja(nnzmx)
         !! Column indices of nonzero Jacobian elements.
      integer, intent(in) :: ia(neq + 1)
         !! Indices of 1st nonzero element in each row.
      real(rk), intent(inout) :: awk(nnzmx)
         !! Work array.
      integer, intent(inout) :: jwk(nnzmx)
         !! Work array.
      integer, intent(inout) :: iwk(neq + 1)
         !! Work array.
      integer, intent(inout) :: perm(neq)
         !! Array containing the permutation used in reordering the rows and columns of the
         !! Jacobian matrix.
      integer, intent(inout) :: qperm(neq)
         !! Integer array holding the inverse of the permutation in array `perm`.
      integer, intent(inout) :: levels(neq)
         !! Work array used by the [[bfs]] reordering subroutine.
      integer, intent(inout) :: mask(neq)
         !! Work array used by the [[bfs]] reordering subroutine.

      integer :: i, nfirst, nlev, maskval

      ! Copy JAC, JA, and IA to AWK, JWK, and IWK.
      call atob(neq, jac, ja, ia, awk, jwk, iwk)

      ! Perform a Cuthill-McKee reordering of the Jacobian.
      nfirst = 1
      perm(1) = 0
      mask = 1
      maskval = 1
      qperm(1) = 1
      call bfs(neq, jwk, iwk, nfirst, perm, mask, maskval, qperm, levels, nlev)

      ! Reverse the permutation to obtain the reverse Cuthill-McKee reordering.
      call rversp(neq, qperm)

      ! Calculate the inverse of QPERM and put it in PERM.
      do i = 1, neq
         perm(qperm(i)) = i
      end do

      ! Permute rows and columns of Jacobian using PERM.
      call dperm(neq, awk, jwk, iwk, jac, ja, ia, perm, perm, 1)

   end subroutine jreord

end module daskr_ilupre