ImsLinear.f90

Go to the documentation of this file.

MODULE imslinearmodule
 
  use kindmodule, only: dp, i4b
  use constantsmodule, only: linelength, lensolutionname, lenmempath, &
                             izero, dzero, dprec, dsame, &
                             dem8, dem6, dem5, dem4, dem3, dem2, dem1, &
                             dhalf, done, dtwo, &
                             vdebug
  use imslinearbasemodule, only: ims_base_cg, ims_base_bcgs, &
                                 ims_base_pccrs, ims_base_calc_order, &
                                 ims_base_scale, ims_base_pcu, &
                                 ims_base_residual, ims_base_epfact, &
                                 ims_calc_pcdims
  use blockparsermodule, only: blockparsertype
  use matrixbasemodule
  use convergencesummarymodule
  use imslinearsettingsmodule
 
  IMPLICIT NONE
  private
 
  TYPE, PUBLIC :: imslineardatatype
    character(len=LENMEMPATH) :: memorypath !< the path for storing variables in the memory manager
    integer(I4B), POINTER :: iout => null() !< simulation listing file unit
    integer(I4B), POINTER :: iprims => null() !< print flag
    ! input variables (pointing to fields in input structure)
    real(dp), pointer :: dvclose => null() !< dependent variable closure criterion
    real(dp), pointer :: rclose => null() !< residual closure criterion
    integer(I4B), pointer :: icnvgopt => null() !< convergence option
    integer(I4B), pointer :: iter1 => null() !< max. iterations
    integer(I4B), pointer :: ilinmeth => null() !< linear solver method
    integer(I4B), pointer :: iscl => null() !< scaling method
    integer(I4B), pointer :: iord => null() !< reordering method
    integer(I4B), pointer :: north => null() !< number of orthogonalizations
 
    real(dp), pointer :: relax => null() !< relaxation factor
    integer(I4B), pointer :: level => null() !< nr. of preconditioner levels
    real(dp), pointer :: droptol => null() !< drop tolerance for preconditioner
    !
    integer(I4B), POINTER :: ipc => null() !< preconditioner flag
    integer(I4B), POINTER :: iacpc => null() !< preconditioner CRS row pointers
    integer(I4B), POINTER :: niterc => null() !<
    integer(I4B), POINTER :: niabcgs => null() !< size of working vectors for BCGS linear accelerator
    integer(I4B), POINTER :: niapc => null() !< preconditioner number of rows
    integer(I4B), POINTER :: njapc => null() !< preconditioner number of non-zero entries
    real(dp), POINTER :: epfact => null() !< factor for decreasing convergence criteria in seubsequent Picard iterations
    real(dp), POINTER :: l2norm0 => null() !< initial L2 norm
    ! -- ilut variables
    integer(I4B), POINTER :: njlu => null() !< length of jlu work vector
    integer(I4B), POINTER :: njw => null() !< length of jw work vector
    integer(I4B), POINTER :: nwlu => null() !< length of wlu work vector
    ! -- pointers to solution variables
    integer(I4B), POINTER :: neq => null() !< number of equations (rows in matrix)
    integer(I4B), POINTER :: nja => null() !< number of non-zero values in amat
    integer(I4B), dimension(:), pointer, contiguous :: ia => null() !< position of start of each row
    integer(I4B), dimension(:), pointer, contiguous :: ja => null() !< column pointer
    real(dp), dimension(:), pointer, contiguous :: amat => null() !< coefficient matrix
    real(dp), dimension(:), pointer, contiguous :: rhs => null() !< right-hand side of equation
    real(dp), dimension(:), pointer, contiguous :: x => null() !< dependent variable
    ! VECTORS
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: dscale => null() !< scaling factor
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: dscale2 => null() !< unscaling factor
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: iapc => null() !< position of start of each row in preconditioner matrix
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: japc => null() !< preconditioner matrix column pointer
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: apc => null() !< preconditioner coefficient matrix
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: lorder => null() !< reordering mapping
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: iorder => null() !< mapping to restore reordered matrix
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: iaro => null() !< position of start of each row in reordered matrix
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: jaro => null() !< reordered matrix column pointer
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: aro => null() !< reordered coefficient matrix
    ! WORKING ARRAYS
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: iw => null() !< integer working array
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: w => null() !< real working array
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: id => null() !< integer working array
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: d => null() !< real working array
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: p => null() !< real working array
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: q => null() !< real working array
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: z => null() !< real working array
    ! BICGSTAB WORKING ARRAYS
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: t => null() !< BICGSTAB real working array
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: v => null() !< BICGSTAB real working array
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: dhat => null() !< BICGSTAB real working array
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: phat => null() !< BICGSTAB real working array
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: qhat => null() !< rBICGSTAB eal working array
    ! POINTERS FOR USE WITH BOTH ORIGINAL AND RCM ORDERINGS
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: ia0 => null() !< pointer to current CRS row pointers
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: ja0 => null() !< pointer to current CRS column pointers
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: a0 => null() !< pointer to current coefficient matrix
    ! ILUT WORKING ARRAYS
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: jlu => null() !< ilut integer working array
    integer(I4B), POINTER, DIMENSION(:), CONTIGUOUS :: jw => null() !< ilut integer working array
    real(dp), POINTER, DIMENSION(:), CONTIGUOUS :: wlu => null() !< ilut real working array
 
    ! PROCEDURES (METHODS)
  CONTAINS
    PROCEDURE :: imslinear_allocate => imslinear_ar
    procedure :: imslinear_summary
    PROCEDURE :: imslinear_apply => imslinear_ap
    procedure :: imslinear_da => imslinear_da
    procedure, private :: allocate_scalars
    ! -- PRIVATE PROCEDURES
    PROCEDURE, PRIVATE :: set_imslinear_input => imslinear_set_input
  END TYPE imslineardatatype
 
CONTAINS
 
  !> @ brief Allocate storage and read data
    !!
    !!  Allocate storage for linear accelerators and read data
    !!
  !<
  SUBROUTINE imslinear_ar(this, NAME, IOUT, IPRIMS, MXITER, &
                          NEQ, matrix, RHS, X, linear_settings)
    ! -- modules
    use memorymanagermodule, only: mem_allocate
    use memoryhelpermodule, only: create_mem_path
    use simmodule, only: store_error, count_errors, &
                         deprecation_warning
    ! -- dummy variables
    CLASS(imslineardatatype), INTENT(INOUT) :: this !< ImsLinearDataType instance
    CHARACTER(LEN=LENSOLUTIONNAME), INTENT(IN) :: NAME !< solution name
    integer(I4B), INTENT(IN) :: IOUT !< simulation listing file unit
    integer(I4B), TARGET, INTENT(IN) :: IPRIMS !< print option
    integer(I4B), INTENT(IN) :: MXITER !< maximum outer iterations
    integer(I4B), TARGET, INTENT(IN) :: NEQ !< number of equations
    class(matrixbasetype), pointer :: matrix
    real(DP), DIMENSION(NEQ), TARGET, INTENT(INOUT) :: RHS !< right-hand side
    real(DP), DIMENSION(NEQ), TARGET, INTENT(INOUT) :: X !< dependent variables
    type(imslinearsettingstype), pointer :: linear_settings !< the settings form the IMS file
    ! -- local variables
    character(len=LINELENGTH) :: errmsg
    integer(I4B) :: n
    integer(I4B) :: i0
    integer(I4B) :: iscllen, iolen
 
    !
    ! -- DEFINE NAME
    this%memoryPath = create_mem_path(name, 'IMSLINEAR')
    !
    ! -- SET pointers to IMS settings
    this%DVCLOSE => linear_settings%dvclose
    this%RCLOSE => linear_settings%rclose
    this%ICNVGOPT => linear_settings%icnvgopt
    this%ITER1 => linear_settings%iter1
    this%ILINMETH => linear_settings%ilinmeth
    this%ISCL => linear_settings%iscl
    this%IORD => linear_settings%iord
    this%NORTH => linear_settings%north
    this%RELAX => linear_settings%relax
    this%LEVEL => linear_settings%level
    this%DROPTOL => linear_settings%droptol
    !
    ! -- SET POINTERS TO SOLUTION STORAGE
    this%IPRIMS => iprims
    this%NEQ => neq
    call matrix%get_aij(this%IA, this%JA, this%AMAT)
    call mem_allocate(this%NJA, 'NJA', this%memoryPath)
    this%NJA = size(this%AMAT)
    this%RHS => rhs
    this%X => x
    !
    ! -- ALLOCATE SCALAR VARIABLES
    call this%allocate_scalars()
    !
    ! -- initialize iout
    this%iout = iout
    !
    ! -- DEFAULT VALUES
    this%IPC = ipc_unknown
    !
    this%IACPC = 0
    !
    ! -- PRINT A MESSAGE IDENTIFYING IMSLINEAR SOLVER PACKAGE
    write (iout, 2000)
02000 FORMAT(1x, /1x, 'IMSLINEAR -- UNSTRUCTURED LINEAR SOLUTION', &
           ' PACKAGE, VERSION 8, 04/28/2017')
    !
    ! -- DETERMINE PRECONDITIONER
    this%IPC = resolve_ipc(this%LEVEL, this%RELAX)
    !
    ! -- ERROR CHECKING FOR OPTIONS
    IF (this%ISCL < 0) this%ISCL = 0
    IF (this%ISCL > 2) THEN
      WRITE (errmsg, '(A)') 'IMSLINEAR7AR ISCL MUST BE <= 2'
      call store_error(errmsg)
    END IF
    IF (this%IORD < 0) this%IORD = 0
    IF (this%IORD > 2) THEN
      WRITE (errmsg, '(A)') 'IMSLINEAR7AR IORD MUST BE <= 2'
      call store_error(errmsg)
    END IF
    IF (this%NORTH < 0) THEN
      WRITE (errmsg, '(A)') 'IMSLINEAR7AR NORTH MUST >= 0'
      call store_error(errmsg)
    END IF
    IF (this%RCLOSE == dzero) THEN
      IF (this%ICNVGOPT /= 3) THEN
        WRITE (errmsg, '(A)') 'IMSLINEAR7AR RCLOSE MUST > 0.0'
        call store_error(errmsg)
      END IF
    END IF
    IF (this%RELAX < dzero) THEN
      WRITE (errmsg, '(A)') 'IMSLINEAR7AR RELAX MUST BE >= 0.0'
      call store_error(errmsg)
    END IF
    IF (this%RELAX > done) THEN
      WRITE (errmsg, '(A)') 'IMSLINEAR7AR RELAX MUST BE <= 1.0'
      call store_error(errmsg)
    END IF
    !
    ! -- INITIALIZE IMSLINEAR VARIABLES
    this%NITERC = 0
    !
    ! -- ALLOCATE AND INITIALIZE MEMORY FOR IMSLINEAR
    iscllen = 1
    IF (this%ISCL .NE. 0) iscllen = neq
    CALL mem_allocate(this%DSCALE, iscllen, 'DSCALE', trim(this%memoryPath))
    CALL mem_allocate(this%DSCALE2, iscllen, 'DSCALE2', trim(this%memoryPath))
    !
    ! -- allocate and initialize the preconditioner work arrays
    call precond_allocate(this)
    !
    ! -- ALLOCATE SPACE FOR PERMUTATION VECTOR
    i0 = 1
    iolen = 1
    IF (this%IORD .NE. 0) THEN
      i0 = this%NEQ
      iolen = this%NJA
    END IF
    CALL mem_allocate(this%LORDER, i0, 'LORDER', trim(this%memoryPath))
    CALL mem_allocate(this%IORDER, i0, 'IORDER', trim(this%memoryPath))
    CALL mem_allocate(this%IARO, i0 + 1, 'IARO', trim(this%memoryPath))
    CALL mem_allocate(this%JARO, iolen, 'JARO', trim(this%memoryPath))
    CALL mem_allocate(this%ARO, iolen, 'ARO', trim(this%memoryPath))
    !
    ! -- ALLOCATE WORKING VECTORS FOR IMSLINEAR SOLVER
    CALL mem_allocate(this%ID, this%NEQ, 'ID', trim(this%memoryPath))
    CALL mem_allocate(this%D, this%NEQ, 'D', trim(this%memoryPath))
    CALL mem_allocate(this%P, this%NEQ, 'P', trim(this%memoryPath))
    CALL mem_allocate(this%Q, this%NEQ, 'Q', trim(this%memoryPath))
    CALL mem_allocate(this%Z, this%NEQ, 'Z', trim(this%memoryPath))
    !
    ! -- ALLOCATE MEMORY FOR BCGS WORKING ARRAYS
    this%NIABCGS = 1
    IF (this%ILINMETH == 2) THEN
      this%NIABCGS = this%NEQ
    END IF
    CALL mem_allocate(this%T, this%NIABCGS, 'T', trim(this%memoryPath))
    CALL mem_allocate(this%V, this%NIABCGS, 'V', trim(this%memoryPath))
    CALL mem_allocate(this%DHAT, this%NIABCGS, 'DHAT', trim(this%memoryPath))
    CALL mem_allocate(this%PHAT, this%NIABCGS, 'PHAT', trim(this%memoryPath))
    CALL mem_allocate(this%QHAT, this%NIABCGS, 'QHAT', trim(this%memoryPath))
    !
    ! -- INITIALIZE IMSLINEAR VECTORS
    DO n = 1, iscllen
      this%DSCALE(n) = done
      this%DSCALE2(n) = done
    END DO
    !
    ! -- WORKING VECTORS
    DO n = 1, this%NEQ
      this%ID(n) = izero
      this%D(n) = dzero
      this%P(n) = dzero
      this%Q(n) = dzero
      this%Z(n) = dzero
    END DO
    DO n = 1, this%NIAPC
      this%IW(n) = izero
      this%W(n) = dzero
    END DO
    !
    ! -- BCGS WORKING VECTORS
    DO n = 1, this%NIABCGS
      this%T(n) = dzero
      this%V(n) = dzero
      this%DHAT(n) = dzero
      this%PHAT(n) = dzero
      this%QHAT(n) = dzero
    END DO
    !
    ! -- ILUT AND MILUT WORKING VECTORS
    DO n = 1, this%NJLU
      this%JLU(n) = izero
    END DO
    DO n = 1, this%NJW
      this%JW(n) = izero
    END DO
    DO n = 1, this%NWLU
      this%WLU(n) = dzero
    END DO
    !
    ! -- REORDERING VECTORS
    DO n = 1, i0 + 1
      this%IARO(n) = izero
    END DO
    DO n = 1, iolen
      this%JARO(n) = izero
      this%ARO(n) = dzero
    END DO
    !
    ! -- REVERSE CUTHILL MCKEE AND MINIMUM DEGREE ORDERING
    IF (this%IORD .NE. 0) THEN
      CALL ims_base_calc_order(this%IORD, this%NEQ, this%NJA, this%IA, &
                               this%JA, this%LORDER, this%IORDER)
    END IF
    !
    ! -- ALLOCATE MEMORY FOR STORING ITERATION CONVERGENCE DATA
  end SUBROUTINE imslinear_ar
 
  !> @ brief Write summary of settings
    !!
    !!  Write summary of linear accelerator settings.
    !!
  !<
  subroutine imslinear_summary(this, mxiter)
    ! -- dummy variables
    class(imslineardatatype), intent(inout) :: this !< ImsLinearDataType instance
    integer(I4B), intent(in) :: mxiter !< maximum number of outer iterations
    ! -- local variables
    CHARACTER(LEN=10) :: clin(0:2)
    CHARACTER(LEN=31) :: clintit(0:2)
    CHARACTER(LEN=20) :: cipc(0:4)
    CHARACTER(LEN=20) :: cscale(0:2)
    CHARACTER(LEN=25) :: corder(0:2)
    CHARACTER(LEN=16), DIMENSION(0:4) :: ccnvgopt
    CHARACTER(LEN=15) :: clevel
    CHARACTER(LEN=15) :: cdroptol
    integer(I4B) :: i
    integer(I4B) :: j
    ! -- data
    DATA clin/'UNKNOWN   ', &
             &'CG        ', &
             &'BCGS      '/
    DATA clintit/'             UNKNOWN           ', &
                &'       CONJUGATE-GRADIENT      ', &
                &'BICONJUGATE-GRADIENT STABILIZED'/
    DATA cipc/'UNKNOWN             ', &
             &'INCOMPLETE LU       ', &
             &'MOD. INCOMPLETE LU  ', &
             &'INCOMPLETE LUT      ', &
             &'MOD. INCOMPLETE LUT '/
    DATA cscale/'NO SCALING          ', &
               &'SYMMETRIC SCALING   ', &
               &'L2 NORM SCALING     '/
    DATA corder/'ORIGINAL ORDERING        ', &
               &'RCM ORDERING             ', &
               &'MINIMUM DEGREE ORDERING  '/
    DATA ccnvgopt/'INFINITY NORM   ', &
                 &'INFINITY NORM S ', &
                 &'L2 NORM         ', &
                 &'RELATIVE L2NORM ', &
                 &'L2 NORM W. REL. '/
    ! -- formats
02010 FORMAT(1x, /, 7x, 'SOLUTION BY THE', 1x, a31, 1x, 'METHOD', &
           /, 1x, 66('-'), /, &
           ' MAXIMUM OF ', i0, ' CALLS OF SOLUTION ROUTINE', /, &
           ' MAXIMUM OF ', i0, &
           ' INTERNAL ITERATIONS PER CALL TO SOLUTION ROUTINE', /, &
           ' LINEAR ACCELERATION METHOD            =', 1x, a, /, &
           ' MATRIX PRECONDITIONING TYPE           =', 1x, a, /, &
           ' MATRIX SCALING APPROACH               =', 1x, a, /, &
           ' MATRIX REORDERING APPROACH            =', 1x, a, /, &
           ' NUMBER OF ORTHOGONALIZATIONS          =', 1x, i0, /, &
           ' HEAD CHANGE CRITERION FOR CLOSURE     =', e15.5, /, &
           ' RESIDUAL CHANGE CRITERION FOR CLOSURE =', e15.5, /, &
           ' RESIDUAL CONVERGENCE OPTION           =', 1x, i0, /, &
           ' RESIDUAL CONVERGENCE NORM             =', 1x, a, /, &
           ' RELAXATION FACTOR                     =', e15.5)
02015 FORMAT(' NUMBER OF LEVELS                      =', a15, /, &
           ' DROP TOLERANCE                        =', a15, //)
2030 FORMAT(1x, a20, 1x, 6(i6, 1x))
2040 FORMAT(1x, 20('-'), 1x, 6(6('-'), 1x))
2050 FORMAT(1x, 62('-'),/) !
! -- -----------------------------------------------------------
    !
    ! -- initialize clevel and cdroptol
    clevel = ''
    cdroptol = ''
    !
    ! -- write common variables to all linear accelerators
    write (this%iout, 2010) &
      clintit(this%ILINMETH), mxiter, this%ITER1, &
      clin(this%ILINMETH), cipc(this%IPC), &
      cscale(this%ISCL), corder(this%IORD), &
      this%NORTH, this%DVCLOSE, this%RCLOSE, &
      this%ICNVGOPT, ccnvgopt(this%ICNVGOPT), &
      this%RELAX
    if (this%level > 0) then
      write (clevel, '(i15)') this%level
    end if
    if (this%droptol > dzero) then
      write (cdroptol, '(e15.5)') this%droptol
    end if
    IF (this%level > 0 .or. this%droptol > dzero) THEN
      write (this%iout, 2015) trim(adjustl(clevel)), &
        trim(adjustl(cdroptol))
    ELSE
      write (this%iout, '(//)')
    END IF
 
    if (this%iord /= 0) then
      !
      ! -- WRITE SUMMARY OF REORDERING INFORMATION TO LIST FILE
      if (this%iprims == 2) then
        DO i = 1, this%neq, 6
          write (this%iout, 2030) 'ORIGINAL NODE      :', &
            (j, j=i, min(i + 5, this%neq))
          write (this%iout, 2040)
          write (this%iout, 2030) 'REORDERED INDEX    :', &
            (this%lorder(j), j=i, min(i + 5, this%neq))
          write (this%iout, 2030) 'REORDERED NODE     :', &
            (this%iorder(j), j=i, min(i + 5, this%neq))
          write (this%iout, 2050)
        END DO
      END IF
    end if
  end subroutine imslinear_summary
 
  !> @ brief Allocate and initialize scalars
    !!
    !!  Allocate and initialize linear accelerator scalars
    !!
  !<
  subroutine allocate_scalars(this)
    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy variables
    class(imslineardatatype), intent(inout) :: this !< ImsLinearDataType instance
    !
    ! -- allocate scalars
    call mem_allocate(this%iout, 'IOUT', this%memoryPath)
    call mem_allocate(this%ipc, 'IPC', this%memoryPath)
    call mem_allocate(this%iacpc, 'IACPC', this%memoryPath)
    call mem_allocate(this%niterc, 'NITERC', this%memoryPath)
    call mem_allocate(this%niabcgs, 'NIABCGS', this%memoryPath)
    call mem_allocate(this%niapc, 'NIAPC', this%memoryPath)
    call mem_allocate(this%njapc, 'NJAPC', this%memoryPath)
    call mem_allocate(this%epfact, 'EPFACT', this%memoryPath)
    call mem_allocate(this%l2norm0, 'L2NORM0', this%memoryPath)
    call mem_allocate(this%njlu, 'NJLU', this%memoryPath)
    call mem_allocate(this%njw, 'NJW', this%memoryPath)
    call mem_allocate(this%nwlu, 'NWLU', this%memoryPath)
    !
    ! -- initialize scalars
    this%iout = 0
    this%ipc = 0
    this%iacpc = 0
    this%niterc = 0
    this%niabcgs = 0
    this%niapc = 0
    this%njapc = 0
    this%epfact = dzero
    this%l2norm0 = 0
    this%njlu = 0
    this%njw = 0
    this%nwlu = 0
  end subroutine allocate_scalars
 
  !> @ brief Deallocate memory
    !!
    !!  Deallocate linear accelerator memory.
    !!
  !<
  subroutine imslinear_da(this)
    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy variables
    class(imslineardatatype), intent(inout) :: this !< linear datatype instance
    !
    ! -- arrays
    call mem_deallocate(this%dscale)
    call mem_deallocate(this%dscale2)
    call precond_destroy(this)
    call mem_deallocate(this%lorder)
    call mem_deallocate(this%iorder)
    call mem_deallocate(this%iaro)
    call mem_deallocate(this%jaro)
    call mem_deallocate(this%aro)
    call mem_deallocate(this%id)
    call mem_deallocate(this%d)
    call mem_deallocate(this%p)
    call mem_deallocate(this%q)
    call mem_deallocate(this%z)
    call mem_deallocate(this%t)
    call mem_deallocate(this%v)
    call mem_deallocate(this%dhat)
    call mem_deallocate(this%phat)
    call mem_deallocate(this%qhat)
    !
    ! -- scalars
    call mem_deallocate(this%iout)
    call mem_deallocate(this%ipc)
    call mem_deallocate(this%iacpc)
    call mem_deallocate(this%niterc)
    call mem_deallocate(this%niabcgs)
    call mem_deallocate(this%niapc)
    call mem_deallocate(this%njapc)
    call mem_deallocate(this%epfact)
    call mem_deallocate(this%l2norm0)
    call mem_deallocate(this%njlu)
    call mem_deallocate(this%njw)
    call mem_deallocate(this%nwlu)
    call mem_deallocate(this%NJA)
    !
    ! -- nullify pointers
    nullify (this%iprims)
    nullify (this%neq)
    nullify (this%nja)
    nullify (this%ia)
    nullify (this%ja)
    nullify (this%amat)
    nullify (this%rhs)
    nullify (this%x)
  end subroutine imslinear_da
 
  !> @brief Allocate and initialize the preconditioner work arrays
  !!
  !! Determines the preconditioner array dimensions from the current
  !! NEQ/NJA/LEVEL/IPC and allocates the ILU-family work arrays
  !! (IAPC/JAPC/APC/IW/W/JLU/JW/WLU); for ILU0/MILU0 the IAPC/JAPC sparsity is
  !! generated. Separated from imslinear_ar so the preconditioner can be
  !! (re)allocated at runtime.
  !<
  subroutine precond_allocate(this)
    use memorymanagermodule, only: mem_allocate
    ! -- dummy variables
    class(imslineardatatype), intent(inout) :: this !< ImsLinearDataType instance
    ! -- local variables
    integer(I4B) :: n
    !
    ! -- release any existing preconditioner arrays so this routine is
    !    idempotent and can be called again to reallocate at runtime
    if (associated(this%IAPC)) call precond_destroy(this)
    !
    ! -- determine dimensions for preconditioning arrays
    call ims_calc_pcdims(this%NEQ, this%NJA, this%IA, this%LEVEL, this%IPC, &
                         this%NIAPC, this%NJAPC, this%NJLU, this%NJW, this%NWLU)
    !
    ! -- allocate base preconditioner vectors
    call mem_allocate(this%IAPC, this%NIAPC + 1, 'IAPC', trim(this%memoryPath))
    call mem_allocate(this%JAPC, this%NJAPC, 'JAPC', trim(this%memoryPath))
    call mem_allocate(this%APC, this%NJAPC, 'APC', trim(this%memoryPath))
    !
    ! -- allocate ILU0/MILU0 non-zero row entry vectors
    call mem_allocate(this%IW, this%NIAPC, 'IW', trim(this%memoryPath))
    call mem_allocate(this%W, this%NIAPC, 'W', trim(this%memoryPath))
    !
    ! -- allocate ILUT vectors
    call mem_allocate(this%JLU, this%NJLU, 'JLU', trim(this%memoryPath))
    call mem_allocate(this%JW, this%NJW, 'JW', trim(this%memoryPath))
    call mem_allocate(this%WLU, this%NWLU, 'WLU', trim(this%memoryPath))
    !
    ! -- generate IAPC and JAPC for ILU0 and MILU0
    if (this%IPC == ipc_ilu0 .or. this%IPC == ipc_milu0) then
      call ims_base_pccrs(this%NEQ, this%NJA, this%IA, this%JA, &
                          this%IAPC, this%JAPC)
    end if
    !
    ! -- initialize preconditioner values
    do n = 1, this%NJAPC
      this%APC(n) = dzero
    end do
  end subroutine precond_allocate
 
  !> @brief Deallocate the preconditioner work arrays
  !!
  !! Preconditioner-only teardown, separated from imslinear_da so the
  !! preconditioner can be destroyed and reallocated at runtime without
  !! disturbing the solver work arrays.
  !<
  subroutine precond_destroy(this)
    use memorymanagerextmodule, only: memorystore_release
    ! -- dummy variables
    class(imslineardatatype), intent(inout) :: this !< ImsLinearDataType instance
    !
    ! -- mem_deallocate is a no-op in the memory manager, so releasing the named
    !    store entries (and nullifying the aliases below) is what actually frees
    !    these arrays and lets the preconditioner be reallocated at runtime
    call memorystore_release('IAPC', trim(this%memoryPath))
    call memorystore_release('JAPC', trim(this%memoryPath))
    call memorystore_release('APC', trim(this%memoryPath))
    call memorystore_release('IW', trim(this%memoryPath))
    call memorystore_release('W', trim(this%memoryPath))
    call memorystore_release('JLU', trim(this%memoryPath))
    call memorystore_release('JW', trim(this%memoryPath))
    call memorystore_release('WLU', trim(this%memoryPath))
    !
    nullify (this%IAPC)
    nullify (this%JAPC)
    nullify (this%APC)
    nullify (this%IW)
    nullify (this%W)
    nullify (this%JLU)
    nullify (this%JW)
    nullify (this%WLU)
  end subroutine precond_destroy
 
  !> @ brief Set default settings
    !!
    !!  Set default linear accelerator settings.
    !!
  !<
  SUBROUTINE imslinear_set_input(this, IFDPARAM)
    ! -- dummy variables
    CLASS(imslineardatatype), INTENT(INOUT) :: this !< ImsLinearDataType instance
    integer(I4B), INTENT(IN) :: IFDPARAM !< complexity option
    ! -- code
    SELECT CASE (ifdparam)
      !
      ! -- Simple option
    CASE (1)
      this%ITER1 = 50
      this%ILINMETH = 1
      this%IPC = ipc_ilu0
      this%ISCL = 0
      this%IORD = 0
      this%DVCLOSE = dem3
      this%RCLOSE = dem1
      this%RELAX = dzero
      this%LEVEL = 0
      this%DROPTOL = dzero
      this%NORTH = 0
      !
      ! -- Moderate
    CASE (2)
      this%ITER1 = 100
      this%ILINMETH = 2
      this%IPC = ipc_milu0
      this%ISCL = 0
      this%IORD = 0
      this%DVCLOSE = dem2
      this%RCLOSE = dem1
      this%RELAX = 0.97d0
      this%LEVEL = 0
      this%DROPTOL = dzero
      this%NORTH = 0
      !
      ! -- Complex
    CASE (3)
      this%ITER1 = 500
      this%ILINMETH = 2
      this%IPC = ipc_ilut
      this%ISCL = 0
      this%IORD = 0
      this%DVCLOSE = dem1
      this%RCLOSE = dem1
      this%RELAX = dzero
      this%LEVEL = 5
      this%DROPTOL = dem4
      this%NORTH = 2
    END SELECT
  end SUBROUTINE imslinear_set_input
 
  !> @ brief Base linear accelerator subroutine
    !!
    !!  Base linear accelerator subroutine that scales and reorders
    !!  the system of equations, if necessary, updates the preconditioner,
    !!  and calls the appropriate linear accelerator.
    !!
  !<
  SUBROUTINE imslinear_ap(this, ICNVG, KSTP, KITER, IN_ITER, &
                          NCONV, CONVNMOD, CONVMODSTART, &
                          CACCEL, summary)
    ! -- modules
    USE simmodule
    ! -- dummy variables
    CLASS(imslineardatatype), INTENT(INOUT) :: this !< ImsLinearDataType instance
    integer(I4B), INTENT(INOUT) :: ICNVG !< convergence flag (1) non-convergence (0)
    integer(I4B), INTENT(IN) :: KSTP !< time step number
    integer(I4B), INTENT(IN) :: KITER !< outer iteration number
    integer(I4B), INTENT(INOUT) :: IN_ITER !< inner iteration number
    ! -- convergence information dummy variables
    integer(I4B), INTENT(IN) :: NCONV !<
    integer(I4B), INTENT(IN) :: CONVNMOD !<
    integer(I4B), DIMENSION(CONVNMOD + 1), INTENT(INOUT) :: CONVMODSTART !<
    character(len=31), DIMENSION(NCONV), INTENT(INOUT) :: CACCEL !<
    type(convergencesummarytype), pointer, intent(in) :: summary !< Convergence summary report
    ! -- local variables
    integer(I4B) :: n
    integer(I4B) :: innerit
    integer(I4B) :: irc
    integer(I4B) :: itmax
    real(DP) :: dnrm2
    !
    ! -- set epfact based on timestep
    this%EPFACT = ims_base_epfact(this%ICNVGOPT, kstp)
    !
    ! -- SCALE PROBLEM
    IF (this%ISCL .NE. 0) THEN
      CALL ims_base_scale(0, this%ISCL, &
                          this%NEQ, this%NJA, this%IA, this%JA, &
                          this%AMAT, this%X, this%RHS, &
                          this%DSCALE, this%DSCALE2)
    END IF
    !
    ! -- PERMUTE ROWS, COLUMNS, AND RHS
    IF (this%IORD /= 0) THEN
      CALL dperm(this%NEQ, this%AMAT, this%JA, this%IA, &
                 this%ARO, this%JARO, this%IARO, &
                 this%LORDER, this%ID, 1)
      CALL dvperm(this%NEQ, this%X, this%LORDER)
      CALL dvperm(this%NEQ, this%RHS, this%LORDER)
      this%IA0 => this%IARO
      this%JA0 => this%JARO
      this%A0 => this%ARO
    ELSE
      this%IA0 => this%IA
      this%JA0 => this%JA
      this%A0 => this%AMAT
    END IF
    !
    ! -- UPDATE PRECONDITIONER
    CALL ims_base_pcu(this%iout, this%NJA, this%NEQ, this%NIAPC, this%NJAPC, &
                      this%IPC, this%RELAX, this%A0, this%IA0, this%JA0, &
                      this%APC, this%IAPC, this%JAPC, this%IW, this%W, &
                      this%LEVEL, this%DROPTOL, this%NJLU, this%NJW, &
                      this%NWLU, this%JLU, this%JW, this%WLU)
    !
    ! -- INITIALIZE SOLUTION VARIABLE AND ARRAYS
    IF (kiter == 1) then
      this%NITERC = 0
      summary%iter_cnt = 0
    end if
    irc = 1
    icnvg = 0
    DO n = 1, this%NEQ
      this%D(n) = dzero
      this%P(n) = dzero
      this%Q(n) = dzero
      this%Z(n) = dzero
    END DO
    !
    ! -- CALCULATE INITIAL RESIDUAL
    call ims_base_residual(this%NEQ, this%NJA, this%X, this%RHS, this%D, &
                           this%A0, this%IA0, this%JA0)
    this%L2NORM0 = dnrm2(this%NEQ, this%D, 1)
    !
    ! -- CHECK FOR EXACT SOLUTION
    itmax = this%ITER1
    IF (this%L2NORM0 == dzero) THEN
      itmax = 0
      icnvg = 1
    END IF
    !
    ! -- SOLUTION BY THE CONJUGATE GRADIENT METHOD
    IF (this%ILINMETH == 1) THEN
      CALL ims_base_cg(icnvg, itmax, innerit, &
                       this%NEQ, this%NJA, this%NIAPC, this%NJAPC, &
                       this%IPC, this%ICNVGOPT, this%NORTH, &
                       this%DVCLOSE, this%RCLOSE, this%L2NORM0, &
                       this%EPFACT, this%IA0, this%JA0, this%A0, &
                       this%IAPC, this%JAPC, this%APC, &
                       this%X, this%RHS, this%D, this%P, this%Q, this%Z, &
                       this%NJLU, this%IW, this%JLU, &
                       nconv, convnmod, convmodstart, &
                       caccel, summary)
      !
      ! -- SOLUTION BY THE BICONJUGATE GRADIENT STABILIZED METHOD
    ELSE IF (this%ILINMETH == 2) THEN
      CALL ims_base_bcgs(icnvg, itmax, innerit, &
                         this%NEQ, this%NJA, this%NIAPC, this%NJAPC, &
                         this%IPC, this%ICNVGOPT, this%NORTH, &
                         this%ISCL, this%DSCALE, &
                         this%DVCLOSE, this%RCLOSE, this%L2NORM0, &
                         this%EPFACT, this%IA0, this%JA0, this%A0, &
                         this%IAPC, this%JAPC, this%APC, &
                         this%X, this%RHS, this%D, this%P, this%Q, &
                         this%T, this%V, this%DHAT, this%PHAT, this%QHAT, &
                         this%NJLU, this%IW, this%JLU, &
                         nconv, convnmod, convmodstart, &
                         caccel, summary)
    END IF
    !
    ! -- BACK PERMUTE AMAT, SOLUTION, AND RHS
    IF (this%IORD /= 0) THEN
      CALL dperm(this%NEQ, this%A0, this%JA0, this%IA0, &
                 this%AMAT, this%JA, this%IA, &
                 this%IORDER, this%ID, 1)
      CALL dvperm(this%NEQ, this%X, this%IORDER)
      CALL dvperm(this%NEQ, this%RHS, this%IORDER)
    END IF
    !
    ! -- UNSCALE PROBLEM
    IF (this%ISCL .NE. 0) THEN
      CALL ims_base_scale(1, this%ISCL, &
                          this%NEQ, this%NJA, this%IA, this%JA, &
                          this%AMAT, this%X, this%RHS, &
                          this%DSCALE, this%DSCALE2)
    END IF
    !
    ! -- SET IMS INNER ITERATION NUMBER (IN_ITER) TO NUMBER OF
    !       IMSLINEAR INNER ITERATIONS (innerit)
    in_iter = innerit
  end SUBROUTINE imslinear_ap
 
END MODULE imslinearmodule