23 character(len=LENMEMPATH) :: memorypath
24 integer(I4B),
POINTER :: iout => null()
25 integer(I4B),
POINTER :: iprims => null()
27 real(dp),
pointer :: dvclose => null()
28 real(dp),
pointer :: rclose => null()
29 integer(I4B),
pointer :: icnvgopt => null()
30 integer(I4B),
pointer :: iter1 => null()
31 integer(I4B),
pointer :: ilinmeth => null()
32 integer(I4B),
pointer :: iscl => null()
33 integer(I4B),
pointer :: iord => null()
34 integer(I4B),
pointer :: north => null()
36 real(dp),
pointer :: relax => null()
37 integer(I4B),
pointer :: level => null()
38 real(dp),
pointer :: droptol => null()
40 integer(I4B),
POINTER :: ipc => null()
41 integer(I4B),
POINTER :: iacpc => null()
42 integer(I4B),
POINTER :: niterc => null()
43 integer(I4B),
POINTER :: niabcgs => null()
44 integer(I4B),
POINTER :: niapc => null()
45 integer(I4B),
POINTER :: njapc => null()
46 real(dp),
POINTER :: epfact => null()
47 real(dp),
POINTER :: l2norm0 => null()
49 integer(I4B),
POINTER :: njlu => null()
50 integer(I4B),
POINTER :: njw => null()
51 integer(I4B),
POINTER :: nwlu => null()
53 integer(I4B),
POINTER :: neq => null()
54 integer(I4B),
POINTER :: nja => null()
55 integer(I4B),
dimension(:),
pointer,
contiguous :: ia => null()
56 integer(I4B),
dimension(:),
pointer,
contiguous :: ja => null()
57 real(dp),
dimension(:),
pointer,
contiguous :: amat => null()
58 real(dp),
dimension(:),
pointer,
contiguous :: rhs => null()
59 real(dp),
dimension(:),
pointer,
contiguous :: x => null()
61 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: dscale => null()
62 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: dscale2 => null()
63 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: iapc => null()
64 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: japc => null()
65 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: apc => null()
66 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: lorder => null()
67 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: iorder => null()
68 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: iaro => null()
69 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: jaro => null()
70 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: aro => null()
72 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: iw => null()
73 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: w => null()
74 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: id => null()
75 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: d => null()
76 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: p => null()
77 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: q => null()
78 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: z => null()
80 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: t => null()
81 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: v => null()
82 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: dhat => null()
83 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: phat => null()
84 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: qhat => null()
86 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: ia0 => null()
87 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: ja0 => null()
88 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: a0 => null()
90 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: jlu => null()
91 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: jw => null()
92 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: wlu => null()
113 NEQ, matrix, RHS, X, linear_settings)
121 CHARACTER(LEN=LENSOLUTIONNAME),
INTENT(IN) :: NAME
122 integer(I4B),
INTENT(IN) :: IOUT
123 integer(I4B),
TARGET,
INTENT(IN) :: IPRIMS
124 integer(I4B),
INTENT(IN) :: MXITER
125 integer(I4B),
TARGET,
INTENT(IN) :: NEQ
127 real(DP),
DIMENSION(NEQ),
TARGET,
INTENT(INOUT) :: RHS
128 real(DP),
DIMENSION(NEQ),
TARGET,
INTENT(INOUT) :: X
131 character(len=LINELENGTH) :: errmsg
134 integer(I4B) :: iscllen, iolen
141 this%DVCLOSE => linear_settings%dvclose
142 this%RCLOSE => linear_settings%rclose
143 this%ICNVGOPT => linear_settings%icnvgopt
144 this%ITER1 => linear_settings%iter1
145 this%ILINMETH => linear_settings%ilinmeth
146 this%ISCL => linear_settings%iscl
147 this%IORD => linear_settings%iord
148 this%NORTH => linear_settings%north
149 this%RELAX => linear_settings%relax
150 this%LEVEL => linear_settings%level
151 this%DROPTOL => linear_settings%droptol
154 this%IPRIMS => iprims
156 call matrix%get_aij(this%IA, this%JA, this%AMAT)
158 this%NJA =
size(this%AMAT)
163 call this%allocate_scalars()
175 02000
FORMAT(1x, /1x,
'IMSLINEAR -- UNSTRUCTURED LINEAR SOLUTION', &
176 ' PACKAGE, VERSION 8, 04/28/2017')
179 IF (this%LEVEL > 0 .OR. this%DROPTOL >
dzero)
THEN
184 IF (this%RELAX >
dzero)
THEN
185 this%IPC = this%IPC + 1
189 IF (this%ISCL < 0) this%ISCL = 0
190 IF (this%ISCL > 2)
THEN
191 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR ISCL MUST BE <= 2'
194 IF (this%IORD < 0) this%IORD = 0
195 IF (this%IORD > 2)
THEN
196 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR IORD MUST BE <= 2'
199 IF (this%NORTH < 0)
THEN
200 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR NORTH MUST >= 0'
203 IF (this%RCLOSE ==
dzero)
THEN
204 IF (this%ICNVGOPT /= 3)
THEN
205 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR RCLOSE MUST > 0.0'
209 IF (this%RELAX <
dzero)
THEN
210 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR RELAX MUST BE >= 0.0'
213 IF (this%RELAX >
done)
THEN
214 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR RELAX MUST BE <= 1.0'
223 IF (this%ISCL .NE. 0) iscllen = neq
224 CALL mem_allocate(this%DSCALE, iscllen,
'DSCALE', trim(this%memoryPath))
225 CALL mem_allocate(this%DSCALE2, iscllen,
'DSCALE2', trim(this%memoryPath))
228 call ims_calc_pcdims(this%NEQ, this%NJA, this%IA, this%LEVEL, this%IPC, &
229 this%NIAPC, this%NJAPC, this%NJLU, this%NJW, this%NWLU)
232 CALL mem_allocate(this%IAPC, this%NIAPC + 1,
'IAPC', trim(this%memoryPath))
233 CALL mem_allocate(this%JAPC, this%NJAPC,
'JAPC', trim(this%memoryPath))
234 CALL mem_allocate(this%APC, this%NJAPC,
'APC', trim(this%memoryPath))
237 CALL mem_allocate(this%IW, this%NIAPC,
'IW', trim(this%memoryPath))
238 CALL mem_allocate(this%W, this%NIAPC,
'W', trim(this%memoryPath))
241 CALL mem_allocate(this%JLU, this%NJLU,
'JLU', trim(this%memoryPath))
242 CALL mem_allocate(this%JW, this%NJW,
'JW', trim(this%memoryPath))
243 CALL mem_allocate(this%WLU, this%NWLU,
'WLU', trim(this%memoryPath))
246 IF (this%IPC == 1 .OR. this%IPC == 2)
THEN
248 this%IAPC, this%JAPC)
254 IF (this%IORD .NE. 0)
THEN
258 CALL mem_allocate(this%LORDER, i0,
'LORDER', trim(this%memoryPath))
259 CALL mem_allocate(this%IORDER, i0,
'IORDER', trim(this%memoryPath))
260 CALL mem_allocate(this%IARO, i0 + 1,
'IARO', trim(this%memoryPath))
261 CALL mem_allocate(this%JARO, iolen,
'JARO', trim(this%memoryPath))
262 CALL mem_allocate(this%ARO, iolen,
'ARO', trim(this%memoryPath))
265 CALL mem_allocate(this%ID, this%NEQ,
'ID', trim(this%memoryPath))
266 CALL mem_allocate(this%D, this%NEQ,
'D', trim(this%memoryPath))
267 CALL mem_allocate(this%P, this%NEQ,
'P', trim(this%memoryPath))
268 CALL mem_allocate(this%Q, this%NEQ,
'Q', trim(this%memoryPath))
269 CALL mem_allocate(this%Z, this%NEQ,
'Z', trim(this%memoryPath))
273 IF (this%ILINMETH == 2)
THEN
274 this%NIABCGS = this%NEQ
276 CALL mem_allocate(this%T, this%NIABCGS,
'T', trim(this%memoryPath))
277 CALL mem_allocate(this%V, this%NIABCGS,
'V', trim(this%memoryPath))
278 CALL mem_allocate(this%DHAT, this%NIABCGS,
'DHAT', trim(this%memoryPath))
279 CALL mem_allocate(this%PHAT, this%NIABCGS,
'PHAT', trim(this%memoryPath))
280 CALL mem_allocate(this%QHAT, this%NIABCGS,
'QHAT', trim(this%memoryPath))
284 this%DSCALE(n) =
done
285 this%DSCALE2(n) =
done
305 DO n = 1, this%NIABCGS
334 IF (this%IORD .NE. 0)
THEN
336 this%JA, this%LORDER, this%IORDER)
350 integer(I4B),
intent(in) :: mxiter
352 CHARACTER(LEN=10) :: clin(0:2)
353 CHARACTER(LEN=31) :: clintit(0:2)
354 CHARACTER(LEN=20) :: cipc(0:4)
355 CHARACTER(LEN=20) :: cscale(0:2)
356 CHARACTER(LEN=25) :: corder(0:2)
357 CHARACTER(LEN=16),
DIMENSION(0:4) :: ccnvgopt
358 CHARACTER(LEN=15) :: clevel
359 CHARACTER(LEN=15) :: cdroptol
363 DATA clin/
'UNKNOWN ', &
366 DATA clintit/
' UNKNOWN ', &
367 &
' CONJUGATE-GRADIENT ', &
368 &
'BICONJUGATE-GRADIENT STABILIZED'/
369 DATA cipc/
'UNKNOWN ', &
371 &
'MOD. INCOMPLETE LU ', &
372 &
'INCOMPLETE LUT ', &
373 &
'MOD. INCOMPLETE LUT '/
374 DATA cscale/
'NO SCALING ', &
375 &
'SYMMETRIC SCALING ', &
377 DATA corder/
'ORIGINAL ORDERING ', &
379 &
'MINIMUM DEGREE ORDERING '/
380 DATA ccnvgopt/
'INFINITY NORM ', &
381 &
'INFINITY NORM S ', &
383 &
'RELATIVE L2NORM ', &
386 02010
FORMAT(1x, /, 7x,
'SOLUTION BY THE', 1x, a31, 1x,
'METHOD', &
388 ' MAXIMUM OF ', i0,
' CALLS OF SOLUTION ROUTINE', /, &
389 ' MAXIMUM OF ', i0, &
390 ' INTERNAL ITERATIONS PER CALL TO SOLUTION ROUTINE', /, &
391 ' LINEAR ACCELERATION METHOD =', 1x, a, /, &
392 ' MATRIX PRECONDITIONING TYPE =', 1x, a, /, &
393 ' MATRIX SCALING APPROACH =', 1x, a, /, &
394 ' MATRIX REORDERING APPROACH =', 1x, a, /, &
395 ' NUMBER OF ORTHOGONALIZATIONS =', 1x, i0, /, &
396 ' HEAD CHANGE CRITERION FOR CLOSURE =', e15.5, /, &
397 ' RESIDUAL CHANGE CRITERION FOR CLOSURE =', e15.5, /, &
398 ' RESIDUAL CONVERGENCE OPTION =', 1x, i0, /, &
399 ' RESIDUAL CONVERGENCE NORM =', 1x, a, /, &
400 ' RELAXATION FACTOR =', e15.5)
401 02015
FORMAT(
' NUMBER OF LEVELS =', a15, /, &
402 ' DROP TOLERANCE =', a15, //)
403 2030
FORMAT(1x, a20, 1x, 6(i6, 1x))
404 2040
FORMAT(1x, 20(
'-'), 1x, 6(6(
'-'), 1x))
405 2050
FORMAT(1x, 62(
'-'),/)
413 write (this%iout, 2010) &
414 clintit(this%ILINMETH), mxiter, this%ITER1, &
415 clin(this%ILINMETH), cipc(this%IPC), &
416 cscale(this%ISCL), corder(this%IORD), &
417 this%NORTH, this%DVCLOSE, this%RCLOSE, &
418 this%ICNVGOPT, ccnvgopt(this%ICNVGOPT), &
420 if (this%level > 0)
then
421 write (clevel,
'(i15)') this%level
423 if (this%droptol >
dzero)
then
424 write (cdroptol,
'(e15.5)') this%droptol
426 IF (this%level > 0 .or. this%droptol >
dzero)
THEN
427 write (this%iout, 2015) trim(adjustl(clevel)), &
428 trim(adjustl(cdroptol))
430 write (this%iout,
'(//)')
433 if (this%iord /= 0)
then
436 if (this%iprims == 2)
then
437 DO i = 1, this%neq, 6
438 write (this%iout, 2030)
'ORIGINAL NODE :', &
439 (j, j=i, min(i + 5, this%neq))
440 write (this%iout, 2040)
441 write (this%iout, 2030)
'REORDERED INDEX :', &
442 (this%lorder(j), j=i, min(i + 5, this%neq))
443 write (this%iout, 2030)
'REORDERED NODE :', &
444 (this%iorder(j), j=i, min(i + 5, this%neq))
445 write (this%iout, 2050)
466 call mem_allocate(this%niterc,
'NITERC', this%memoryPath)
467 call mem_allocate(this%niabcgs,
'NIABCGS', this%memoryPath)
470 call mem_allocate(this%epfact,
'EPFACT', this%memoryPath)
471 call mem_allocate(this%l2norm0,
'L2NORM0', this%memoryPath)
545 nullify (this%iprims)
563 integer(I4B),
INTENT(IN) :: IFDPARAM
565 SELECT CASE (ifdparam)
619 NCONV, CONVNMOD, CONVMODSTART, &
625 integer(I4B),
INTENT(INOUT) :: ICNVG
626 integer(I4B),
INTENT(IN) :: KSTP
627 integer(I4B),
INTENT(IN) :: KITER
628 integer(I4B),
INTENT(INOUT) :: IN_ITER
630 integer(I4B),
INTENT(IN) :: NCONV
631 integer(I4B),
INTENT(IN) :: CONVNMOD
632 integer(I4B),
DIMENSION(CONVNMOD + 1),
INTENT(INOUT) :: CONVMODSTART
633 character(len=31),
DIMENSION(NCONV),
INTENT(INOUT) :: CACCEL
637 integer(I4B) :: innerit
639 integer(I4B) :: itmax
646 IF (this%ISCL .NE. 0)
THEN
648 this%NEQ, this%NJA, this%IA, this%JA, &
649 this%AMAT, this%X, this%RHS, &
650 this%DSCALE, this%DSCALE2)
654 IF (this%IORD /= 0)
THEN
655 CALL dperm(this%NEQ, this%AMAT, this%JA, this%IA, &
656 this%ARO, this%JARO, this%IARO, &
657 this%LORDER, this%ID, 1)
658 CALL dvperm(this%NEQ, this%X, this%LORDER)
659 CALL dvperm(this%NEQ, this%RHS, this%LORDER)
660 this%IA0 => this%IARO
661 this%JA0 => this%JARO
670 CALL ims_base_pcu(this%iout, this%NJA, this%NEQ, this%NIAPC, this%NJAPC, &
671 this%IPC, this%RELAX, this%A0, this%IA0, this%JA0, &
672 this%APC, this%IAPC, this%JAPC, this%IW, this%W, &
673 this%LEVEL, this%DROPTOL, this%NJLU, this%NJW, &
674 this%NWLU, this%JLU, this%JW, this%WLU)
692 this%A0, this%IA0, this%JA0)
693 this%L2NORM0 = dnrm2(this%NEQ, this%D, 1)
697 IF (this%L2NORM0 ==
dzero)
THEN
703 IF (this%ILINMETH == 1)
THEN
705 this%NEQ, this%NJA, this%NIAPC, this%NJAPC, &
706 this%IPC, this%ICNVGOPT, this%NORTH, &
707 this%DVCLOSE, this%RCLOSE, this%L2NORM0, &
708 this%EPFACT, this%IA0, this%JA0, this%A0, &
709 this%IAPC, this%JAPC, this%APC, &
710 this%X, this%RHS, this%D, this%P, this%Q, this%Z, &
711 this%NJLU, this%IW, this%JLU, &
712 nconv, convnmod, convmodstart, &
716 ELSE IF (this%ILINMETH == 2)
THEN
718 this%NEQ, this%NJA, this%NIAPC, this%NJAPC, &
719 this%IPC, this%ICNVGOPT, this%NORTH, &
720 this%ISCL, this%DSCALE, &
721 this%DVCLOSE, this%RCLOSE, this%L2NORM0, &
722 this%EPFACT, this%IA0, this%JA0, this%A0, &
723 this%IAPC, this%JAPC, this%APC, &
724 this%X, this%RHS, this%D, this%P, this%Q, &
725 this%T, this%V, this%DHAT, this%PHAT, this%QHAT, &
726 this%NJLU, this%IW, this%JLU, &
727 nconv, convnmod, convmodstart, &
732 IF (this%IORD /= 0)
THEN
733 CALL dperm(this%NEQ, this%A0, this%JA0, this%IA0, &
734 this%AMAT, this%JA, this%IA, &
735 this%IORDER, this%ID, 1)
736 CALL dvperm(this%NEQ, this%X, this%IORDER)
737 CALL dvperm(this%NEQ, this%RHS, this%IORDER)
741 IF (this%ISCL .NE. 0)
THEN
743 this%NEQ, this%NJA, this%IA, this%JA, &
744 this%AMAT, this%X, this%RHS, &
745 this%DSCALE, this%DSCALE2)
This module contains block parser methods.
This module contains simulation constants.
real(dp), parameter dsame
real constant for values that are considered the same based on machine precision
integer(i4b), parameter linelength
maximum length of a standard line
integer(i4b), parameter lensolutionname
maximum length of the solution name
real(dp), parameter dem8
real constant 1e-8
real(dp), parameter dem1
real constant 1e-1
real(dp), parameter dhalf
real constant 1/2
real(dp), parameter dem3
real constant 1e-3
integer(i4b), parameter izero
integer constant zero
real(dp), parameter dem4
real constant 1e-4
real(dp), parameter dem6
real constant 1e-6
real(dp), parameter dzero
real constant zero
real(dp), parameter dem5
real constant 1e-5
real(dp), parameter dprec
real constant machine precision
@ vdebug
write debug output
real(dp), parameter dem2
real constant 1e-2
real(dp), parameter dtwo
real constant 2
integer(i4b), parameter lenmempath
maximum length of the memory path
real(dp), parameter done
real constant 1
This module contains the IMS linear accelerator subroutines.
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
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
subroutine ims_base_pccrs(NEQ, NJA, IA, JA, IAPC, JAPC)
@ brief Generate CRS pointers for the preconditioner
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
subroutine ims_calc_pcdims(neq, nja, ia, level, ipc, niapc, njapc, njlu, njw, nwlu)
subroutine ims_base_scale(IOPT, ISCL, NEQ, NJA, IA, JA, AMAT, X, B, DSCALE, DSCALE2)
@ brief Scale the coefficient matrix
real(dp) function ims_base_epfact(icnvgopt, kstp)
Function returning EPFACT.
subroutine ims_base_calc_order(IORD, NEQ, NJA, IA, JA, LORDER, IORDER)
@ brief Calculate LORDER AND IORDER
subroutine ims_base_residual(NEQ, NJA, X, B, D, A, IA, JA)
Calculate residual.
subroutine allocate_scalars(this)
@ brief Allocate and initialize scalars
subroutine imslinear_summary(this, mxiter)
@ brief Write summary of settings
subroutine imslinear_set_input(this, IFDPARAM)
@ brief Set default settings
subroutine imslinear_da(this)
@ brief Deallocate memory
subroutine imslinear_ar(this, NAME, IOUT, IPRIMS, MXITER, NEQ, matrix, RHS, X, linear_settings)
@ brief Allocate storage and read data
subroutine imslinear_ap(this, ICNVG, KSTP, KITER, IN_ITER, NCONV, CONVNMOD, CONVMODSTART, CACCEL, summary)
@ brief Base linear accelerator subroutine
This module defines variable data types.
character(len=lenmempath) function create_mem_path(component, subcomponent, context)
returns the path to the memory object
This module contains simulation methods.
subroutine, public store_error(msg, terminate)
Store an error message.
integer(i4b) function, public count_errors()
Return number of errors.
subroutine, public deprecation_warning(cblock, cvar, cver, endmsg, iunit)
Store deprecation warning message.
subroutine dvperm(n, x, perm)
subroutine dperm(nrow, a, ja, ia, ao, jao, iao, perm, qperm, job)
This structure stores the generic convergence info for a solution.