|
2329
|
1 SUBROUTINE DDAINI (X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR, |
|
|
2 + IPAR, PHI, DELTA, E, WM, IWM, HMIN, UROUND, NONNEG, NTEMP) |
|
|
3 C***BEGIN PROLOGUE DDAINI |
|
|
4 C***SUBSIDIARY |
|
|
5 C***PURPOSE Initialization routine for DDASSL. |
|
|
6 C***LIBRARY SLATEC (DASSL) |
|
|
7 C***TYPE DOUBLE PRECISION (SDAINI-S, DDAINI-D) |
|
|
8 C***AUTHOR PETZOLD, LINDA R., (LLNL) |
|
|
9 C***DESCRIPTION |
|
|
10 C----------------------------------------------------------------- |
|
|
11 C DDAINI TAKES ONE STEP OF SIZE H OR SMALLER |
|
|
12 C WITH THE BACKWARD EULER METHOD, TO |
|
|
13 C FIND YPRIME. X AND Y ARE UPDATED TO BE CONSISTENT WITH THE |
|
|
14 C NEW STEP. A MODIFIED DAMPED NEWTON ITERATION IS USED TO |
|
|
15 C SOLVE THE CORRECTOR ITERATION. |
|
|
16 C |
|
|
17 C THE INITIAL GUESS FOR YPRIME IS USED IN THE |
|
|
18 C PREDICTION, AND IN FORMING THE ITERATION |
|
|
19 C MATRIX, BUT IS NOT INVOLVED IN THE |
|
|
20 C ERROR TEST. THIS MAY HAVE TROUBLE |
|
|
21 C CONVERGING IF THE INITIAL GUESS IS NO |
|
|
22 C GOOD, OR IF G(X,Y,YPRIME) DEPENDS |
|
|
23 C NONLINEARLY ON YPRIME. |
|
|
24 C |
|
|
25 C THE PARAMETERS REPRESENT: |
|
|
26 C X -- INDEPENDENT VARIABLE |
|
|
27 C Y -- SOLUTION VECTOR AT X |
|
|
28 C YPRIME -- DERIVATIVE OF SOLUTION VECTOR |
|
|
29 C NEQ -- NUMBER OF EQUATIONS |
|
|
30 C H -- STEPSIZE. IMDER MAY USE A STEPSIZE |
|
|
31 C SMALLER THAN H. |
|
|
32 C WT -- VECTOR OF WEIGHTS FOR ERROR |
|
|
33 C CRITERION |
|
|
34 C IDID -- COMPLETION CODE WITH THE FOLLOWING MEANINGS |
|
|
35 C IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY |
|
|
36 C IDID=-12 -- DDAINI FAILED TO FIND YPRIME |
|
|
37 C RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS |
|
|
38 C THAT ARE NOT ALTERED BY DDAINI |
|
|
39 C PHI -- WORK SPACE FOR DDAINI |
|
|
40 C DELTA,E -- WORK SPACE FOR DDAINI |
|
|
41 C WM,IWM -- REAL AND INTEGER ARRAYS STORING |
|
|
42 C MATRIX INFORMATION |
|
|
43 C |
|
|
44 C----------------------------------------------------------------- |
|
|
45 C***ROUTINES CALLED DDAJAC, DDANRM, DDASLV |
|
|
46 C***REVISION HISTORY (YYMMDD) |
|
|
47 C 830315 DATE WRITTEN |
|
|
48 C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch) |
|
|
49 C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format. |
|
|
50 C 901026 Added explicit declarations for all variables and minor |
|
|
51 C cosmetic changes to prologue. (FNF) |
|
|
52 C 901030 Minor corrections to declarations. (FNF) |
|
|
53 C***END PROLOGUE DDAINI |
|
|
54 C |
|
|
55 INTEGER NEQ, IDID, IPAR(*), IWM(*), NONNEG, NTEMP |
|
|
56 DOUBLE PRECISION |
|
|
57 * X, Y(*), YPRIME(*), H, WT(*), RPAR(*), PHI(NEQ,*), DELTA(*), |
|
|
58 * E(*), WM(*), HMIN, UROUND |
|
|
59 EXTERNAL RES, JAC |
|
|
60 C |
|
|
61 EXTERNAL DDAJAC, DDANRM, DDASLV |
|
|
62 DOUBLE PRECISION DDANRM |
|
|
63 C |
|
|
64 INTEGER I, IER, IRES, JCALC, LNJE, LNRE, M, MAXIT, MJAC, NCF, |
|
|
65 * NEF, NSF |
|
|
66 DOUBLE PRECISION |
|
|
67 * CJ, DAMP, DELNRM, ERR, OLDNRM, R, RATE, S, XOLD, YNORM |
|
|
68 LOGICAL CONVGD |
|
|
69 C |
|
|
70 PARAMETER (LNRE=12) |
|
|
71 PARAMETER (LNJE=13) |
|
|
72 C |
|
|
73 DATA MAXIT/10/,MJAC/5/ |
|
|
74 DATA DAMP/0.75D0/ |
|
|
75 C |
|
|
76 C |
|
|
77 C--------------------------------------------------- |
|
|
78 C BLOCK 1. |
|
|
79 C INITIALIZATIONS. |
|
|
80 C--------------------------------------------------- |
|
|
81 C |
|
|
82 C***FIRST EXECUTABLE STATEMENT DDAINI |
|
|
83 IDID=1 |
|
|
84 NEF=0 |
|
|
85 NCF=0 |
|
|
86 NSF=0 |
|
|
87 XOLD=X |
|
|
88 YNORM=DDANRM(NEQ,Y,WT,RPAR,IPAR) |
|
|
89 C |
|
|
90 C SAVE Y AND YPRIME IN PHI |
|
|
91 DO 100 I=1,NEQ |
|
|
92 PHI(I,1)=Y(I) |
|
|
93 100 PHI(I,2)=YPRIME(I) |
|
|
94 C |
|
|
95 C |
|
|
96 C---------------------------------------------------- |
|
|
97 C BLOCK 2. |
|
|
98 C DO ONE BACKWARD EULER STEP. |
|
|
99 C---------------------------------------------------- |
|
|
100 C |
|
|
101 C SET UP FOR START OF CORRECTOR ITERATION |
|
|
102 200 CJ=1.0D0/H |
|
|
103 X=X+H |
|
|
104 C |
|
|
105 C PREDICT SOLUTION AND DERIVATIVE |
|
|
106 DO 250 I=1,NEQ |
|
|
107 250 Y(I)=Y(I)+H*YPRIME(I) |
|
|
108 C |
|
|
109 JCALC=-1 |
|
|
110 M=0 |
|
|
111 CONVGD=.TRUE. |
|
|
112 C |
|
|
113 C |
|
|
114 C CORRECTOR LOOP. |
|
|
115 300 IWM(LNRE)=IWM(LNRE)+1 |
|
|
116 IRES=0 |
|
|
117 C |
|
|
118 CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR) |
|
|
119 IF (IRES.LT.0) GO TO 430 |
|
|
120 C |
|
|
121 C |
|
|
122 C EVALUATE THE ITERATION MATRIX |
|
|
123 IF (JCALC.NE.-1) GO TO 310 |
|
|
124 IWM(LNJE)=IWM(LNJE)+1 |
|
|
125 JCALC=0 |
|
|
126 CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H, |
|
|
127 * IER,WT,E,WM,IWM,RES,IRES, |
|
|
128 * UROUND,JAC,RPAR,IPAR,NTEMP) |
|
|
129 C |
|
|
130 S=1000000.D0 |
|
|
131 IF (IRES.LT.0) GO TO 430 |
|
|
132 IF (IER.NE.0) GO TO 430 |
|
|
133 NSF=0 |
|
|
134 C |
|
|
135 C |
|
|
136 C |
|
|
137 C MULTIPLY RESIDUAL BY DAMPING FACTOR |
|
|
138 310 CONTINUE |
|
|
139 DO 320 I=1,NEQ |
|
|
140 320 DELTA(I)=DELTA(I)*DAMP |
|
|
141 C |
|
|
142 C COMPUTE A NEW ITERATE (BACK SUBSTITUTION) |
|
|
143 C STORE THE CORRECTION IN DELTA |
|
|
144 C |
|
|
145 CALL DDASLV(NEQ,DELTA,WM,IWM) |
|
|
146 C |
|
|
147 C UPDATE Y AND YPRIME |
|
|
148 DO 330 I=1,NEQ |
|
|
149 Y(I)=Y(I)-DELTA(I) |
|
|
150 330 YPRIME(I)=YPRIME(I)-CJ*DELTA(I) |
|
|
151 C |
|
|
152 C TEST FOR CONVERGENCE OF THE ITERATION. |
|
|
153 C |
|
|
154 DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR) |
|
|
155 IF (DELNRM.LE.100.D0*UROUND*YNORM) |
|
|
156 * GO TO 400 |
|
|
157 C |
|
|
158 IF (M.GT.0) GO TO 340 |
|
|
159 OLDNRM=DELNRM |
|
|
160 GO TO 350 |
|
|
161 C |
|
|
162 340 RATE=(DELNRM/OLDNRM)**(1.0D0/M) |
|
|
163 IF (RATE.GT.0.90D0) GO TO 430 |
|
|
164 S=RATE/(1.0D0-RATE) |
|
|
165 C |
|
|
166 350 IF (S*DELNRM .LE. 0.33D0) GO TO 400 |
|
|
167 C |
|
|
168 C |
|
|
169 C THE CORRECTOR HAS NOT YET CONVERGED. UPDATE |
|
|
170 C M AND AND TEST WHETHER THE MAXIMUM |
|
|
171 C NUMBER OF ITERATIONS HAVE BEEN TRIED. |
|
|
172 C EVERY MJAC ITERATIONS, GET A NEW |
|
|
173 C ITERATION MATRIX. |
|
|
174 C |
|
|
175 M=M+1 |
|
|
176 IF (M.GE.MAXIT) GO TO 430 |
|
|
177 C |
|
|
178 IF ((M/MJAC)*MJAC.EQ.M) JCALC=-1 |
|
|
179 GO TO 300 |
|
|
180 C |
|
|
181 C |
|
|
182 C THE ITERATION HAS CONVERGED. |
|
|
183 C CHECK NONNEGATIVITY CONSTRAINTS |
|
|
184 400 IF (NONNEG.EQ.0) GO TO 450 |
|
|
185 DO 410 I=1,NEQ |
|
|
186 410 DELTA(I)=MIN(Y(I),0.0D0) |
|
|
187 C |
|
|
188 DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR) |
|
|
189 IF (DELNRM.GT.0.33D0) GO TO 430 |
|
|
190 C |
|
|
191 DO 420 I=1,NEQ |
|
|
192 Y(I)=Y(I)-DELTA(I) |
|
|
193 420 YPRIME(I)=YPRIME(I)-CJ*DELTA(I) |
|
|
194 GO TO 450 |
|
|
195 C |
|
|
196 C |
|
|
197 C EXITS FROM CORRECTOR LOOP. |
|
|
198 430 CONVGD=.FALSE. |
|
|
199 450 IF (.NOT.CONVGD) GO TO 600 |
|
|
200 C |
|
|
201 C |
|
|
202 C |
|
|
203 C----------------------------------------------------- |
|
|
204 C BLOCK 3. |
|
|
205 C THE CORRECTOR ITERATION CONVERGED. |
|
|
206 C DO ERROR TEST. |
|
|
207 C----------------------------------------------------- |
|
|
208 C |
|
|
209 DO 510 I=1,NEQ |
|
|
210 510 E(I)=Y(I)-PHI(I,1) |
|
|
211 ERR=DDANRM(NEQ,E,WT,RPAR,IPAR) |
|
|
212 C |
|
|
213 IF (ERR.LE.1.0D0) RETURN |
|
|
214 C |
|
|
215 C |
|
|
216 C |
|
|
217 C-------------------------------------------------------- |
|
|
218 C BLOCK 4. |
|
|
219 C THE BACKWARD EULER STEP FAILED. RESTORE X, Y |
|
|
220 C AND YPRIME TO THEIR ORIGINAL VALUES. |
|
|
221 C REDUCE STEPSIZE AND TRY AGAIN, IF |
|
|
222 C POSSIBLE. |
|
|
223 C--------------------------------------------------------- |
|
|
224 C |
|
|
225 600 CONTINUE |
|
|
226 X = XOLD |
|
|
227 DO 610 I=1,NEQ |
|
|
228 Y(I)=PHI(I,1) |
|
|
229 610 YPRIME(I)=PHI(I,2) |
|
|
230 C |
|
|
231 IF (CONVGD) GO TO 640 |
|
|
232 IF (IER.EQ.0) GO TO 620 |
|
|
233 NSF=NSF+1 |
|
|
234 H=H*0.25D0 |
|
|
235 IF (NSF.LT.3.AND.ABS(H).GE.HMIN) GO TO 690 |
|
|
236 IDID=-12 |
|
|
237 RETURN |
|
|
238 620 IF (IRES.GT.-2) GO TO 630 |
|
|
239 IDID=-12 |
|
|
240 RETURN |
|
|
241 630 NCF=NCF+1 |
|
|
242 H=H*0.25D0 |
|
|
243 IF (NCF.LT.10.AND.ABS(H).GE.HMIN) GO TO 690 |
|
|
244 IDID=-12 |
|
|
245 RETURN |
|
|
246 C |
|
|
247 640 NEF=NEF+1 |
|
|
248 R=0.90D0/(2.0D0*ERR+0.0001D0) |
|
|
249 R=MAX(0.1D0,MIN(0.5D0,R)) |
|
|
250 H=H*R |
|
|
251 IF (ABS(H).GE.HMIN.AND.NEF.LT.10) GO TO 690 |
|
|
252 IDID=-12 |
|
|
253 RETURN |
|
|
254 690 GO TO 200 |
|
|
255 C |
|
|
256 C-------------END OF SUBROUTINE DDAINI---------------------- |
|
|
257 END |
