setup_ilupre Subroutine

   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