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()
35 real(dp),
pointer :: relax => null()
36 integer(I4B),
pointer :: level => null()
37 real(dp),
pointer :: droptol => null()
39 integer(I4B),
POINTER :: ipc => null()
40 integer(I4B),
POINTER :: iacpc => null()
41 integer(I4B),
POINTER :: niterc => null()
42 integer(I4B),
POINTER :: niabcgs => null()
43 integer(I4B),
POINTER :: niapc => null()
44 integer(I4B),
POINTER :: njapc => null()
45 real(dp),
POINTER :: epfact => null()
46 real(dp),
POINTER :: l2norm0 => null()
48 integer(I4B),
POINTER :: njlu => null()
49 integer(I4B),
POINTER :: njw => null()
50 integer(I4B),
POINTER :: nwlu => null()
52 integer(I4B),
POINTER :: neq => null()
53 integer(I4B),
POINTER :: nja => null()
54 integer(I4B),
dimension(:),
pointer,
contiguous :: ia => null()
55 integer(I4B),
dimension(:),
pointer,
contiguous :: ja => null()
56 real(dp),
dimension(:),
pointer,
contiguous :: amat => null()
57 real(dp),
dimension(:),
pointer,
contiguous :: rhs => null()
58 real(dp),
dimension(:),
pointer,
contiguous :: x => null()
60 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: dscale => null()
61 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: dscale2 => null()
62 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: iapc => null()
63 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: japc => null()
64 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: apc => null()
65 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: lorder => null()
66 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: iorder => null()
67 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: iaro => null()
68 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: jaro => null()
69 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: aro => null()
71 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: iw => null()
72 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: w => null()
73 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: id => null()
74 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: d => null()
75 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: p => null()
76 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: q => null()
77 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: z => null()
79 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: t => null()
80 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: v => null()
81 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: dhat => null()
82 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: phat => null()
83 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: qhat => null()
85 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: ia0 => null()
86 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: ja0 => null()
87 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: a0 => null()
89 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: jlu => null()
90 integer(I4B),
POINTER,
DIMENSION(:),
CONTIGUOUS :: jw => null()
91 real(dp),
POINTER,
DIMENSION(:),
CONTIGUOUS :: wlu => null()
112 NEQ, matrix, RHS, X, linear_settings)
120 CHARACTER(LEN=LENSOLUTIONNAME),
INTENT(IN) :: NAME
121 integer(I4B),
INTENT(IN) :: IOUT
122 integer(I4B),
TARGET,
INTENT(IN) :: IPRIMS
123 integer(I4B),
INTENT(IN) :: MXITER
124 integer(I4B),
TARGET,
INTENT(IN) :: NEQ
126 real(DP),
DIMENSION(NEQ),
TARGET,
INTENT(INOUT) :: RHS
127 real(DP),
DIMENSION(NEQ),
TARGET,
INTENT(INOUT) :: X
130 character(len=LINELENGTH) :: errmsg
133 integer(I4B) :: iscllen, iolen
140 this%DVCLOSE => linear_settings%dvclose
141 this%RCLOSE => linear_settings%rclose
142 this%ICNVGOPT => linear_settings%icnvgopt
143 this%ITER1 => linear_settings%iter1
144 this%ILINMETH => linear_settings%ilinmeth
145 this%iSCL => linear_settings%iscl
146 this%IORD => linear_settings%iord
147 this%NORTH => linear_settings%north
148 this%RELAX => linear_settings%relax
149 this%LEVEL => linear_settings%level
150 this%DROPTOL => linear_settings%droptol
153 this%IPRIMS => iprims
155 call matrix%get_aij(this%IA, this%JA, this%AMAT)
157 this%NJA =
size(this%AMAT)
162 call this%allocate_scalars()
174 02000
FORMAT(1x, /1x,
'IMSLINEAR -- UNSTRUCTURED LINEAR SOLUTION', &
175 ' PACKAGE, VERSION 8, 04/28/2017')
178 IF (this%LEVEL > 0 .OR. this%DROPTOL >
dzero)
THEN
183 IF (this%RELAX >
dzero)
THEN
184 this%IPC = this%IPC + 1
188 IF (this%ISCL < 0) this%ISCL = 0
189 IF (this%ISCL > 2)
THEN
190 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR ISCL MUST BE <= 2'
193 IF (this%IORD < 0) this%IORD = 0
194 IF (this%IORD > 2)
THEN
195 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR IORD MUST BE <= 2'
198 IF (this%NORTH < 0)
THEN
199 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR NORTH MUST >= 0'
202 IF (this%RCLOSE ==
dzero)
THEN
203 IF (this%ICNVGOPT /= 3)
THEN
204 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR RCLOSE MUST > 0.0'
208 IF (this%RELAX <
dzero)
THEN
209 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR RELAX MUST BE >= 0.0'
212 IF (this%RELAX >
done)
THEN
213 WRITE (errmsg,
'(A)')
'IMSLINEAR7AR RELAX MUST BE <= 1.0'
222 IF (this%ISCL .NE. 0) iscllen = neq
223 CALL mem_allocate(this%DSCALE, iscllen,
'DSCALE', trim(this%memoryPath))
224 CALL mem_allocate(this%DSCALE2, iscllen,
'DSCALE2', trim(this%memoryPath))
227 call ims_calc_pcdims(this%NEQ, this%NJA, this%IA, this%LEVEL, this%IPC, &
228 this%NIAPC, this%NJAPC, this%NJLU, this%NJW, this%NWLU)
231 CALL mem_allocate(this%IAPC, this%NIAPC + 1,
'IAPC', trim(this%memoryPath))
232 CALL mem_allocate(this%JAPC, this%NJAPC,
'JAPC', trim(this%memoryPath))
233 CALL mem_allocate(this%APC, this%NJAPC,
'APC', trim(this%memoryPath))
236 CALL mem_allocate(this%IW, this%NIAPC,
'IW', trim(this%memoryPath))
237 CALL mem_allocate(this%W, this%NIAPC,
'W', trim(this%memoryPath))
240 CALL mem_allocate(this%JLU, this%NJLU,
'JLU', trim(this%memoryPath))
241 CALL mem_allocate(this%JW, this%NJW,
'JW', trim(this%memoryPath))
242 CALL mem_allocate(this%WLU, this%NWLU,
'WLU', trim(this%memoryPath))
245 IF (this%IPC == 1 .OR. this%IPC == 2)
THEN
247 this%IAPC, this%JAPC)
253 IF (this%IORD .NE. 0)
THEN
257 CALL mem_allocate(this%LORDER, i0,
'LORDER', trim(this%memoryPath))
258 CALL mem_allocate(this%IORDER, i0,
'IORDER', trim(this%memoryPath))
259 CALL mem_allocate(this%IARO, i0 + 1,
'IARO', trim(this%memoryPath))
260 CALL mem_allocate(this%JARO, iolen,
'JARO', trim(this%memoryPath))
261 CALL mem_allocate(this%ARO, iolen,
'ARO', trim(this%memoryPath))
264 CALL mem_allocate(this%ID, this%NEQ,
'ID', trim(this%memoryPath))
265 CALL mem_allocate(this%D, this%NEQ,
'D', trim(this%memoryPath))
266 CALL mem_allocate(this%P, this%NEQ,
'P', trim(this%memoryPath))
267 CALL mem_allocate(this%Q, this%NEQ,
'Q', trim(this%memoryPath))
268 CALL mem_allocate(this%Z, this%NEQ,
'Z', trim(this%memoryPath))
272 IF (this%ILINMETH == 2)
THEN
273 this%NIABCGS = this%NEQ
275 CALL mem_allocate(this%T, this%NIABCGS,
'T', trim(this%memoryPath))
276 CALL mem_allocate(this%V, this%NIABCGS,
'V', trim(this%memoryPath))
277 CALL mem_allocate(this%DHAT, this%NIABCGS,
'DHAT', trim(this%memoryPath))
278 CALL mem_allocate(this%PHAT, this%NIABCGS,
'PHAT', trim(this%memoryPath))
279 CALL mem_allocate(this%QHAT, this%NIABCGS,
'QHAT', trim(this%memoryPath))
283 this%DSCALE(n) =
done
284 this%DSCALE2(n) =
done
304 DO n = 1, this%NIABCGS
333 IF (this%IORD .NE. 0)
THEN
335 this%JA, this%LORDER, this%IORDER)
352 integer(I4B),
intent(in) :: mxiter
354 CHARACTER(LEN=10) :: clin(0:2)
355 CHARACTER(LEN=31) :: clintit(0:2)
356 CHARACTER(LEN=20) :: cipc(0:4)
357 CHARACTER(LEN=20) :: cscale(0:2)
358 CHARACTER(LEN=25) :: corder(0:2)
359 CHARACTER(LEN=16),
DIMENSION(0:4) :: ccnvgopt
360 CHARACTER(LEN=15) :: clevel
361 CHARACTER(LEN=15) :: cdroptol
365 DATA clin/
'UNKNOWN ', &
368 DATA clintit/
' UNKNOWN ', &
369 &
' CONJUGATE-GRADIENT ', &
370 &
'BICONJUGATE-GRADIENT STABILIZED'/
371 DATA cipc/
'UNKNOWN ', &
373 &
'MOD. INCOMPLETE LU ', &
374 &
'INCOMPLETE LUT ', &
375 &
'MOD. INCOMPLETE LUT '/
376 DATA cscale/
'NO SCALING ', &
377 &
'SYMMETRIC SCALING ', &
379 DATA corder/
'ORIGINAL ORDERING ', &
381 &
'MINIMUM DEGREE ORDERING '/
382 DATA ccnvgopt/
'INFINITY NORM ', &
383 &
'INFINITY NORM S ', &
385 &
'RELATIVE L2NORM ', &
388 02010
FORMAT(1x, /, 7x,
'SOLUTION BY THE', 1x, a31, 1x,
'METHOD', &
390 ' MAXIMUM OF ', i0,
' CALLS OF SOLUTION ROUTINE', /, &
391 ' MAXIMUM OF ', i0, &
392 ' INTERNAL ITERATIONS PER CALL TO SOLUTION ROUTINE', /, &
393 ' LINEAR ACCELERATION METHOD =', 1x, a, /, &
394 ' MATRIX PRECONDITIONING TYPE =', 1x, a, /, &
395 ' MATRIX SCALING APPROACH =', 1x, a, /, &
396 ' MATRIX REORDERING APPROACH =', 1x, a, /, &
397 ' NUMBER OF ORTHOGONALIZATIONS =', 1x, i0, /, &
398 ' HEAD CHANGE CRITERION FOR CLOSURE =', e15.5, /, &
399 ' RESIDUAL CHANGE CRITERION FOR CLOSURE =', e15.5, /, &
400 ' RESIDUAL CONVERGENCE OPTION =', 1x, i0, /, &
401 ' RESIDUAL CONVERGENCE NORM =', 1x, a, /, &
402 ' RELAXATION FACTOR =', e15.5)
403 02015
FORMAT(
' NUMBER OF LEVELS =', a15, /, &
404 ' DROP TOLERANCE =', a15, //)
405 2030
FORMAT(1x, a20, 1x, 6(i6, 1x))
406 2040
FORMAT(1x, 20(
'-'), 1x, 6(6(
'-'), 1x))
407 2050
FORMAT(1x, 62(
'-'),/)
415 write (this%iout, 2010) &
416 clintit(this%ILINMETH), mxiter, this%ITER1, &
417 clin(this%ILINMETH), cipc(this%IPC), &
418 cscale(this%ISCL), corder(this%IORD), &
419 this%NORTH, this%DVCLOSE, this%RCLOSE, &
420 this%ICNVGOPT, ccnvgopt(this%ICNVGOPT), &
422 if (this%level > 0)
then
423 write (clevel,
'(i15)') this%level
425 if (this%droptol >
dzero)
then
426 write (cdroptol,
'(e15.5)') this%droptol
428 IF (this%level > 0 .or. this%droptol >
dzero)
THEN
429 write (this%iout, 2015) trim(adjustl(clevel)), &
430 trim(adjustl(cdroptol))
432 write (this%iout,
'(//)')
435 if (this%iord /= 0)
then
438 if (this%iprims == 2)
then
439 DO i = 1, this%neq, 6
440 write (this%iout, 2030)
'ORIGINAL NODE :', &
441 (j, j=i, min(i + 5, this%neq))
442 write (this%iout, 2040)
443 write (this%iout, 2030)
'REORDERED INDEX :', &
444 (this%lorder(j), j=i, min(i + 5, this%neq))
445 write (this%iout, 2030)
'REORDERED NODE :', &
446 (this%iorder(j), j=i, min(i + 5, this%neq))
447 write (this%iout, 2050)
471 call mem_allocate(this%niterc,
'NITERC', this%memoryPath)
472 call mem_allocate(this%niabcgs,
'NIABCGS', this%memoryPath)
475 call mem_allocate(this%epfact,
'EPFACT', this%memoryPath)
476 call mem_allocate(this%l2norm0,
'L2NORM0', this%memoryPath)
553 nullify (this%iprims)
574 integer(I4B),
INTENT(IN) :: IFDPARAM
576 SELECT CASE (ifdparam)
633 NCONV, CONVNMOD, CONVMODSTART, &
639 integer(I4B),
INTENT(INOUT) :: ICNVG
640 integer(I4B),
INTENT(IN) :: KSTP
641 integer(I4B),
INTENT(IN) :: KITER
642 integer(I4B),
INTENT(INOUT) :: IN_ITER
644 integer(I4B),
INTENT(IN) :: NCONV
645 integer(I4B),
INTENT(IN) :: CONVNMOD
646 integer(I4B),
DIMENSION(CONVNMOD + 1),
INTENT(INOUT) :: CONVMODSTART
647 character(len=31),
DIMENSION(NCONV),
INTENT(INOUT) :: CACCEL
651 integer(I4B) :: innerit
653 integer(I4B) :: itmax
660 IF (this%ISCL .NE. 0)
THEN
662 this%NEQ, this%NJA, this%IA, this%JA, &
663 this%AMAT, this%X, this%RHS, &
664 this%DSCALE, this%DSCALE2)
668 IF (this%IORD /= 0)
THEN
669 CALL dperm(this%NEQ, this%AMAT, this%JA, this%IA, &
670 this%ARO, this%JARO, this%IARO, &
671 this%LORDER, this%ID, 1)
672 CALL dvperm(this%NEQ, this%X, this%LORDER)
673 CALL dvperm(this%NEQ, this%RHS, this%LORDER)
674 this%IA0 => this%IARO
675 this%JA0 => this%JARO
684 CALL ims_base_pcu(this%iout, this%NJA, this%NEQ, this%NIAPC, this%NJAPC, &
685 this%IPC, this%RELAX, this%A0, this%IA0, this%JA0, &
686 this%APC, this%IAPC, this%JAPC, this%IW, this%W, &
687 this%LEVEL, this%DROPTOL, this%NJLU, this%NJW, &
688 this%NWLU, this%JLU, this%JW, this%WLU)
706 this%A0, this%IA0, this%JA0)
707 this%L2NORM0 = dnrm2(this%NEQ, this%D, 1)
711 IF (this%L2NORM0 ==
dzero)
THEN
717 IF (this%ILINMETH == 1)
THEN
719 this%NEQ, this%NJA, this%NIAPC, this%NJAPC, &
720 this%IPC, this%ICNVGOPT, this%NORTH, &
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, this%Z, &
725 this%NJLU, this%IW, this%JLU, &
726 nconv, convnmod, convmodstart, &
730 ELSE IF (this%ILINMETH == 2)
THEN
732 this%NEQ, this%NJA, this%NIAPC, this%NJAPC, &
733 this%IPC, this%ICNVGOPT, this%NORTH, &
734 this%ISCL, this%DSCALE, &
735 this%DVCLOSE, this%RCLOSE, this%L2NORM0, &
736 this%EPFACT, this%IA0, this%JA0, this%A0, &
737 this%IAPC, this%JAPC, this%APC, &
738 this%X, this%RHS, this%D, this%P, this%Q, &
739 this%T, this%V, this%DHAT, this%PHAT, this%QHAT, &
740 this%NJLU, this%IW, this%JLU, &
741 nconv, convnmod, convmodstart, &
746 IF (this%IORD /= 0)
THEN
747 CALL dperm(this%NEQ, this%A0, this%JA0, this%IA0, &
748 this%AMAT, this%JA, this%IA, &
749 this%IORDER, this%ID, 1)
750 CALL dvperm(this%NEQ, this%X, this%IORDER)
751 CALL dvperm(this%NEQ, this%RHS, this%IORDER)
755 IF (this%ISCL .NE. 0)
THEN
757 this%NEQ, this%NJA, this%IA, this%JA, &
758 this%AMAT, this%X, this%RHS, &
759 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.