MODFLOW 6  version 6.8.0.dev0
USGS Modular Hydrologic Model
Functions/Subroutines | Variables
imslinearbasemodule Module Reference

This module contains the IMS linear accelerator subroutines. More...

Functions/Subroutines

subroutine ims_base_cg (ICNVG, ITMAX, INNERIT, NEQ, NJA, NIAPC, NJAPC, IPC, ICNVGOPT, NORTH, DVCLOSE, RCLOSE, L2NORM0, EPFACT, IA0, JA0, A0, IAPC, JAPC, APC, X, B, D, P, Q, Z, NJLU, IW, JLU, NCONV, CONVNMOD, CONVMODSTART, CACCEL, summary)
 @ brief Preconditioned Conjugate Gradient linear accelerator More...
 
subroutine ims_base_bcgs (ICNVG, ITMAX, INNERIT, NEQ, NJA, NIAPC, NJAPC, IPC, ICNVGOPT, NORTH, ISCL, DSCALE, DVCLOSE, RCLOSE, L2NORM0, EPFACT, IA0, JA0, A0, IAPC, JAPC, APC, X, B, D, P, Q, T, V, DHAT, PHAT, QHAT, NJLU, IW, JLU, NCONV, CONVNMOD, CONVMODSTART, CACCEL, summary)
 @ brief Preconditioned BiConjugate Gradient Stabilized linear accelerator More...
 
subroutine ims_base_calc_order (IORD, NEQ, NJA, IA, JA, LORDER, IORDER)
 @ brief Calculate LORDER AND IORDER More...
 
subroutine ims_base_scale (IOPT, ISCL, NEQ, NJA, IA, JA, AMAT, X, B, DSCALE, DSCALE2)
 @ brief Scale the coefficient matrix More...
 
subroutine ims_base_pcu (IOUT, NJA, NEQ, NIAPC, NJAPC, IPC, RELAX, AMAT, IA, JA, APC, IAPC, JAPC, IW, W, LEVEL, DROPTOL, NJLU, NJW, NWLU, JLU, JW, WLU)
 @ brief Update the preconditioner More...
 
subroutine ims_base_pcjac (NJA, NEQ, AMAT, APC, IA, JA)
 @ brief Jacobi preconditioner More...
 
subroutine ims_base_jaca (NEQ, A, D1, D2)
 @ brief Apply the Jacobi preconditioner More...
 
subroutine ims_base_pcilu0 (NJA, NEQ, AMAT, IA, JA, APC, IAPC, JAPC, IW, W, RELAX, IPCFLAG, DELTA)
 @ brief Update the ILU0 preconditioner More...
 
subroutine ims_base_ilu0a (NJA, NEQ, APC, IAPC, JAPC, R, D)
 @ brief Apply the ILU0 and MILU0 preconditioners More...
 
subroutine ims_base_testcnvg (Icnvgopt, Icnvg, Iiter, Dvmax, Rmax, Rmax0, Epfact, Dvclose, Rclose)
 @ brief Test for solver convergence More...
 
subroutine ims_calc_pcdims (neq, nja, ia, level, ipc, niapc, njapc, njlu, njw, nwlu)
 
subroutine ims_base_pccrs (NEQ, NJA, IA, JA, IAPC, JAPC)
 @ brief Generate CRS pointers for the preconditioner More...
 
subroutine ims_base_isort (NVAL, IARRAY)
 In-place sorting for an integer array. More...
 
subroutine ims_base_residual (NEQ, NJA, X, B, D, A, IA, JA)
 Calculate residual. More...
 
real(dp) function ims_base_epfact (icnvgopt, kstp)
 Function returning EPFACT. More...
 

Variables

type(blockparsertype), private parser
 

Detailed Description

This module contains the IMS linear accelerator subroutines used by a MODFLOW 6 solution.

Function/Subroutine Documentation

◆ ims_base_bcgs()

subroutine imslinearbasemodule::ims_base_bcgs ( integer(i4b), intent(inout)  ICNVG,
integer(i4b), intent(in)  ITMAX,
integer(i4b), intent(inout)  INNERIT,
integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NIAPC,
integer(i4b), intent(in)  NJAPC,
integer(i4b), intent(in)  IPC,
integer(i4b), intent(in)  ICNVGOPT,
integer(i4b), intent(in)  NORTH,
integer(i4b), intent(in)  ISCL,
real(dp), dimension(neq), intent(in)  DSCALE,
real(dp), intent(in)  DVCLOSE,
real(dp), intent(in)  RCLOSE,
real(dp), intent(in)  L2NORM0,
real(dp), intent(in)  EPFACT,
integer(i4b), dimension(neq + 1), intent(in)  IA0,
integer(i4b), dimension(nja), intent(in)  JA0,
real(dp), dimension(nja), intent(in)  A0,
integer(i4b), dimension(niapc + 1), intent(in)  IAPC,
integer(i4b), dimension(njapc), intent(in)  JAPC,
real(dp), dimension(njapc), intent(in)  APC,
real(dp), dimension(neq), intent(inout)  X,
real(dp), dimension(neq), intent(in)  B,
real(dp), dimension(neq), intent(inout)  D,
real(dp), dimension(neq), intent(inout)  P,
real(dp), dimension(neq), intent(inout)  Q,
real(dp), dimension(neq), intent(inout)  T,
real(dp), dimension(neq), intent(inout)  V,
real(dp), dimension(neq), intent(inout)  DHAT,
real(dp), dimension(neq), intent(inout)  PHAT,
real(dp), dimension(neq), intent(inout)  QHAT,
integer(i4b), intent(in)  NJLU,
integer(i4b), dimension(niapc), intent(in)  IW,
integer(i4b), dimension(njlu), intent(in)  JLU,
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), intent(in), pointer  summary 
)

Apply the Preconditioned BiConjugate Gradient Stabilized linear accelerator to the current coefficient matrix, right-hand side, using the currentdependent-variable.

Parameters
[in,out]icnvgconvergence flag (1) non-convergence (0)
[in]itmaxmaximum number of inner iterations
[in,out]inneritinner iteration count
[in]neqnumber of equations
[in]njanumber of non-zero entries
[in]niapcpreconditioner number of rows
[in]njapcpreconditioner number of non-zero entries
[in]ipcpreconditioner option
[in]icnvgoptflow convergence criteria option
[in]northorthogonalization frequency
[in]isclscaling option
[in]dscalescaling vector
[in]dvclosedependent-variable closure criteria
[in]rcloseflow closure criteria
[in]l2norm0initial L-2 norm for system of equations
[in]epfactfactor for decreasing flow convergence criteria for subsequent Picard iterations
[in]ia0CRS row pointers
[in]ja0CRS column pointers
[in]a0coefficient matrix
[in]iapcpreconditioner CRS row pointers
[in]japcpreconditioner CRS column pointers
[in]apcpreconditioner matrix
[in,out]xdependent-variable vector
[in]bright-hand side vector
[in,out]dpreconditioner working vector
[in,out]ppreconditioner working vector
[in,out]qpreconditioner working vector
[in,out]tpreconditioner working vector
[in,out]vpreconditioner working vector
[in,out]dhatBCGS preconditioner working vector
[in,out]phatBCGS preconditioner working vector
[in,out]qhatBCGS preconditioner working vector
[in]njlupreconditioner length of JLU vector
[in]iwpreconditioner integer working vector
[in]jlupreconditioner JLU working vector
[in]nconvmaximum number of inner iterations in a time step (maxiter * maxinner)
[in]convnmodnumber of models in the solution
[in,out]convmodstartpointer to the start of each model in the convmod* arrays
[in,out]caccelconvergence string
[in]summaryConvergence summary report

Definition at line 250 of file ImsLinearBase.f90.

260  ! -- dummy variables
261  integer(I4B), INTENT(INOUT) :: ICNVG !< convergence flag (1) non-convergence (0)
262  integer(I4B), INTENT(IN) :: ITMAX !< maximum number of inner iterations
263  integer(I4B), INTENT(INOUT) :: INNERIT !< inner iteration count
264  integer(I4B), INTENT(IN) :: NEQ !< number of equations
265  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
266  integer(I4B), INTENT(IN) :: NIAPC !< preconditioner number of rows
267  integer(I4B), INTENT(IN) :: NJAPC !< preconditioner number of non-zero entries
268  integer(I4B), INTENT(IN) :: IPC !< preconditioner option
269  integer(I4B), INTENT(IN) :: ICNVGOPT !< flow convergence criteria option
270  integer(I4B), INTENT(IN) :: NORTH !< orthogonalization frequency
271  integer(I4B), INTENT(IN) :: ISCL !< scaling option
272  real(DP), DIMENSION(NEQ), INTENT(IN) :: DSCALE !< scaling vector
273  real(DP), INTENT(IN) :: DVCLOSE !< dependent-variable closure criteria
274  real(DP), INTENT(IN) :: RCLOSE !< flow closure criteria
275  real(DP), INTENT(IN) :: L2NORM0 !< initial L-2 norm for system of equations
276  real(DP), INTENT(IN) :: EPFACT !< factor for decreasing flow convergence criteria for subsequent Picard iterations
277  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA0 !< CRS row pointers
278  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA0 !< CRS column pointers
279  real(DP), DIMENSION(NJA), INTENT(IN) :: A0 !< coefficient matrix
280  integer(I4B), DIMENSION(NIAPC + 1), INTENT(IN) :: IAPC !< preconditioner CRS row pointers
281  integer(I4B), DIMENSION(NJAPC), INTENT(IN) :: JAPC !< preconditioner CRS column pointers
282  real(DP), DIMENSION(NJAPC), INTENT(IN) :: APC !< preconditioner matrix
283  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: X !< dependent-variable vector
284  real(DP), DIMENSION(NEQ), INTENT(IN) :: B !< right-hand side vector
285  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: D !< preconditioner working vector
286  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: P !< preconditioner working vector
287  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: Q !< preconditioner working vector
288  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: T !< preconditioner working vector
289  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: V !< preconditioner working vector
290  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: DHAT !< BCGS preconditioner working vector
291  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: PHAT !< BCGS preconditioner working vector
292  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: QHAT !< BCGS preconditioner working vector
293  ! -- ILUT dummy variables
294  integer(I4B), INTENT(IN) :: NJLU !< preconditioner length of JLU vector
295  integer(I4B), DIMENSION(NIAPC), INTENT(IN) :: IW !< preconditioner integer working vector
296  integer(I4B), DIMENSION(NJLU), INTENT(IN) :: JLU !< preconditioner JLU working vector
297  ! -- convergence information dummy variables
298  integer(I4B), INTENT(IN) :: NCONV !< maximum number of inner iterations in a time step (maxiter * maxinner)
299  integer(I4B), INTENT(IN) :: CONVNMOD !< number of models in the solution
300  integer(I4B), DIMENSION(CONVNMOD + 1), INTENT(INOUT) :: CONVMODSTART !< pointer to the start of each model in the convmod* arrays
301  character(len=31), DIMENSION(NCONV), INTENT(INOUT) :: CACCEL !< convergence string
302  type(ConvergenceSummaryType), pointer, intent(in) :: summary !< Convergence summary report
303  ! -- local variables
304  LOGICAL :: LORTH
305  logical :: lsame
306  character(len=15) :: cval1, cval2
307  integer(I4B) :: n
308  integer(I4B) :: iiter
309  integer(I4B) :: xloc, rloc
310  integer(I4B) :: im, im0, im1
311  real(DP) :: ddot
312  real(DP) :: tv
313  real(DP) :: deltax
314  real(DP) :: rmax
315  real(DP) :: l2norm
316  real(DP) :: rcnvg
317  real(DP) :: alpha, alpha0
318  real(DP) :: beta
319  real(DP) :: rho, rho0
320  real(DP) :: omega, omega0
321  real(DP) :: numerator, denominator
322  !
323  ! -- initialize local variables
324  innerit = 0
325  alpha = dzero
326  alpha0 = dzero
327  beta = dzero
328  rho = dzero
329  rho0 = dzero
330  omega = dzero
331  omega0 = dzero
332  !
333  ! -- SAVE INITIAL RESIDUAL
334  DO n = 1, neq
335  dhat(n) = d(n)
336  END DO
337  !
338  ! -- INNER ITERATION
339  inner: DO iiter = 1, itmax
340  innerit = innerit + 1
341  summary%iter_cnt = summary%iter_cnt + 1
342  !
343  ! -- CALCULATE rho
344  rho = ddot(neq, dhat, 1, d, 1)
345  !
346  ! -- COMPUTE DIRECTIONAL VECTORS
347  IF (iiter == 1) THEN
348  DO n = 1, neq
349  p(n) = d(n)
350  END DO
351  ELSE
352  beta = (rho / rho0) * (alpha0 / omega0)
353  DO n = 1, neq
354  p(n) = d(n) + beta * (p(n) - omega0 * v(n))
355  END DO
356  END IF
357  !
358  ! -- APPLY PRECONDITIONER TO UPDATE PHAT
359  SELECT CASE (ipc)
360  !
361  ! -- ILU0 AND MILU0
362  CASE (ipc_ilu0, ipc_milu0)
363  CALL ims_base_ilu0a(nja, neq, apc, iapc, japc, p, phat)
364  !
365  ! -- ILUT AND MILUT
366  CASE (ipc_ilut, ipc_milut)
367  CALL lusol(neq, p, phat, apc, jlu, iw)
368  END SELECT
369  !
370  ! -- COMPUTE ITERATES
371  !
372  ! -- UPDATE V WITH A AND PHAT
373  call amux(neq, phat, v, a0, ja0, ia0)
374  !
375  ! -- UPDATE alpha WITH DHAT AND V
376  denominator = ddot(neq, dhat, 1, v, 1)
377  denominator = denominator + sign(dprec, denominator)
378  alpha = rho / denominator
379  !
380  ! -- UPDATE Q
381  DO n = 1, neq
382  q(n) = d(n) - alpha * v(n)
383  END DO
384  !
385  ! ! -- CALCULATE INFINITY NORM OF Q - TEST FOR TERMINATION
386  ! ! TERMINATE IF rmax IS LESS THAN MACHINE PRECISION (DPREC)
387  ! rmax = DZERO
388  ! DO n = 1, NEQ
389  ! tv = Q(n)
390  ! IF (ISCL.NE.0 ) tv = tv / DSCALE(n)
391  ! IF (ABS(tv) > ABS(rmax) ) rmax = tv
392  ! END DO
393  ! IF (ABS(rmax).LE.DPREC) THEN
394  ! deltax = DZERO
395  ! DO n = 1, NEQ
396  ! tv = alpha * PHAT(n)
397  ! IF (ISCL.NE.0) THEN
398  ! tv = tv * DSCALE(n)
399  ! END IF
400  ! X(n) = X(n) + tv
401  ! IF (ABS(tv) > ABS(deltax) ) deltax = tv
402  ! END DO
403  ! CALL IMSLINEARSUB_TESTCNVG(ICNVGOPT, ICNVG, INNERIT, &
404  ! deltax, rmax, &
405  ! rmax, EPFACT, DVCLOSE, RCLOSE )
406  ! IF (ICNVG.NE.0 ) EXIT INNER
407  ! END IF
408  !
409  ! -- APPLY PRECONDITIONER TO UPDATE QHAT
410  SELECT CASE (ipc)
411  !
412  ! -- ILU0 AND MILU0
413  CASE (ipc_ilu0, ipc_milu0)
414  CALL ims_base_ilu0a(nja, neq, apc, iapc, japc, q, qhat)
415  !
416  ! -- ILUT AND MILUT
417  CASE (ipc_ilut, ipc_milut)
418  CALL lusol(neq, q, qhat, apc, jlu, iw)
419  END SELECT
420  !
421  ! -- UPDATE T WITH A AND QHAT
422  call amux(neq, qhat, t, a0, ja0, ia0)
423  !
424  ! -- UPDATE omega
425  numerator = ddot(neq, t, 1, q, 1)
426  denominator = ddot(neq, t, 1, t, 1)
427  denominator = denominator + sign(dprec, denominator)
428  omega = numerator / denominator
429  !
430  ! -- UPDATE X AND RESIDUAL
431  deltax = dzero
432  rmax = dzero
433  l2norm = dzero
434  DO im = 1, convnmod
435  summary%dvmax(im) = dzero
436  summary%rmax(im) = dzero
437  END DO
438  im = 1
439  im0 = convmodstart(1)
440  im1 = convmodstart(2)
441  DO n = 1, neq
442  !
443  ! -- determine current model index
444  if (n == im1) then
445  im = im + 1
446  im0 = convmodstart(im)
447  im1 = convmodstart(im + 1)
448  end if
449  !
450  ! -- X AND DX
451  tv = alpha * phat(n) + omega * qhat(n)
452  x(n) = x(n) + tv
453  IF (iscl .NE. 0) THEN
454  tv = tv * dscale(n)
455  END IF
456  IF (abs(tv) > abs(deltax)) THEN
457  deltax = tv
458  xloc = n
459  END IF
460  IF (abs(tv) > abs(summary%dvmax(im))) THEN
461  summary%dvmax(im) = tv
462  summary%locdv(im) = n
463  END IF
464  !
465  ! -- RESIDUAL
466  tv = q(n) - omega * t(n)
467  d(n) = tv
468  IF (iscl .NE. 0) THEN
469  tv = tv / dscale(n)
470  END IF
471  IF (abs(tv) > abs(rmax)) THEN
472  rmax = tv
473  rloc = n
474  END IF
475  IF (abs(tv) > abs(summary%rmax(im))) THEN
476  summary%rmax(im) = tv
477  summary%locr(im) = n
478  END IF
479  l2norm = l2norm + tv * tv
480  END DO
481  l2norm = sqrt(l2norm)
482  !
483  ! -- SAVE SOLVER convergence information dummy variables
484  IF (nconv > 1) THEN !<
485  n = summary%iter_cnt
486  WRITE (cval1, '(g15.7)') alpha
487  WRITE (cval2, '(g15.7)') omega
488  caccel(n) = trim(adjustl(cval1))//','//trim(adjustl(cval2))
489  summary%itinner(n) = iiter
490  DO im = 1, convnmod
491  summary%convdvmax(im, n) = summary%dvmax(im)
492  summary%convlocdv(im, n) = summary%locdv(im)
493  summary%convrmax(im, n) = summary%rmax(im)
494  summary%convlocr(im, n) = summary%locr(im)
495  END DO
496  END IF
497  !
498  ! -- TEST FOR SOLVER CONVERGENCE
499  IF (icnvgopt == 2 .OR. icnvgopt == 3 .OR. icnvgopt == 4) THEN
500  rcnvg = l2norm
501  ELSE
502  rcnvg = rmax
503  END IF
504  CALL ims_base_testcnvg(icnvgopt, icnvg, innerit, &
505  deltax, rcnvg, &
506  l2norm0, epfact, dvclose, rclose)
507  !
508  ! -- CHECK FOR EXACT SOLUTION
509  IF (rcnvg == dzero) icnvg = 1
510  !
511  ! -- CHECK FOR STANDARD CONVERGENCE
512  IF (icnvg .NE. 0) EXIT inner
513  !
514  ! -- CHECK THAT CURRENT AND PREVIOUS rho, alpha, AND omega ARE
515  ! DIFFERENT
516  lsame = is_close(rho, rho0)
517  IF (lsame) THEN
518  EXIT inner
519  END IF
520  lsame = is_close(alpha, alpha0)
521  IF (lsame) THEN
522  EXIT inner
523  END IF
524  lsame = is_close(omega, omega0)
525  IF (lsame) THEN
526  EXIT inner
527  END IF
528  !
529  ! -- RECALCULATE THE RESIDUAL
530  IF (north > 0) THEN
531  lorth = mod(iiter + 1, north) == 0
532  IF (lorth) THEN
533  call ims_base_residual(neq, nja, x, b, d, a0, ia0, ja0)
534  END IF
535  END IF
536  !
537  ! -- exit inner if rho or omega are zero
538  if (rho * omega == dzero) then
539  exit inner
540  end if
541  !
542  ! -- SAVE CURRENT INNER ITERATES
543  rho0 = rho
544  alpha0 = alpha
545  omega0 = omega
546  END DO inner
547  !
548  ! -- RESET ICNVG
549  IF (icnvg < 0) icnvg = 0
ddot
real(kind=8) function ddot(n, dx, incx, dy, incy)
Definition: blas1_d.f90:296
lusol
subroutine lusol(n, y, x, alu, jlu, ju)
Definition: ilut.f90:466
amux
subroutine amux(n, x, y, a, ja, ia)
Definition: sparsekit.f90:2
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_calc_order()

subroutine imslinearbasemodule::ims_base_calc_order ( integer(i4b), intent(in)  IORD,
integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA,
integer(i4b), dimension(neq), intent(inout)  LORDER,
integer(i4b), dimension(neq), intent(inout)  IORDER 
)

Calculate LORDER and IORDER for reordering.

Parameters
[in]iordreordering option
[in]neqnumber of rows
[in]njanumber of non-zero entries
[in]iarow pointer
[in]jacolumn pointer
[in,out]lorderreorder vector
[in,out]iorderinverse of reorder vector

Definition at line 557 of file ImsLinearBase.f90.

558  ! -- modules
559  use simmodule, only: store_error, count_errors
560  ! -- dummy variables
561  integer(I4B), INTENT(IN) :: IORD !< reordering option
562  integer(I4B), INTENT(IN) :: NEQ !< number of rows
563  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
564  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< row pointer
565  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< column pointer
566  integer(I4B), DIMENSION(NEQ), INTENT(INOUT) :: LORDER !< reorder vector
567  integer(I4B), DIMENSION(NEQ), INTENT(INOUT) :: IORDER !< inverse of reorder vector
568  ! -- local variables
569  character(len=LINELENGTH) :: errmsg
570  integer(I4B) :: n
571  integer(I4B) :: nsp
572  integer(I4B), DIMENSION(:), ALLOCATABLE :: iwork0
573  integer(I4B), DIMENSION(:), ALLOCATABLE :: iwork1
574  integer(I4B) :: iflag
575  !
576  ! -- initialize lorder and iorder
577  DO n = 1, neq
578  lorder(n) = izero
579  iorder(n) = izero
580  END DO
581  ! ALLOCATE (iwork0(NEQ))
582  SELECT CASE (iord)
583  CASE (1)
584  CALL genrcm(neq, nja, ia, ja, lorder)
585  CASE (2)
586  nsp = 3 * neq + 4 * nja
587  allocate (iwork0(neq))
588  allocate (iwork1(nsp))
589  CALL ims_odrv(neq, nja, nsp, ia, ja, lorder, iwork0, &
590  iwork1, iflag)
591  IF (iflag .NE. 0) THEN
592  write (errmsg, '(A,1X,A)') &
593  'IMSLINEARSUB_CALC_ORDER error creating minimum degree ', &
594  'order permutation '
595  call store_error(errmsg)
596  END IF
597  !
598  ! -- DEALLOCATE TEMPORARY STORAGE
599  deallocate (iwork0, iwork1)
600  END SELECT
601  !
602  ! -- GENERATE INVERSE OF LORDER
603  DO n = 1, neq
604  iorder(lorder(n)) = n
605  END DO
606  !
607  ! -- terminate if errors occurred
608  if (count_errors() > 0) then
609  call parser%StoreErrorUnit()
610  end if
simmodule
This module contains simulation methods.
Definition: Sim.f90:10
simmodule::store_error
subroutine, public store_error(msg, terminate)
Store an error message.
Definition: Sim.f90:92
simmodule::count_errors
integer(i4b) function, public count_errors()
Return number of errors.
Definition: Sim.f90:59
genrcm
subroutine genrcm(node_num, adj_num, adj_row, adj, perm)
Definition: rcm.f90:1004
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_cg()

subroutine imslinearbasemodule::ims_base_cg ( integer(i4b), intent(inout)  ICNVG,
integer(i4b), intent(in)  ITMAX,
integer(i4b), intent(inout)  INNERIT,
integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NIAPC,
integer(i4b), intent(in)  NJAPC,
integer(i4b), intent(in)  IPC,
integer(i4b), intent(in)  ICNVGOPT,
integer(i4b), intent(in)  NORTH,
real(dp), intent(in)  DVCLOSE,
real(dp), intent(in)  RCLOSE,
real(dp), intent(in)  L2NORM0,
real(dp), intent(in)  EPFACT,
integer(i4b), dimension(neq + 1), intent(in)  IA0,
integer(i4b), dimension(nja), intent(in)  JA0,
real(dp), dimension(nja), intent(in)  A0,
integer(i4b), dimension(niapc + 1), intent(in)  IAPC,
integer(i4b), dimension(njapc), intent(in)  JAPC,
real(dp), dimension(njapc), intent(in)  APC,
real(dp), dimension(neq), intent(inout)  X,
real(dp), dimension(neq), intent(inout)  B,
real(dp), dimension(neq), intent(inout)  D,
real(dp), dimension(neq), intent(inout)  P,
real(dp), dimension(neq), intent(inout)  Q,
real(dp), dimension(neq), intent(inout)  Z,
integer(i4b), intent(in)  NJLU,
integer(i4b), dimension(niapc), intent(in)  IW,
integer(i4b), dimension(njlu), intent(in)  JLU,
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), intent(in), pointer  summary 
)

Apply the Preconditioned Conjugate Gradient linear accelerator to the current coefficient matrix, right-hand side, using the current dependent-variable.

Parameters
[in,out]icnvgconvergence flag (1) non-convergence (0)
[in]itmaxmaximum number of inner iterations
[in,out]inneritinner iteration count
[in]neqnumber of equations
[in]njanumber of non-zero entries
[in]niapcpreconditioner number of rows
[in]njapcpreconditioner number of non-zero entries
[in]ipcpreconditioner option
[in]icnvgoptflow convergence criteria option
[in]northorthogonalization frequency
[in]dvclosedependent-variable closure criteria
[in]rcloseflow closure criteria
[in]l2norm0initial L-2 norm for system of equations
[in]epfactfactor for decreasing flow convergence criteria for subsequent Picard iterations
[in]ia0CRS row pointers
[in]ja0CRS column pointers
[in]a0coefficient matrix
[in]iapcpreconditioner CRS row pointers
[in]japcpreconditioner CRS column pointers
[in]apcpreconditioner matrix
[in,out]xdependent-variable vector
[in,out]bright-hand side vector
[in,out]dworking vector
[in,out]pworking vector
[in,out]qworking vector
[in,out]zworking vector
[in]njlupreconditioner length of JLU vector
[in]iwpreconditioner integer working vector
[in]jlupreconditioner JLU working vector
[in]nconvmaximum number of inner iterations in a time step (maxiter * maxinner)
[in]convnmodnumber of models in the solution
[in,out]convmodstartpointer to the start of each model in the convmod* arrays
[in,out]caccelconvergence string
[in]summaryConvergence summary report

Definition at line 31 of file ImsLinearBase.f90.

40  ! -- dummy variables
41  integer(I4B), INTENT(INOUT) :: ICNVG !< convergence flag (1) non-convergence (0)
42  integer(I4B), INTENT(IN) :: ITMAX !< maximum number of inner iterations
43  integer(I4B), INTENT(INOUT) :: INNERIT !< inner iteration count
44  integer(I4B), INTENT(IN) :: NEQ !< number of equations
45  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
46  integer(I4B), INTENT(IN) :: NIAPC !< preconditioner number of rows
47  integer(I4B), INTENT(IN) :: NJAPC !< preconditioner number of non-zero entries
48  integer(I4B), INTENT(IN) :: IPC !< preconditioner option
49  integer(I4B), INTENT(IN) :: ICNVGOPT !< flow convergence criteria option
50  integer(I4B), INTENT(IN) :: NORTH !< orthogonalization frequency
51  real(DP), INTENT(IN) :: DVCLOSE !< dependent-variable closure criteria
52  real(DP), INTENT(IN) :: RCLOSE !< flow closure criteria
53  real(DP), INTENT(IN) :: L2NORM0 !< initial L-2 norm for system of equations
54  real(DP), INTENT(IN) :: EPFACT !< factor for decreasing flow convergence criteria for subsequent Picard iterations
55  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA0 !< CRS row pointers
56  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA0 !< CRS column pointers
57  real(DP), DIMENSION(NJA), INTENT(IN) :: A0 !< coefficient matrix
58  integer(I4B), DIMENSION(NIAPC + 1), INTENT(IN) :: IAPC !< preconditioner CRS row pointers
59  integer(I4B), DIMENSION(NJAPC), INTENT(IN) :: JAPC !< preconditioner CRS column pointers
60  real(DP), DIMENSION(NJAPC), INTENT(IN) :: APC !< preconditioner matrix
61  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: X !< dependent-variable vector
62  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: B !< right-hand side vector
63  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: D !< working vector
64  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: P !< working vector
65  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: Q !< working vector
66  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: Z !< working vector
67  ! -- ILUT dummy variables
68  integer(I4B), INTENT(IN) :: NJLU !< preconditioner length of JLU vector
69  integer(I4B), DIMENSION(NIAPC), INTENT(IN) :: IW !< preconditioner integer working vector
70  integer(I4B), DIMENSION(NJLU), INTENT(IN) :: JLU !< preconditioner JLU working vector
71  ! -- convergence information dummy variables dummy variables
72  integer(I4B), INTENT(IN) :: NCONV !< maximum number of inner iterations in a time step (maxiter * maxinner)
73  integer(I4B), INTENT(IN) :: CONVNMOD !< number of models in the solution
74  integer(I4B), DIMENSION(CONVNMOD + 1), INTENT(INOUT) :: CONVMODSTART !< pointer to the start of each model in the convmod* arrays
75  character(len=31), DIMENSION(NCONV), INTENT(INOUT) :: CACCEL !< convergence string
76  type(ConvergenceSummaryType), pointer, intent(in) :: summary !< Convergence summary report
77  ! -- local variables
78  LOGICAL :: lorth
79  logical :: lsame
80  character(len=31) :: cval
81  integer(I4B) :: n
82  integer(I4B) :: iiter
83  integer(I4B) :: xloc, rloc
84  integer(I4B) :: im, im0, im1
85  real(DP) :: ddot
86  real(DP) :: tv
87  real(DP) :: deltax
88  real(DP) :: rmax
89  real(DP) :: l2norm
90  real(DP) :: rcnvg
91  real(DP) :: denominator
92  real(DP) :: alpha, beta
93  real(DP) :: rho, rho0
94  !
95  ! -- initialize local variables
96  rho0 = dzero
97  rho = dzero
98  innerit = 0
99  !
100  ! -- INNER ITERATION
101  inner: DO iiter = 1, itmax
102  innerit = innerit + 1
103  summary%iter_cnt = summary%iter_cnt + 1
104  !
105  ! -- APPLY PRECONDITIONER
106  SELECT CASE (ipc)
107  !
108  ! -- ILU0 AND MILU0
109  CASE (ipc_ilu0, ipc_milu0)
110  CALL ims_base_ilu0a(nja, neq, apc, iapc, japc, d, z)
111  !
112  ! -- ILUT AND MILUT
113  CASE (ipc_ilut, ipc_milut)
114  CALL lusol(neq, d, z, apc, jlu, iw)
115  END SELECT
116  rho = ddot(neq, d, 1, z, 1)
117  !
118  ! -- COMPUTE DIRECTIONAL VECTORS
119  IF (iiter == 1) THEN
120  DO n = 1, neq
121  p(n) = z(n)
122  END DO
123  ELSE
124  beta = rho / rho0
125  DO n = 1, neq
126  p(n) = z(n) + beta * p(n)
127  END DO
128  END IF
129  !
130  ! -- COMPUTE ITERATES
131  !
132  ! -- UPDATE Q
133  call amux(neq, p, q, a0, ja0, ia0)
134  denominator = ddot(neq, p, 1, q, 1)
135  denominator = denominator + sign(dprec, denominator)
136  alpha = rho / denominator
137  !
138  ! -- UPDATE X AND RESIDUAL
139  deltax = dzero
140  rmax = dzero
141  l2norm = dzero
142  DO im = 1, convnmod
143  summary%locdv(im) = 0
144  summary%dvmax(im) = dzero
145  summary%locr(im) = 0
146  summary%rmax(im) = dzero
147  END DO
148  im = 1
149  im0 = convmodstart(1)
150  im1 = convmodstart(2)
151  DO n = 1, neq
152  !
153  ! -- determine current model index
154  if (n == im1) then
155  im = im + 1
156  im0 = convmodstart(im)
157  im1 = convmodstart(im + 1)
158  end if
159  !
160  ! -- identify deltax and rmax
161  tv = alpha * p(n)
162  x(n) = x(n) + tv
163  IF (abs(tv) > abs(deltax)) THEN
164  deltax = tv
165  xloc = n
166  END IF
167  IF (abs(tv) > abs(summary%dvmax(im))) THEN
168  summary%dvmax(im) = tv
169  summary%locdv(im) = n
170  END IF
171  tv = d(n)
172  tv = tv - alpha * q(n)
173  d(n) = tv
174  IF (abs(tv) > abs(rmax)) THEN
175  rmax = tv
176  rloc = n
177  END IF
178  IF (abs(tv) > abs(summary%rmax(im))) THEN
179  summary%rmax(im) = tv
180  summary%locr(im) = n
181  END IF
182  l2norm = l2norm + tv * tv
183  END DO
184  l2norm = sqrt(l2norm)
185  !
186  ! -- SAVE SOLVER convergence information dummy variables
187  IF (nconv > 1) THEN !<
188  n = summary%iter_cnt
189  WRITE (cval, '(g15.7)') alpha
190  caccel(n) = cval
191  summary%itinner(n) = iiter
192  DO im = 1, convnmod
193  summary%convlocdv(im, n) = summary%locdv(im)
194  summary%convlocr(im, n) = summary%locr(im)
195  summary%convdvmax(im, n) = summary%dvmax(im)
196  summary%convrmax(im, n) = summary%rmax(im)
197  END DO
198  END IF
199  !
200  ! -- TEST FOR SOLVER CONVERGENCE
201  IF (icnvgopt == 2 .OR. icnvgopt == 3 .OR. icnvgopt == 4) THEN
202  rcnvg = l2norm
203  ELSE
204  rcnvg = rmax
205  END IF
206  CALL ims_base_testcnvg(icnvgopt, icnvg, innerit, &
207  deltax, rcnvg, &
208  l2norm0, epfact, dvclose, rclose)
209  !
210  ! -- CHECK FOR EXACT SOLUTION
211  IF (rcnvg == dzero) icnvg = 1
212  !
213  ! -- CHECK FOR STANDARD CONVERGENCE
214  IF (icnvg .NE. 0) EXIT inner
215  !
216  ! -- CHECK THAT CURRENT AND PREVIOUS rho ARE DIFFERENT
217  lsame = is_close(rho, rho0)
218  IF (lsame) THEN
219  EXIT inner
220  END IF
221  !
222  ! -- RECALCULATE THE RESIDUAL
223  IF (north > 0) THEN
224  lorth = mod(iiter + 1, north) == 0
225  IF (lorth) THEN
226  call ims_base_residual(neq, nja, x, b, d, a0, ia0, ja0)
227  END IF
228  END IF
229  !
230  ! -- exit inner if rho is zero
231  if (rho == dzero) then
232  exit inner
233  end if
234  !
235  ! -- SAVE CURRENT INNER ITERATES
236  rho0 = rho
237  END DO inner
238  !
239  ! -- RESET ICNVG
240  IF (icnvg < 0) icnvg = 0
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_epfact()

real(dp) function imslinearbasemodule::ims_base_epfact ( integer(i4b)  icnvgopt,
integer(i4b)  kstp 
)
Parameters
icnvgoptIMS convergence option
kstptime step number
Returns
factor for decreasing convergence criteria in subsequent Picard iterations

Definition at line 1317 of file ImsLinearBase.f90.

1318  integer(I4B) :: icnvgopt !< IMS convergence option
1319  integer(I4B) :: kstp !< time step number
1320  real(DP) :: epfact !< factor for decreasing convergence criteria in subsequent Picard iterations
1321 
1322  if (icnvgopt == 2) then
1323  if (kstp == 1) then
1324  epfact = 0.01
1325  else
1326  epfact = 0.10
1327  end if
1328  else if (icnvgopt == 4) then
1329  epfact = dem4
1330  else
1331  epfact = done
1332  end if
1333 
Here is the caller graph for this function:

◆ ims_base_ilu0a()

subroutine imslinearbasemodule::ims_base_ilu0a ( integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NEQ,
real(dp), dimension(nja), intent(in)  APC,
integer(i4b), dimension(neq + 1), intent(in)  IAPC,
integer(i4b), dimension(nja), intent(in)  JAPC,
real(dp), dimension(neq), intent(in)  R,
real(dp), dimension(neq), intent(inout)  D 
)

Apply the ILU0 and MILU0 preconditioners to the passed vector (R).

Parameters
[in]njanumber of non-zero entries
[in]neqnumber of equations
[in]apcILU0/MILU0 preconditioner matrix
[in]iapcILU0/MILU0 preconditioner CRS row pointers
[in]japcILU0/MILU0 preconditioner CRS column pointers
[in]rinput vector
[in,out]doutput vector after applying APC to R

Definition at line 1050 of file ImsLinearBase.f90.

1051  ! -- dummy variables
1052  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
1053  integer(I4B), INTENT(IN) :: NEQ !< number of equations
1054  real(DP), DIMENSION(NJA), INTENT(IN) :: APC !< ILU0/MILU0 preconditioner matrix
1055  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IAPC !< ILU0/MILU0 preconditioner CRS row pointers
1056  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JAPC !< ILU0/MILU0 preconditioner CRS column pointers
1057  real(DP), DIMENSION(NEQ), INTENT(IN) :: R !< input vector
1058  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: D !< output vector after applying APC to R
1059  ! -- local variables
1060  integer(I4B) :: ic0, ic1
1061  integer(I4B) :: iu
1062  integer(I4B) :: jcol
1063  integer(I4B) :: j, n
1064  real(DP) :: tv
1065  !
1066  ! -- FORWARD SOLVE - APC * D = R
1067  forward: DO n = 1, neq
1068  tv = r(n)
1069  ic0 = iapc(n)
1070  ic1 = iapc(n + 1) - 1
1071  iu = japc(n) - 1
1072  lower: DO j = ic0, iu
1073  jcol = japc(j)
1074  tv = tv - apc(j) * d(jcol)
1075  END DO lower
1076  d(n) = tv
1077  END DO forward
1078  !
1079  ! -- BACKWARD SOLVE - D = D / U
1080  backward: DO n = neq, 1, -1
1081  ic0 = iapc(n)
1082  ic1 = iapc(n + 1) - 1
1083  iu = japc(n)
1084  tv = d(n)
1085  upper: DO j = iu, ic1
1086  jcol = japc(j)
1087  tv = tv - apc(j) * d(jcol)
1088  END DO upper
1089  !
1090  ! -- COMPUTE D FOR DIAGONAL - D = D / U
1091  d(n) = tv * apc(n)
1092  END DO backward
Here is the caller graph for this function:

◆ ims_base_isort()

subroutine imslinearbasemodule::ims_base_isort ( integer(i4b), intent(in)  NVAL,
integer(i4b), dimension(nval), intent(inout)  IARRAY 
)

Subroutine sort an integer array in-place.

Parameters
[in]nvallength of the integer array
[in,out]iarrayinteger array to be sorted

Definition at line 1269 of file ImsLinearBase.f90.

1270  ! -- dummy variables
1271  integer(I4B), INTENT(IN) :: NVAL !< length of the integer array
1272  integer(I4B), DIMENSION(NVAL), INTENT(INOUT) :: IARRAY !< integer array to be sorted
1273  ! -- local variables
1274  integer(I4B) :: i, j, itemp
1275  ! -- code
1276  DO i = 1, nval - 1
1277  DO j = i + 1, nval
1278  if (iarray(i) > iarray(j)) then
1279  itemp = iarray(j)
1280  iarray(j) = iarray(i)
1281  iarray(i) = itemp
1282  END IF
1283  END DO
1284  END DO
Here is the caller graph for this function:

◆ ims_base_jaca()

subroutine imslinearbasemodule::ims_base_jaca ( integer(i4b), intent(in)  NEQ,
real(dp), dimension(neq), intent(in)  A,
real(dp), dimension(neq), intent(in)  D1,
real(dp), dimension(neq), intent(inout)  D2 
)

Apply the Jacobi preconditioner and return the resultant vector.

Parameters
[in]neqnumber of equations
[in]aJacobi preconditioner
[in]d1input vector
[in,out]d2resultant vector

Definition at line 908 of file ImsLinearBase.f90.

909  ! -- dummy variables
910  integer(I4B), INTENT(IN) :: NEQ !< number of equations
911  real(DP), DIMENSION(NEQ), INTENT(IN) :: A !< Jacobi preconditioner
912  real(DP), DIMENSION(NEQ), INTENT(IN) :: D1 !< input vector
913  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: D2 !< resultant vector
914  ! -- local variables
915  integer(I4B) :: n
916  real(DP) :: tv
917  ! -- code
918  DO n = 1, neq
919  tv = a(n) * d1(n)
920  d2(n) = tv
921  END DO

◆ ims_base_pccrs()

subroutine imslinearbasemodule::ims_base_pccrs ( integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA,
integer(i4b), dimension(neq + 1), intent(inout)  IAPC,
integer(i4b), dimension(nja), intent(inout)  JAPC 
)

Generate the CRS row and column pointers for the preconditioner. JAPC(1:NEQ) hHas the position of the upper entry for a row, JAPC(NEQ+1:NJA) is the column position for entry, APC(1:NEQ) is the preconditioned inverse of the diagonal, and APC(NEQ+1:NJA) are the preconditioned entries for off diagonals.

Definition at line 1208 of file ImsLinearBase.f90.

1210  ! -- dummy variables
1211  integer(I4B), INTENT(IN) :: NEQ !<
1212  integer(I4B), INTENT(IN) :: NJA !<
1213  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !<
1214  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !<
1215  integer(I4B), DIMENSION(NEQ + 1), INTENT(INOUT) :: IAPC !<
1216  integer(I4B), DIMENSION(NJA), INTENT(INOUT) :: JAPC !<
1217  ! -- local variables
1218  integer(I4B) :: n, j
1219  integer(I4B) :: i0, i1
1220  integer(I4B) :: nlen
1221  integer(I4B) :: ic, ip
1222  integer(I4B) :: jcol
1223  integer(I4B), DIMENSION(:), ALLOCATABLE :: iarr
1224  ! -- code
1225  ip = neq + 1
1226  DO n = 1, neq
1227  i0 = ia(n)
1228  i1 = ia(n + 1) - 1
1229  nlen = i1 - i0
1230  ALLOCATE (iarr(nlen))
1231  ic = 0
1232  DO j = i0, i1
1233  jcol = ja(j)
1234  IF (jcol == n) cycle
1235  ic = ic + 1
1236  iarr(ic) = jcol
1237  END DO
1238  CALL ims_base_isort(nlen, iarr)
1239  iapc(n) = ip
1240  DO j = 1, nlen
1241  jcol = iarr(j)
1242  japc(ip) = jcol
1243  ip = ip + 1
1244  END DO
1245  DEALLOCATE (iarr)
1246  END DO
1247  iapc(neq + 1) = nja + 1
1248  !
1249  ! -- POSITION OF THE FIRST UPPER ENTRY FOR ROW
1250  DO n = 1, neq
1251  i0 = iapc(n)
1252  i1 = iapc(n + 1) - 1
1253  japc(n) = iapc(n + 1)
1254  DO j = i0, i1
1255  jcol = japc(j)
1256  IF (jcol > n) THEN
1257  japc(n) = j
1258  EXIT
1259  END IF
1260  END DO
1261  END DO
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_pcilu0()

subroutine imslinearbasemodule::ims_base_pcilu0 ( integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NEQ,
real(dp), dimension(nja), intent(in)  AMAT,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA,
real(dp), dimension(nja), intent(inout)  APC,
integer(i4b), dimension(neq + 1), intent(inout)  IAPC,
integer(i4b), dimension(nja), intent(inout)  JAPC,
integer(i4b), dimension(neq), intent(inout)  IW,
real(dp), dimension(neq), intent(inout)  W,
real(dp), intent(in)  RELAX,
integer(i4b), intent(inout)  IPCFLAG,
real(dp), intent(in)  DELTA 
)

Update the ILU0 preconditioner using the current coefficient matrix.

Parameters
[in]njanumber of non-zero entries
[in]neqnumber of equations
[in]amatcoefficient matrix
[in]iaCRS row pointers
[in]jaCRS column pointers
[in,out]apcpreconditioned matrix
[in,out]iapcpreconditioner CRS row pointers
[in,out]japcpreconditioner CRS column pointers
[in,out]iwpreconditioner integer work vector
[in,out]wpreconditioner work vector
[in]relaxMILU0 preconditioner relaxation factor
[in,out]ipcflagpreconditioner error flag
[in]deltafactor used to correct non-diagonally dominant matrices

Definition at line 929 of file ImsLinearBase.f90.

932  ! -- dummy variables
933  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
934  integer(I4B), INTENT(IN) :: NEQ !< number of equations
935  real(DP), DIMENSION(NJA), INTENT(IN) :: AMAT !< coefficient matrix
936  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< CRS row pointers
937  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< CRS column pointers
938  real(DP), DIMENSION(NJA), INTENT(INOUT) :: APC !< preconditioned matrix
939  integer(I4B), DIMENSION(NEQ + 1), INTENT(INOUT) :: IAPC !< preconditioner CRS row pointers
940  integer(I4B), DIMENSION(NJA), INTENT(INOUT) :: JAPC !< preconditioner CRS column pointers
941  integer(I4B), DIMENSION(NEQ), INTENT(INOUT) :: IW !< preconditioner integer work vector
942  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: W !< preconditioner work vector
943  real(DP), INTENT(IN) :: RELAX !< MILU0 preconditioner relaxation factor
944  integer(I4B), INTENT(INOUT) :: IPCFLAG !< preconditioner error flag
945  real(DP), INTENT(IN) :: DELTA !< factor used to correct non-diagonally dominant matrices
946  ! -- local variables
947  integer(I4B) :: ic0, ic1
948  integer(I4B) :: iic0, iic1
949  integer(I4B) :: iu, iiu
950  integer(I4B) :: j, n
951  integer(I4B) :: jj
952  integer(I4B) :: jcol, jw
953  integer(I4B) :: jjcol
954  real(DP) :: drelax
955  real(DP) :: sd1
956  real(DP) :: tl
957  real(DP) :: rs
958  real(DP) :: d
959  !
960  ! -- initialize local variables
961  drelax = relax
962  DO n = 1, neq
963  iw(n) = 0
964  w(n) = dzero
965  END DO
966  main: DO n = 1, neq
967  ic0 = ia(n)
968  ic1 = ia(n + 1) - 1
969  DO j = ic0, ic1
970  jcol = ja(j)
971  iw(jcol) = 1
972  w(jcol) = w(jcol) + amat(j)
973  END DO
974  ic0 = iapc(n)
975  ic1 = iapc(n + 1) - 1
976  iu = japc(n)
977  rs = dzero
978  lower: DO j = ic0, iu - 1
979  jcol = japc(j)
980  iic0 = iapc(jcol)
981  iic1 = iapc(jcol + 1) - 1
982  iiu = japc(jcol)
983  tl = w(jcol) * apc(jcol)
984  w(jcol) = tl
985  DO jj = iiu, iic1
986  jjcol = japc(jj)
987  jw = iw(jjcol)
988  IF (jw .NE. 0) THEN
989  w(jjcol) = w(jjcol) - tl * apc(jj)
990  ELSE
991  rs = rs + tl * apc(jj)
992  END IF
993  END DO
994  END DO lower
995  !
996  ! -- DIAGONAL - CALCULATE INVERSE OF DIAGONAL FOR SOLUTION
997  d = w(n)
998  tl = (done + delta) * d - (drelax * rs)
999  !
1000  ! -- ENSURE THAT THE SIGN OF THE DIAGONAL HAS NOT CHANGED AND IS
1001  sd1 = sign(d, tl)
1002  IF (sd1 .NE. d) THEN
1003  !
1004  ! -- USE SMALL VALUE IF DIAGONAL SCALING IS NOT EFFECTIVE FOR
1005  ! PIVOTS THAT CHANGE THE SIGN OF THE DIAGONAL
1006  IF (ipcflag > 1) THEN
1007  tl = sign(dem6, d)
1008  !
1009  ! -- DIAGONAL SCALING CONTINUES TO BE EFFECTIVE
1010  ELSE
1011  ipcflag = 1
1012  EXIT main
1013  END IF
1014  END IF
1015  IF (abs(tl) == dzero) THEN
1016  !
1017  ! -- USE SMALL VALUE IF DIAGONAL SCALING IS NOT EFFECTIVE FOR
1018  ! ZERO PIVOTS
1019  IF (ipcflag > 1) THEN
1020  tl = sign(dem6, d)
1021  !
1022  ! -- DIAGONAL SCALING CONTINUES TO BE EFFECTIVE FOR ELIMINATING
1023  ELSE
1024  ipcflag = 1
1025  EXIT main
1026  END IF
1027  END IF
1028  apc(n) = done / tl
1029  !
1030  ! -- RESET POINTER FOR IW TO ZERO
1031  iw(n) = 0
1032  w(n) = dzero
1033  DO j = ic0, ic1
1034  jcol = japc(j)
1035  apc(j) = w(jcol)
1036  iw(jcol) = 0
1037  w(jcol) = dzero
1038  END DO
1039  END DO main
1040  !
1041  ! -- RESET IPCFLAG IF SUCCESSFUL COMPLETION OF MAIN
1042  ipcflag = 0
Here is the caller graph for this function:

◆ ims_base_pcjac()

subroutine imslinearbasemodule::ims_base_pcjac ( integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NEQ,
real(dp), dimension(nja), intent(in)  AMAT,
real(dp), dimension(neq), intent(inout)  APC,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA 
)

Calculate the Jacobi preconditioner (inverse of the diagonal) using the current coefficient matrix.

Parameters
[in]njanumber of non-zero entries
[in]neqnumber of equations
[in]amatcoefficient matrix
[in,out]apcpreconditioner matrix
[in]iaCRS row pointers
[in]jaCRS column pointers

Definition at line 873 of file ImsLinearBase.f90.

874  ! -- dummy variables
875  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
876  integer(I4B), INTENT(IN) :: NEQ !< number of equations
877  real(DP), DIMENSION(NJA), INTENT(IN) :: AMAT !< coefficient matrix
878  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: APC !< preconditioner matrix
879  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< CRS row pointers
880  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< CRS column pointers
881  ! -- local variables
882  integer(I4B) :: i, n
883  integer(I4B) :: ic0, ic1
884  integer(I4B) :: id
885  real(DP) :: tv
886  ! -- code
887  DO n = 1, neq
888  ic0 = ia(n)
889  ic1 = ia(n + 1) - 1
890  id = ia(n)
891  DO i = ic0, ic1
892  IF (ja(i) == n) THEN
893  id = i
894  EXIT
895  END IF
896  END DO
897  tv = amat(id)
898  IF (abs(tv) > dzero) tv = done / tv
899  apc(n) = tv
900  END DO

◆ ims_base_pcu()

subroutine imslinearbasemodule::ims_base_pcu ( integer(i4b), intent(in)  IOUT,
integer(i4b), intent(in)  NJA,
integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NIAPC,
integer(i4b), intent(in)  NJAPC,
integer(i4b), intent(in)  IPC,
real(dp), intent(in)  RELAX,
real(dp), dimension(nja), intent(in)  AMAT,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA,
real(dp), dimension(njapc), intent(inout)  APC,
integer(i4b), dimension(niapc + 1), intent(inout)  IAPC,
integer(i4b), dimension(njapc), intent(inout)  JAPC,
integer(i4b), dimension(niapc), intent(inout)  IW,
real(dp), dimension(niapc), intent(inout)  W,
integer(i4b), intent(in)  LEVEL,
real(dp), intent(in)  DROPTOL,
integer(i4b), intent(in)  NJLU,
integer(i4b), intent(in)  NJW,
integer(i4b), intent(in)  NWLU,
integer(i4b), dimension(njlu), intent(inout)  JLU,
integer(i4b), dimension(njw), intent(inout)  JW,
real(dp), dimension(nwlu), intent(inout)  WLU 
)

Update the preconditioner using the current coefficient matrix.

Parameters
[in]ioutsimulation listing file unit
[in]njanumber of non-zero entries
[in]neqnumber of equations
[in]niapcpreconditioner number of rows
[in]njapcpreconditioner number of non-zero entries
[in]ipcpreconditioner (1) ILU0 (2) MILU0 (3) ILUT (4) MILUT
[in]relaxpreconditioner relaxation factor for MILU0 and MILUT
[in]amatcoefficient matrix
[in]iaCRS row pointers
[in]jaCRS column pointers
[in,out]apcpreconditioner matrix
[in,out]iapcpreconditioner CRS row pointers
[in,out]japcpreconditioner CRS column pointers
[in,out]iwpreconditioner integer work vector
[in,out]wpreconditioner work vector
[in]levelnumber of levels of fill for ILUT and MILUT
[in]droptoldrop tolerance
[in]njlulength of JLU working vector
[in]njwlength of JW working vector
[in]nwlulength of WLU working vector
[in,out]jluILUT/MILUT JLU working vector
[in,out]jwILUT/MILUT JW working vector
[in,out]wluILUT/MILUT WLU working vector

Definition at line 762 of file ImsLinearBase.f90.

765  ! -- modules
766  use simmodule, only: store_error, count_errors
767  ! -- dummy variables
768  integer(I4B), INTENT(IN) :: IOUT !< simulation listing file unit
769  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
770  integer(I4B), INTENT(IN) :: NEQ !< number of equations
771  integer(I4B), INTENT(IN) :: NIAPC !< preconditioner number of rows
772  integer(I4B), INTENT(IN) :: NJAPC !< preconditioner number of non-zero entries
773  integer(I4B), INTENT(IN) :: IPC !< preconditioner (1) ILU0 (2) MILU0 (3) ILUT (4) MILUT
774  real(DP), INTENT(IN) :: RELAX !< preconditioner relaxation factor for MILU0 and MILUT
775  real(DP), DIMENSION(NJA), INTENT(IN) :: AMAT !< coefficient matrix
776  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< CRS row pointers
777  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< CRS column pointers
778  real(DP), DIMENSION(NJAPC), INTENT(INOUT) :: APC !< preconditioner matrix
779  integer(I4B), DIMENSION(NIAPC + 1), INTENT(INOUT) :: IAPC !< preconditioner CRS row pointers
780  integer(I4B), DIMENSION(NJAPC), INTENT(INOUT) :: JAPC !< preconditioner CRS column pointers
781  integer(I4B), DIMENSION(NIAPC), INTENT(INOUT) :: IW !< preconditioner integer work vector
782  real(DP), DIMENSION(NIAPC), INTENT(INOUT) :: W !< preconditioner work vector
783  ! -- ILUT dummy variables
784  integer(I4B), INTENT(IN) :: LEVEL !< number of levels of fill for ILUT and MILUT
785  real(DP), INTENT(IN) :: DROPTOL !< drop tolerance
786  integer(I4B), INTENT(IN) :: NJLU !< length of JLU working vector
787  integer(I4B), INTENT(IN) :: NJW !< length of JW working vector
788  integer(I4B), INTENT(IN) :: NWLU !< length of WLU working vector
789  integer(I4B), DIMENSION(NJLU), INTENT(INOUT) :: JLU !< ILUT/MILUT JLU working vector
790  integer(I4B), DIMENSION(NJW), INTENT(INOUT) :: JW !< ILUT/MILUT JW working vector
791  real(DP), DIMENSION(NWLU), INTENT(INOUT) :: WLU !< ILUT/MILUT WLU working vector
792  ! -- local variables
793  character(len=LINELENGTH) :: errmsg
794  character(len=100), dimension(5), parameter :: cerr = &
795  ["Elimination process has generated a row in L or U whose length is > n.", &
796  &"The matrix L overflows the array al. ", &
797  &"The matrix U overflows the array alu. ", &
798  &"Illegal value for lfil. ", &
799  &"Zero row encountered. "]
800  integer(i4b) :: ipcflag
801  integer(I4B) :: icount
802  integer(I4B) :: ierr
803  real(DP) :: delta
804  ! -- formats
805 2000 FORMAT(/, ' MATRIX IS SEVERELY NON-DIAGONALLY DOMINANT.', &
806  /, ' ADDED SMALL VALUE TO PIVOT ', i0, ' TIMES IN', &
807  ' IMSLINEARSUB_PCU.')
808  !
809  ! -- initialize local variables
810  ipcflag = 0
811  icount = 0
812  delta = dzero
813  pcscale: DO
814  SELECT CASE (ipc)
815  !
816  ! -- ILU0 AND MILU0
817  CASE (ipc_ilu0, ipc_milu0)
818  CALL ims_base_pcilu0(nja, neq, amat, ia, ja, &
819  apc, iapc, japc, iw, w, &
820  relax, ipcflag, delta)
821  !
822  ! -- ILUT AND MILUT
823  CASE (ipc_ilut, ipc_milut)
824  ierr = 0
825  CALL ilut(neq, amat, ja, ia, level, droptol, &
826  apc, jlu, iw, njapc, wlu, jw, ierr, &
827  relax, ipcflag, delta)
828  if (ierr /= 0) then
829  if (ierr > 0) then
830  write (errmsg, '(a,1x,i0,1x,a)') &
831  'ILUT: zero pivot encountered at step number', ierr, '.'
832  else
833  write (errmsg, '(a,1x,a)') 'ILUT:', cerr(-ierr)
834  end if
835  call store_error(errmsg)
836  call parser%StoreErrorUnit()
837  end if
838  !
839  ! -- ADDITIONAL PRECONDITIONERS
840  CASE DEFAULT
841  ipcflag = 0
842  END SELECT
843  IF (ipcflag < 1) THEN
844  EXIT pcscale
845  END IF
846  delta = 1.5d0 * delta + dem3
847  ipcflag = 0
848  IF (delta > dhalf) THEN
849  delta = dhalf
850  ipcflag = 2
851  END IF
852  icount = icount + 1
853  !
854  ! -- terminate pcscale loop if not making progress
855  if (icount > 10) then
856  exit pcscale
857  end if
858 
859  END DO pcscale
860  !
861  ! -- write error message if small value added to pivot
862  if (icount > 0) then
863  write (iout, 2000) icount
864  end if
ilut
subroutine ilut(n, a, ja, ia, lfil, droptol, alu, jlu, ju, iwk, w, jw, ierr, relax, izero, delta)
Definition: ilut.f90:51
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_residual()

subroutine imslinearbasemodule::ims_base_residual ( integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
real(dp), dimension(neq), intent(in)  X,
real(dp), dimension(neq), intent(in)  B,
real(dp), dimension(neq), intent(inout)  D,
real(dp), dimension(nja), intent(in)  A,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA 
)

Subroutine to calculate the residual.

Parameters
[in]neqlength of vectors
[in]njalength of coefficient matrix
[in]xdependent variable
[in]bright-hand side
[in,out]dresidual
[in]acoefficient matrix
[in]iaCRS row pointers
[in]jaCRS column pointers

Definition at line 1292 of file ImsLinearBase.f90.

1293  ! -- dummy variables
1294  integer(I4B), INTENT(IN) :: NEQ !< length of vectors
1295  integer(I4B), INTENT(IN) :: NJA !< length of coefficient matrix
1296  real(DP), DIMENSION(NEQ), INTENT(IN) :: X !< dependent variable
1297  real(DP), DIMENSION(NEQ), INTENT(IN) :: B !< right-hand side
1298  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: D !< residual
1299  real(DP), DIMENSION(NJA), INTENT(IN) :: A !< coefficient matrix
1300  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< CRS row pointers
1301  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< CRS column pointers
1302  ! -- local variables
1303  integer(I4B) :: n
1304  ! -- code
1305  !
1306  ! -- calculate matrix-vector product
1307  call amux(neq, x, d, a, ja, ia)
1308  !
1309  ! -- subtract matrix-vector product from right-hand side
1310  DO n = 1, neq
1311  d(n) = b(n) - d(n)
1312  END DO
Here is the call graph for this function:
Here is the caller graph for this function:

◆ ims_base_scale()

subroutine imslinearbasemodule::ims_base_scale ( integer(i4b), intent(in)  IOPT,
integer(i4b), intent(in)  ISCL,
integer(i4b), intent(in)  NEQ,
integer(i4b), intent(in)  NJA,
integer(i4b), dimension(neq + 1), intent(in)  IA,
integer(i4b), dimension(nja), intent(in)  JA,
real(dp), dimension(nja), intent(inout)  AMAT,
real(dp), dimension(neq), intent(inout)  X,
real(dp), dimension(neq), intent(inout)  B,
real(dp), dimension(neq), intent(inout)  DSCALE,
real(dp), dimension(neq), intent(inout)  DSCALE2 
)

Scale the coefficient matrix (AMAT), the right-hand side (B), and the estimate of the dependent variable (X).

Parameters
[in]ioptflag to scale (0) or unscale the system of equations
[in]isclscaling option (1) symmetric (2) L-2 norm
[in]neqnumber of equations
[in]njanumber of non-zero entries
[in]iaCRS row pointer
[in]jaCRS column pointer
[in,out]amatcoefficient matrix
[in,out]xdependent variable
[in,out]bright-hand side
[in,out]dscalefirst scaling vector
[in,out]dscale2second scaling vector

Definition at line 620 of file ImsLinearBase.f90.

622  ! -- dummy variables
623  integer(I4B), INTENT(IN) :: IOPT !< flag to scale (0) or unscale the system of equations
624  integer(I4B), INTENT(IN) :: ISCL !< scaling option (1) symmetric (2) L-2 norm
625  integer(I4B), INTENT(IN) :: NEQ !< number of equations
626  integer(I4B), INTENT(IN) :: NJA !< number of non-zero entries
627  integer(I4B), DIMENSION(NEQ + 1), INTENT(IN) :: IA !< CRS row pointer
628  integer(I4B), DIMENSION(NJA), INTENT(IN) :: JA !< CRS column pointer
629  real(DP), DIMENSION(NJA), INTENT(INOUT) :: AMAT !< coefficient matrix
630  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: X !< dependent variable
631  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: B !< right-hand side
632  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: DSCALE !< first scaling vector
633  real(DP), DIMENSION(NEQ), INTENT(INOUT) :: DSCALE2 !< second scaling vector
634  ! -- local variables
635  integer(I4B) :: i, n
636  integer(I4B) :: id, jc
637  integer(I4B) :: i0, i1
638  real(DP) :: v, c1, c2
639  !
640  ! -- SCALE SCALE AMAT, X, AND B
641  IF (iopt == 0) THEN
642  !
643  ! -- SYMMETRIC SCALING
644  SELECT CASE (iscl)
645  CASE (1)
646  DO n = 1, neq
647  id = ia(n)
648  v = amat(id)
649  c1 = done / sqrt(abs(v))
650  dscale(n) = c1
651  dscale2(n) = c1
652  END DO
653  !
654  ! -- SCALE AMAT -- AMAT = DSCALE(row) * AMAT(i) * DSCALE2(col)
655  DO n = 1, neq
656  c1 = dscale(n)
657  i0 = ia(n)
658  i1 = ia(n + 1) - 1
659  DO i = i0, i1
660  jc = ja(i)
661  c2 = dscale2(jc)
662  amat(i) = c1 * amat(i) * c2
663  END DO
664  END DO
665  !
666  ! -- L-2 NORM SCALING
667  CASE (2)
668  !
669  ! -- SCALE EACH ROW SO THAT THE L-2 NORM IS 1
670  DO n = 1, neq
671  c1 = dzero
672  i0 = ia(n)
673  i1 = ia(n + 1) - 1
674  DO i = i0, i1
675  c1 = c1 + amat(i) * amat(i)
676  END DO
677  c1 = sqrt(c1)
678  IF (c1 == dzero) THEN
679  c1 = done
680  ELSE
681  c1 = done / c1
682  END IF
683  dscale(n) = c1
684  !
685  ! -- INITIAL SCALING OF AMAT -- AMAT = DSCALE(row) * AMAT(i)
686  DO i = i0, i1
687  amat(i) = c1 * amat(i)
688  END DO
689  END DO
690  !
691  ! -- SCALE EACH COLUMN SO THAT THE L-2 NORM IS 1
692  DO n = 1, neq
693  dscale2(n) = dzero
694  END DO
695  c2 = dzero
696  DO n = 1, neq
697  i0 = ia(n)
698  i1 = ia(n + 1) - 1
699  DO i = i0, i1
700  jc = ja(i)
701  c2 = amat(i)
702  dscale2(jc) = dscale2(jc) + c2 * c2
703  END DO
704  END DO
705  DO n = 1, neq
706  c2 = dscale2(n)
707  IF (c2 == dzero) THEN
708  c2 = done
709  ELSE
710  c2 = done / sqrt(c2)
711  END IF
712  dscale2(n) = c2
713  END DO
714  !
715  ! -- FINAL SCALING OF AMAT -- AMAT = DSCALE2(col) * AMAT(i)
716  DO n = 1, neq
717  i0 = ia(n)
718  i1 = ia(n + 1) - 1
719  DO i = i0, i1
720  jc = ja(i)
721  c2 = dscale2(jc)
722  amat(i) = c2 * amat(i)
723  END DO
724  END DO
725  END SELECT
726  !
727  ! -- SCALE X AND B
728  DO n = 1, neq
729  c1 = dscale(n)
730  c2 = dscale2(n)
731  x(n) = x(n) / c2
732  b(n) = b(n) * c1
733  END DO
734  !
735  ! -- UNSCALE SCALE AMAT, X, AND B
736  ELSE
737  DO n = 1, neq
738  c1 = dscale(n)
739  i0 = ia(n)
740  i1 = ia(n + 1) - 1
741  !
742  ! -- UNSCALE AMAT
743  DO i = i0, i1
744  jc = ja(i)
745  c2 = dscale2(jc)
746  amat(i) = (done / c1) * amat(i) * (done / c2)
747  END DO
748  !
749  ! -- UNSCALE X AND B
750  c2 = dscale2(n)
751  x(n) = x(n) * c2
752  b(n) = b(n) / c1
753  END DO
754  END IF
Here is the caller graph for this function:

◆ ims_base_testcnvg()

subroutine imslinearbasemodule::ims_base_testcnvg ( integer(i4b), intent(in)  Icnvgopt,
integer(i4b), intent(inout)  Icnvg,
integer(i4b), intent(in)  Iiter,
real(dp), intent(in)  Dvmax,
real(dp), intent(in)  Rmax,
real(dp), intent(in)  Rmax0,
real(dp), intent(in)  Epfact,
real(dp), intent(in)  Dvclose,
real(dp), intent(in)  Rclose 
)

General routine for testing for solver convergence based on the user-specified convergence option (Icnvgopt).

Parameters
[in]icnvgoptconvergence option - see documentation for option
[in,out]icnvgflag indicating if convergence achieved (1) or not (0)
[in]iiterinner iteration number (used for strict convergence option)
[in]dvmaxmaximum dependent-variable change
[in]rmaxmaximum flow change
[in]rmax0initial flow change (initial L2-norm)
[in]epfactfactor for reducing convergence criteria in subsequent Picard iterations
[in]dvcloseMaximum depenendent-variable change allowed
[in]rcloseMaximum flow change allowed

Definition at line 1102 of file ImsLinearBase.f90.

1105  ! -- dummy variables
1106  integer(I4B), INTENT(IN) :: Icnvgopt !< convergence option - see documentation for option
1107  integer(I4B), INTENT(INOUT) :: Icnvg !< flag indicating if convergence achieved (1) or not (0)
1108  integer(I4B), INTENT(IN) :: Iiter !< inner iteration number (used for strict convergence option)
1109  real(DP), INTENT(IN) :: Dvmax !< maximum dependent-variable change
1110  real(DP), INTENT(IN) :: Rmax !< maximum flow change
1111  real(DP), INTENT(IN) :: Rmax0 !< initial flow change (initial L2-norm)
1112  real(DP), INTENT(IN) :: Epfact !< factor for reducing convergence criteria in subsequent Picard iterations
1113  real(DP), INTENT(IN) :: Dvclose !< Maximum depenendent-variable change allowed
1114  real(DP), INTENT(IN) :: Rclose !< Maximum flow change allowed
1115  ! -- code
1116  IF (icnvgopt == 0) THEN
1117  IF (abs(dvmax) <= dvclose .AND. abs(rmax) <= rclose) THEN
1118  icnvg = 1
1119  END IF
1120  ELSE IF (icnvgopt == 1) THEN
1121  IF (abs(dvmax) <= dvclose .AND. abs(rmax) <= rclose) THEN
1122  IF (iiter == 1) THEN
1123  icnvg = 1
1124  ELSE
1125  icnvg = -1
1126  END IF
1127  END IF
1128  ELSE IF (icnvgopt == 2) THEN
1129  IF (abs(dvmax) <= dvclose .OR. rmax <= rclose) THEN
1130  icnvg = 1
1131  ELSE IF (rmax <= rmax0 * epfact) THEN
1132  icnvg = -1
1133  END IF
1134  ELSE IF (icnvgopt == 3) THEN
1135  IF (abs(dvmax) <= dvclose) THEN
1136  icnvg = 1
1137  ELSE IF (rmax <= rmax0 * rclose) THEN
1138  icnvg = -1
1139  END IF
1140  ELSE IF (icnvgopt == 4) THEN
1141  IF (abs(dvmax) <= dvclose .AND. rmax <= rclose) THEN
1142  icnvg = 1
1143  ELSE IF (rmax <= rmax0 * epfact) THEN
1144  icnvg = -1
1145  END IF
1146  END IF
Here is the caller graph for this function:

◆ ims_calc_pcdims()

subroutine imslinearbasemodule::ims_calc_pcdims ( integer(i4b), intent(in)  neq,
integer(i4b), intent(in)  nja,
integer(i4b), dimension(:), intent(in)  ia,
integer(i4b), intent(in)  level,
integer(i4b), intent(in)  ipc,
integer(i4b), intent(inout)  niapc,
integer(i4b), intent(inout)  njapc,
integer(i4b), intent(inout)  njlu,
integer(i4b), intent(inout)  njw,
integer(i4b), intent(inout)  nwlu 
)
Parameters
[in]neqnr. of rows A
[in]njanr. of nonzeros A
[in]iaCSR row pointers A
[in]levelfill level ILU
[in]ipcIMS preconditioner type
[in,out]niapcwork array size
[in,out]njapcwork array size
[in,out]njluwork array size
[in,out]njwwork array size
[in,out]nwluwork array size

Definition at line 1149 of file ImsLinearBase.f90.

1151  integer(I4B), intent(in) :: neq !< nr. of rows A
1152  integer(I4B), intent(in) :: nja !< nr. of nonzeros A
1153  integer(I4B), dimension(:), intent(in) :: ia !< CSR row pointers A
1154  integer(I4B), intent(in) :: level !< fill level ILU
1155  integer(I4B), intent(in) :: ipc !< IMS preconditioner type
1156  integer(I4B), intent(inout) :: niapc !< work array size
1157  integer(I4B), intent(inout) :: njapc !< work array size
1158  integer(I4B), intent(inout) :: njlu !< work array size
1159  integer(I4B), intent(inout) :: njw !< work array size
1160  integer(I4B), intent(inout) :: nwlu !< work array size
1161  ! local
1162  integer(I4B) :: n, i
1163  integer(I4B) :: ijlu, ijw, iwlu, iwk
1164 
1165  ijlu = 1
1166  ijw = 1
1167  iwlu = 1
1168 
1169  ! ILU0 and MILU0
1170  niapc = neq
1171  njapc = nja
1172 
1173  ! ILUT and MILUT
1174  if (ipc == 3 .or. ipc == 4) then
1175  niapc = neq
1176  if (level > 0) then
1177  iwk = neq * (level * 2 + 1)
1178  else
1179  iwk = 0
1180  do n = 1, neq
1181  i = ia(n + 1) - ia(n)
1182  if (i > iwk) then
1183  iwk = i
1184  end if
1185  end do
1186  iwk = neq * iwk
1187  end if
1188  njapc = iwk
1189  ijlu = iwk
1190  ijw = 2 * neq
1191  iwlu = neq + 1
1192  end if
1193 
1194  njlu = ijlu
1195  njw = ijw
1196  nwlu = iwlu
1197 
Here is the caller graph for this function:

Variable Documentation

◆ parser

type(blockparsertype), private imslinearbasemodule::parser
private

Definition at line 20 of file ImsLinearBase.f90.

20  type(BlockParserType), private :: parser