|
3911
|
1 C Work performed under the auspices of the U.S. Department of Energy |
|
|
2 C by Lawrence Livermore National Laboratory under contract number |
|
|
3 C W-7405-Eng-48. |
|
|
4 C |
|
|
5 SUBROUTINE DSLVK (NEQ, Y, TN, YPRIME, SAVR, X, EWT, WM, IWM, |
|
|
6 * RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK, |
|
|
7 * RPAR, IPAR) |
|
|
8 C |
|
|
9 C***BEGIN PROLOGUE DSLVK |
|
|
10 C***REFER TO DDASPK |
|
|
11 C***DATE WRITTEN 890101 (YYMMDD) |
|
|
12 C***REVISION DATE 900926 (YYMMDD) |
|
|
13 C***REVISION DATE 940928 Removed MNEWT and added RHOK in call list. |
|
|
14 C |
|
|
15 C |
|
|
16 C----------------------------------------------------------------------- |
|
|
17 C***DESCRIPTION |
|
|
18 C |
|
|
19 C DSLVK uses a restart algorithm and interfaces to DSPIGM for |
|
|
20 C the solution of the linear system arising from a Newton iteration. |
|
|
21 C |
|
|
22 C In addition to variables described elsewhere, |
|
|
23 C communication with DSLVK uses the following variables.. |
|
|
24 C WM = Real work space containing data for the algorithm |
|
|
25 C (Krylov basis vectors, Hessenberg matrix, etc.). |
|
|
26 C IWM = Integer work space containing data for the algorithm. |
|
|
27 C X = The right-hand side vector on input, and the solution vector |
|
|
28 C on output, of length NEQ. |
|
|
29 C IRES = Error flag from RES. |
|
|
30 C IERSL = Output flag .. |
|
|
31 C IERSL = 0 means no trouble occurred (or user RES routine |
|
|
32 C returned IRES < 0) |
|
|
33 C IERSL = 1 means the iterative method failed to converge |
|
|
34 C (DSPIGM returned IFLAG > 0.) |
|
|
35 C IERSL = -1 means there was a nonrecoverable error in the |
|
|
36 C iterative solver, and an error exit will occur. |
|
|
37 C----------------------------------------------------------------------- |
|
|
38 C***ROUTINES CALLED |
|
|
39 C DSCAL, DCOPY, DSPIGM |
|
|
40 C |
|
|
41 C***END PROLOGUE DSLVK |
|
|
42 C |
|
|
43 INTEGER NEQ, IWM, IRES, IERSL, IPAR |
|
|
44 DOUBLE PRECISION Y, TN, YPRIME, SAVR, X, EWT, WM, CJ, EPLIN, |
|
|
45 1 SQRTN, RSQRTN, RHOK, RPAR |
|
|
46 DIMENSION Y(*), YPRIME(*), SAVR(*), X(*), EWT(*), |
|
|
47 1 WM(*), IWM(*), RPAR(*), IPAR(*) |
|
|
48 C |
|
|
49 INTEGER IFLAG, IRST, NRSTS, NRMAX, LR, LDL, LHES, LGMR, LQ, LV, |
|
|
50 1 LWK, LZ, MAXLP1, NPSL |
|
|
51 INTEGER NLI, NPS, NCFL, NRE, MAXL, KMP, MITER |
|
|
52 EXTERNAL RES, PSOL |
|
|
53 C |
|
|
54 PARAMETER (LNRE=12, LNCFL=16, LNLI=20, LNPS=21) |
|
|
55 PARAMETER (LLOCWP=29, LLCIWP=30) |
|
|
56 PARAMETER (LMITER=23, LMAXL=24, LKMP=25, LNRMAX=26) |
|
|
57 C |
|
|
58 C----------------------------------------------------------------------- |
|
|
59 C IRST is set to 1, to indicate restarting is in effect. |
|
|
60 C NRMAX is the maximum number of restarts. |
|
|
61 C----------------------------------------------------------------------- |
|
|
62 DATA IRST/1/ |
|
|
63 C |
|
|
64 LIWP = IWM(LLCIWP) |
|
|
65 NLI = IWM(LNLI) |
|
|
66 NPS = IWM(LNPS) |
|
|
67 NCFL = IWM(LNCFL) |
|
|
68 NRE = IWM(LNRE) |
|
|
69 LWP = IWM(LLOCWP) |
|
|
70 MAXL = IWM(LMAXL) |
|
|
71 KMP = IWM(LKMP) |
|
|
72 NRMAX = IWM(LNRMAX) |
|
|
73 MITER = IWM(LMITER) |
|
|
74 IERSL = 0 |
|
|
75 IRES = 0 |
|
|
76 C----------------------------------------------------------------------- |
|
|
77 C Use a restarting strategy to solve the linear system |
|
|
78 C P*X = -F. Parse the work vector, and perform initializations. |
|
|
79 C Note that zero is the initial guess for X. |
|
|
80 C----------------------------------------------------------------------- |
|
|
81 MAXLP1 = MAXL + 1 |
|
|
82 LV = 1 |
|
|
83 LR = LV + NEQ*MAXL |
|
|
84 LHES = LR + NEQ + 1 |
|
|
85 LQ = LHES + MAXL*MAXLP1 |
|
|
86 LWK = LQ + 2*MAXL |
|
|
87 LDL = LWK + MIN0(1,MAXL-KMP)*NEQ |
|
|
88 LZ = LDL + NEQ |
|
|
89 CALL DSCAL (NEQ, RSQRTN, EWT, 1) |
|
|
90 CALL DCOPY (NEQ, X, 1, WM(LR), 1) |
|
|
91 DO 110 I = 1,NEQ |
|
|
92 110 X(I) = 0.D0 |
|
|
93 C----------------------------------------------------------------------- |
|
|
94 C Top of loop for the restart algorithm. Initial pass approximates |
|
|
95 C X and sets up a transformed system to perform subsequent restarts |
|
|
96 C to update X. NRSTS is initialized to -1, because restarting |
|
|
97 C does not occur until after the first pass. |
|
|
98 C Update NRSTS; conditionally copy DL to R; call the DSPIGM |
|
|
99 C algorithm to solve A*Z = R; updated counters; update X with |
|
|
100 C the residual solution. |
|
|
101 C Note: if convergence is not achieved after NRMAX restarts, |
|
|
102 C then the linear solver is considered to have failed. |
|
|
103 C----------------------------------------------------------------------- |
|
|
104 NRSTS = -1 |
|
|
105 115 CONTINUE |
|
|
106 NRSTS = NRSTS + 1 |
|
|
107 IF (NRSTS .GT. 0) CALL DCOPY (NEQ, WM(LDL), 1, WM(LR),1) |
|
|
108 CALL DSPIGM (NEQ, TN, Y, YPRIME, SAVR, WM(LR), EWT, MAXL, MAXLP1, |
|
|
109 1 KMP, EPLIN, CJ, RES, IRES, NRES, PSOL, NPSL, WM(LZ), WM(LV), |
|
|
110 2 WM(LHES), WM(LQ), LGMR, WM(LWP), IWM(LIWP), WM(LWK), |
|
|
111 3 WM(LDL), RHOK, IFLAG, IRST, NRSTS, RPAR, IPAR) |
|
|
112 NLI = NLI + LGMR |
|
|
113 NPS = NPS + NPSL |
|
|
114 NRE = NRE + NRES |
|
|
115 DO 120 I = 1,NEQ |
|
|
116 120 X(I) = X(I) + WM(LZ+I-1) |
|
|
117 IF ((IFLAG .EQ. 1) .AND. (NRSTS .LT. NRMAX) .AND. (IRES .EQ. 0)) |
|
|
118 1 GO TO 115 |
|
|
119 C----------------------------------------------------------------------- |
|
|
120 C The restart scheme is finished. Test IRES and IFLAG to see if |
|
|
121 C convergence was not achieved, and set flags accordingly. |
|
|
122 C----------------------------------------------------------------------- |
|
|
123 IF (IRES .LT. 0) THEN |
|
|
124 NCFL = NCFL + 1 |
|
|
125 ELSE IF (IFLAG .NE. 0) THEN |
|
|
126 NCFL = NCFL + 1 |
|
|
127 IF (IFLAG .GT. 0) IERSL = 1 |
|
|
128 IF (IFLAG .LT. 0) IERSL = -1 |
|
|
129 ENDIF |
|
|
130 C----------------------------------------------------------------------- |
|
|
131 C Update IWM with counters, rescale EWT, and return. |
|
|
132 C----------------------------------------------------------------------- |
|
|
133 IWM(LNLI) = NLI |
|
|
134 IWM(LNPS) = NPS |
|
|
135 IWM(LNCFL) = NCFL |
|
|
136 IWM(LNRE) = NRE |
|
|
137 CALL DSCAL (NEQ, SQRTN, EWT, 1) |
|
|
138 RETURN |
|
|
139 C |
|
|
140 C------END OF SUBROUTINE DSLVK------------------------------------------ |
|
|
141 END |
