1 | SUBROUTINE INTEGRATE( TIN, TOUT ) |
---|
2 | |
---|
3 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
4 | INCLUDE 'KPP_ROOT_Parameters.h' |
---|
5 | INCLUDE 'KPP_ROOT_Global.h' |
---|
6 | |
---|
7 | C TIN - Start Time |
---|
8 | KPP_REAL TIN |
---|
9 | C TOUT - End Time |
---|
10 | KPP_REAL TOUT, H |
---|
11 | SAVE H |
---|
12 | DATA H /0.0d0/ |
---|
13 | INTEGER Nfun, Njac, Nstp, Nacc, Nrej, Ndec, Nsol |
---|
14 | SAVE Nstp, Nacc, Nrej |
---|
15 | DATA Nstp /0/ |
---|
16 | DATA Nacc /0/ |
---|
17 | DATA Nrej /0/ |
---|
18 | INTEGER i |
---|
19 | |
---|
20 | PARAMETER (LWORK=5*NVAR*NVAR+12*NVAR+20,LIWORK=3*NVAR+20) |
---|
21 | PARAMETER (LRCONT=4*NVAR+4) |
---|
22 | |
---|
23 | KPP_REAL WORK(LWORK), RPAR(1) |
---|
24 | INTEGER IWORK(LIWORK), IPAR(1) |
---|
25 | COMMON /CONT/ ICONT(4),RCONT(LRCONT) |
---|
26 | EXTERNAL FUNC_CHEM,JAC_CHEM,SOLOUT |
---|
27 | |
---|
28 | |
---|
29 | ITOL=1 ! --- VECTOR TOLERANCES |
---|
30 | IJAC=1 ! --- COMPUTE THE JACOBIAN ANALYTICALLY |
---|
31 | IMAS=0 ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM |
---|
32 | IOUT=0 ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION |
---|
33 | MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX |
---|
34 | MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX |
---|
35 | MLMAS=NVAR ! --- JACOBIAN IS A FULL MATRIX |
---|
36 | MUMAS=NVAR ! --- JACOBIAN IS A FULL MATRIX |
---|
37 | |
---|
38 | DO i=1,20 |
---|
39 | IWORK(i) = 0 |
---|
40 | WORK(i) = 0.D0 |
---|
41 | ENDDO |
---|
42 | |
---|
43 | IWORK(3) = 8 ! Max no. of Newton iterations |
---|
44 | IWORK(4) = 0 ! Starting values for Newton are interpolated (0) or zero (1) |
---|
45 | IWORK(8) = 2 ! Gustaffson (1) or classic(2) controller |
---|
46 | WORK(2) = 0.9 ! Safety factor |
---|
47 | |
---|
48 | CALL RADAU5(NVAR,FUNC_CHEM,TIN,VAR,TOUT,H, |
---|
49 | & RTOL,ATOL,ITOL, |
---|
50 | & JAC_CHEM ,IJAC,MLJAC,MUJAC, |
---|
51 | & FUNC_CHEM ,IMAS,MLMAS,MUMAS, |
---|
52 | & SOLOUT,IOUT, |
---|
53 | & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) |
---|
54 | |
---|
55 | |
---|
56 | Nfun = Nfun + IWORK(14) |
---|
57 | Njac = Njac + IWORK(15) |
---|
58 | Nstp = Nstp + IWORK(16) |
---|
59 | Nacc = Nacc + IWORK(17) |
---|
60 | Nrej = Nrej + IWORK(18) |
---|
61 | Ndec = Ndec + IWORK(19) |
---|
62 | Nsol = Nsol + IWORK(20) |
---|
63 | |
---|
64 | print("('Nstep=',I5,' Nacc=',I6,' Nrej=',I6)"),Nstp, Nacc, Nrej |
---|
65 | |
---|
66 | IF (IDID.LT.0) THEN |
---|
67 | print *,'RADAU: Unsucessfull exit at T=', |
---|
68 | & TIN,' (IDID=',IDID,')' |
---|
69 | ENDIF |
---|
70 | |
---|
71 | RETURN |
---|
72 | END |
---|
73 | |
---|
74 | |
---|
75 | SUBROUTINE RADAU5(N,FCN,X,Y,XEND,H, |
---|
76 | & RelTol,AbsTol,ITOL, |
---|
77 | & JAC ,IJAC,MLJAC,MUJAC, |
---|
78 | & MAS ,IMAS,MLMAS,MUMAS, |
---|
79 | & SOLOUT,IOUT, |
---|
80 | & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) |
---|
81 | C ---------------------------------------------------------- |
---|
82 | C NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) |
---|
83 | C SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS |
---|
84 | C M*Y'=F(X,Y). |
---|
85 | C THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M .NE. I) |
---|
86 | C OR EXPLICIT (M=I). |
---|
87 | C THE METHOD USED IS AN IMPLICIT RUNGE-KUTTA METHOD (RADAU IIA) |
---|
88 | C OF ORDER 5 WITH STEP SIZE CONTROL AND CONTINUOUS OUTPUT. |
---|
89 | C C.F. SECTION IV.8 |
---|
90 | C |
---|
91 | C AUTHORS: E. HAIRER AND G. WANNER |
---|
92 | C UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES |
---|
93 | C CH-1211 GENEVE 24, SWITZERLAND |
---|
94 | C E-MAIL: HAIRER@DIVSUN.UNIGE.CH, WANNER@DIVSUN.UNIGE.CH |
---|
95 | C |
---|
96 | C THIS CODE IS PART OF THE BOOK: |
---|
97 | C E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL |
---|
98 | C EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. |
---|
99 | C SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14, |
---|
100 | C SPRINGER-VERLAG (1991) |
---|
101 | C |
---|
102 | C VERSION OF SEPTEMBER 30, 1995 |
---|
103 | C |
---|
104 | C INPUT PARAMETERS |
---|
105 | C ---------------- |
---|
106 | C N DIMENSION OF THE SYSTEM |
---|
107 | C |
---|
108 | C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE |
---|
109 | C VALUE OF F(X,Y): |
---|
110 | C SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR) |
---|
111 | C KPP_REAL X,Y(N),F(N) |
---|
112 | C F(1)=... ETC. |
---|
113 | C RPAR, IPAR (SEE BELOW) |
---|
114 | C |
---|
115 | C X INITIAL X-VALUE |
---|
116 | C |
---|
117 | C Y(N) INITIAL VALUES FOR Y |
---|
118 | C |
---|
119 | C XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) |
---|
120 | C |
---|
121 | C H INITIAL STEP SIZE GUESS; |
---|
122 | C FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, |
---|
123 | C H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD. |
---|
124 | C THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS |
---|
125 | C QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6). |
---|
126 | C |
---|
127 | C RelTol,AbsTol RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY |
---|
128 | C CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N. |
---|
129 | C |
---|
130 | C ITOL SWITCH FOR RelTol AND AbsTol: |
---|
131 | C ITOL=0: BOTH RelTol AND AbsTol ARE SCALARS. |
---|
132 | C THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF |
---|
133 | C Y(I) BELOW RelTol*ABS(Y(I))+AbsTol |
---|
134 | C ITOL=1: BOTH RelTol AND AbsTol ARE VECTORS. |
---|
135 | C THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW |
---|
136 | C RelTol(I)*ABS(Y(I))+AbsTol(I). |
---|
137 | C |
---|
138 | C JAC NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES |
---|
139 | C THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y |
---|
140 | C (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY |
---|
141 | C A DUMMY SUBROUTINE IN THE CASE IJAC=0). |
---|
142 | C FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM |
---|
143 | C SUBROUTINE JAC(N,X,Y,DFY,LDFY,RPAR,IPAR) |
---|
144 | C KPP_REAL X,Y(N),DFY(LDFY,N) |
---|
145 | C DFY(1,1)= ... |
---|
146 | C LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS |
---|
147 | C FURNISHED BY THE CALLING PROGRAM. |
---|
148 | C IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO |
---|
149 | C BE FULL AND THE PARTIAL DERIVATIVES ARE |
---|
150 | C STORED IN DFY AS |
---|
151 | C DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J) |
---|
152 | C ELSE, THE JACOBIAN IS TAKEN AS BANDED AND |
---|
153 | C THE PARTIAL DERIVATIVES ARE STORED |
---|
154 | C DIAGONAL-WISE AS |
---|
155 | C DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J). |
---|
156 | C |
---|
157 | C IJAC SWITCH FOR THE COMPUTATION OF THE JACOBIAN: |
---|
158 | C IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE |
---|
159 | C DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED. |
---|
160 | C IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC. |
---|
161 | C |
---|
162 | C MLJAC SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN: |
---|
163 | C MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR |
---|
164 | C ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. |
---|
165 | C 0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN |
---|
166 | C MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW |
---|
167 | C THE MAIN DIAGONAL). |
---|
168 | C |
---|
169 | C MUJAC UPPER BANDWITH OF JACOBIAN MATRIX (>= NUMBER OF NON- |
---|
170 | C ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). |
---|
171 | C NEED NOT BE DEFINED IF MLJAC=N. |
---|
172 | C |
---|
173 | C ---- MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS ----- |
---|
174 | C ---- FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): - |
---|
175 | C |
---|
176 | C MAS NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS- |
---|
177 | C MATRIX M. |
---|
178 | C IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY |
---|
179 | C MATRIX AND NEEDS NOT TO BE DEFINED; |
---|
180 | C SUPPLY A DUMMY SUBROUTINE IN THIS CASE. |
---|
181 | C IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM |
---|
182 | C SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR) |
---|
183 | C KPP_REAL AM(LMAS,N) |
---|
184 | C AM(1,1)= .... |
---|
185 | C IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED |
---|
186 | C AS FULL MATRIX LIKE |
---|
187 | C AM(I,J) = M(I,J) |
---|
188 | C ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED |
---|
189 | C DIAGONAL-WISE AS |
---|
190 | C AM(I-J+MUMAS+1,J) = M(I,J). |
---|
191 | C |
---|
192 | C IMAS GIVES INFORMATION ON THE MASS-MATRIX: |
---|
193 | C IMAS=0: M IS SUPPOSED TO BE THE IDENTITY |
---|
194 | C MATRIX, MAS IS NEVER CALLED. |
---|
195 | C IMAS=1: MASS-MATRIX IS SUPPLIED. |
---|
196 | C |
---|
197 | C MLMAS SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX: |
---|
198 | C MLMAS=N: THE FULL MATRIX CASE. THE LINEAR |
---|
199 | C ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. |
---|
200 | C 0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE |
---|
201 | C MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW |
---|
202 | C THE MAIN DIAGONAL). |
---|
203 | C MLMAS IS SUPPOSED TO BE .LE. MLJAC. |
---|
204 | C |
---|
205 | C MUMAS UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON- |
---|
206 | C ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). |
---|
207 | C NEED NOT BE DEFINED IF MLMAS=N. |
---|
208 | C MUMAS IS SUPPOSED TO BE .LE. MUJAC. |
---|
209 | C |
---|
210 | C SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE |
---|
211 | C NUMERICAL SOLUTION DURING INTEGRATION. |
---|
212 | C IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. |
---|
213 | C SUPPLY A DUMMY SUBROUTINE IF IOUT=0. |
---|
214 | C IT MUST HAVE THE FORM |
---|
215 | C SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N, |
---|
216 | C RPAR,IPAR,IRTRN) |
---|
217 | C KPP_REAL X,Y(N),CONT(LRC) |
---|
218 | C .... |
---|
219 | C SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH |
---|
220 | C GRID-POINT "X" (THEREBY THE INITIAL VALUE IS |
---|
221 | C THE FIRST GRID-POINT). |
---|
222 | C "XOLD" IS THE PRECEEDING GRID-POINT. |
---|
223 | C "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN |
---|
224 | C IS SET <0, RADAU5 RETURNS TO THE CALLING PROGRAM. |
---|
225 | C |
---|
226 | C ----- CONTINUOUS OUTPUT: ----- |
---|
227 | C DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION |
---|
228 | C FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH |
---|
229 | C THE FUNCTION |
---|
230 | C >>> CONTR5(I,S,CONT,LRC) <<< |
---|
231 | C WHICH PROVIDES AN APPROXIMATION TO THE I-TH |
---|
232 | C COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE |
---|
233 | C S SHOULD LIE IN THE INTERVAL [XOLD,X]. |
---|
234 | C DO NOT CHANGE THE ENTRIES OF CONT(LRC), IF THE |
---|
235 | C DENSE OUTPUT FUNCTION IS USED. |
---|
236 | C |
---|
237 | C IOUT SWITCH FOR CALLING THE SUBROUTINE SOLOUT: |
---|
238 | C IOUT=0: SUBROUTINE IS NEVER CALLED |
---|
239 | C IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT. |
---|
240 | C |
---|
241 | C WORK ARRAY OF WORKING SPACE OF LENGTH "LWORK". |
---|
242 | C WORK(1), WORK(2),.., WORK(20) SERVE AS PARAMETERS |
---|
243 | C FOR THE CODE. FOR STANDARD USE OF THE CODE |
---|
244 | C WORK(1),..,WORK(20) MUST BE SET TO ZERO BEFORE |
---|
245 | C CALLING. SEE BELOW FOR A MORE SOPHISTICATED USE. |
---|
246 | C WORK(21),..,WORK(LWORK) SERVE AS WORKING SPACE |
---|
247 | C FOR ALL VECTORS AND MATRICES. |
---|
248 | C "LWORK" MUST BE AT LEAST |
---|
249 | C N*(LJAC+LMAS+3*LE+12)+20 |
---|
250 | C WHERE |
---|
251 | C LJAC=N IF MLJAC=N (FULL JACOBIAN) |
---|
252 | C LJAC=MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.) |
---|
253 | C AND |
---|
254 | C LMAS=0 IF IMAS=0 |
---|
255 | C LMAS=N IF IMAS=1 AND MLMAS=N (FULL) |
---|
256 | C LMAS=MLMAS+MUMAS+1 IF MLMAS<N (BANDED MASS-M.) |
---|
257 | C AND |
---|
258 | C LE=N IF MLJAC=N (FULL JACOBIAN) |
---|
259 | C LE=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.) |
---|
260 | C |
---|
261 | C IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE |
---|
262 | C MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM |
---|
263 | C STORAGE REQUIREMENT IS |
---|
264 | C LWORK = 4*N*N+12*N+20. |
---|
265 | C IF IWORK(9)=M1>0 THEN "LWORK" MUST BE AT LEAST |
---|
266 | C N*(LJAC+12)+(N-M1)*(LMAS+3*LE)+20 |
---|
267 | C WHERE IN THE DEFINITIONS OF LJAC, LMAS AND LE THE |
---|
268 | C NUMBER N CAN BE REPLACED BY N-M1. |
---|
269 | C |
---|
270 | C LWORK DECLARED LENGTH OF ARRAY "WORK". |
---|
271 | C |
---|
272 | C IWORK INTEGER WORKING SPACE OF LENGTH "LIWORK". |
---|
273 | C IWORK(1),IWORK(2),...,IWORK(20) SERVE AS PARAMETERS |
---|
274 | C FOR THE CODE. FOR STANDARD USE, SET IWORK(1),.., |
---|
275 | C IWORK(20) TO ZERO BEFORE CALLING. |
---|
276 | C IWORK(21),...,IWORK(LIWORK) SERVE AS WORKING AREA. |
---|
277 | C "LIWORK" MUST BE AT LEAST 3*N+20. |
---|
278 | C |
---|
279 | C LIWORK DECLARED LENGTH OF ARRAY "IWORK". |
---|
280 | C |
---|
281 | C RPAR, IPAR REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH |
---|
282 | C CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING |
---|
283 | C PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES. |
---|
284 | C |
---|
285 | C ---------------------------------------------------------------------- |
---|
286 | C |
---|
287 | C SOPHISTICATED SETTING OF PARAMETERS |
---|
288 | C ----------------------------------- |
---|
289 | C SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK |
---|
290 | C WELL. THEY MAY BE DEFINED BY SETTING WORK(1),... |
---|
291 | C AS WELL AS IWORK(1),... DIFFERENT FROM ZERO. |
---|
292 | C FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES: |
---|
293 | C |
---|
294 | C IWORK(1) IF IWORK(1).NE.0, THE CODE TRANSFORMS THE JACOBIAN |
---|
295 | C MATRIX TO HESSENBERG FORM. THIS IS PARTICULARLY |
---|
296 | C ADVANTAGEOUS FOR LARGE SYSTEMS WITH FULL JACOBIAN. |
---|
297 | C IT DOES NOT WORK FOR BANDED JACOBIAN (MLJAC<N) |
---|
298 | C AND NOT FOR IMPLICIT SYSTEMS (IMAS=1). |
---|
299 | C |
---|
300 | C IWORK(2) THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS. |
---|
301 | C THE DEFAULT VALUE (FOR IWORK(2)=0) IS 100000. |
---|
302 | C |
---|
303 | C IWORK(3) THE MAXIMUM NUMBER OF NEWTON ITERATIONS FOR THE |
---|
304 | C SOLUTION OF THE IMPLICIT SYSTEM IN EACH STEP. |
---|
305 | C THE DEFAULT VALUE (FOR IWORK(3)=0) IS 7. |
---|
306 | C |
---|
307 | C IWORK(4) IF IWORK(4).EQ.0 THE EXTRAPOLATED COLLOCATION SOLUTION |
---|
308 | C IS TAKEN AS STARTING VALUE FOR NEWTON'S METHOD. |
---|
309 | C IF IWORK(4).NE.0 ZERO STARTING VALUES ARE USED. |
---|
310 | C THE LATTER IS RECOMMENDED IF NEWTON'S METHOD HAS |
---|
311 | C DIFFICULTIES WITH CONVERGENCE (THIS IS THE CASE WHEN |
---|
312 | C NSTEP IS LARGER THAN NACCPT + NREJCT; SEE OUTPUT PARAM.). |
---|
313 | C DEFAULT IS IWORK(4)=0. |
---|
314 | C |
---|
315 | C THE FOLLOWING 3 PARAMETERS ARE IMPORTANT FOR |
---|
316 | C DIFFERENTIAL-ALGEBRAIC SYSTEMS OF INDEX > 1. |
---|
317 | C THE FUNCTION-SUBROUTINE SHOULD BE WRITTEN SUCH THAT |
---|
318 | C THE INDEX 1,2,3 VARIABLES APPEAR IN THIS ORDER. |
---|
319 | C IN ESTIMATING THE ERROR THE INDEX 2 VARIABLES ARE |
---|
320 | C MULTIPLIED BY H, THE INDEX 3 VARIABLES BY H**2. |
---|
321 | C |
---|
322 | C IWORK(5) DIMENSION OF THE INDEX 1 VARIABLES (MUST BE > 0). FOR |
---|
323 | C ODE'S THIS EQUALS THE DIMENSION OF THE SYSTEM. |
---|
324 | C DEFAULT IWORK(5)=N. |
---|
325 | C |
---|
326 | C IWORK(6) DIMENSION OF THE INDEX 2 VARIABLES. DEFAULT IWORK(6)=0. |
---|
327 | C |
---|
328 | C IWORK(7) DIMENSION OF THE INDEX 3 VARIABLES. DEFAULT IWORK(7)=0. |
---|
329 | C |
---|
330 | C IWORK(8) SWITCH FOR STEP SIZE STRATEGY |
---|
331 | C IF IWORK(8).EQ.1 MOD. PREDICTIVE CONTROLLER (GUSTAFSSON) |
---|
332 | C IF IWORK(8).EQ.2 CLASSICAL STEP SIZE CONTROL |
---|
333 | C THE DEFAULT VALUE (FOR IWORK(8)=0) IS IWORK(8)=1. |
---|
334 | C THE CHOICE IWORK(8).EQ.1 SEEMS TO PRODUCE SAFER RESULTS; |
---|
335 | C FOR SIMPLE PROBLEMS, THE CHOICE IWORK(8).EQ.2 PRODUCES |
---|
336 | C OFTEN SLIGHTLY FASTER RUNS |
---|
337 | C |
---|
338 | C IF THE DIFFERENTIAL SYSTEM HAS THE SPECIAL STRUCTURE THAT |
---|
339 | C Y(I)' = Y(I+M2) FOR I=1,...,M1, |
---|
340 | C WITH M1 A MULTIPLE OF M2, A SUBSTANTIAL GAIN IN COMPUTERTIME |
---|
341 | C CAN BE ACHIEVED BY SETTING THE FOLLOWING TWO PARAMETERS. E.G., |
---|
342 | C FOR SECOND ORDER SYSTEMS P'=V, V'=G(P,V), WHERE P AND V ARE |
---|
343 | C VECTORS OF DIMENSION N/2, ONE HAS TO PUT M1=M2=N/2. |
---|
344 | C FOR M1>0 SOME OF THE INPUT PARAMETERS HAVE DIFFERENT MEANINGS: |
---|
345 | C - JAC: ONLY THE ELEMENTS OF THE NON-TRIVIAL PART OF THE |
---|
346 | C JACOBIAN HAVE TO BE STORED |
---|
347 | C IF (MLJAC.EQ.N-M1) THE JACOBIAN IS SUPPOSED TO BE FULL |
---|
348 | C DFY(I,J) = PARTIAL F(I+M1) / PARTIAL Y(J) |
---|
349 | C FOR I=1,N-M1 AND J=1,N. |
---|
350 | C ELSE, THE JACOBIAN IS BANDED ( M1 = M2 * MM ) |
---|
351 | C DFY(I-J+MUJAC+1,J+K*M2) = PARTIAL F(I+M1) / PARTIAL Y(J+K*M2) |
---|
352 | C FOR I=1,MLJAC+MUJAC+1 AND J=1,M2 AND K=0,MM. |
---|
353 | C - MLJAC: MLJAC=N-M1: IF THE NON-TRIVIAL PART OF THE JACOBIAN IS FULL |
---|
354 | C 0<=MLJAC<N-M1: IF THE (MM+1) SUBMATRICES (FOR K=0,MM) |
---|
355 | C PARTIAL F(I+M1) / PARTIAL Y(J+K*M2), I,J=1,M2 |
---|
356 | C ARE BANDED, MLJAC IS THE MAXIMAL LOWER BANDWIDTH |
---|
357 | C OF THESE MM+1 SUBMATRICES |
---|
358 | C - MUJAC: MAXIMAL UPPER BANDWIDTH OF THESE MM+1 SUBMATRICES |
---|
359 | C NEED NOT BE DEFINED IF MLJAC=N-M1 |
---|
360 | C - MAS: IF IMAS=0 THIS MATRIX IS ASSUMED TO BE THE IDENTITY AND |
---|
361 | C NEED NOT BE DEFINED. SUPPLY A DUMMY SUBROUTINE IN THIS CASE. |
---|
362 | C IT IS ASSUMED THAT ONLY THE ELEMENTS OF RIGHT LOWER BLOCK OF |
---|
363 | C DIMENSION N-M1 DIFFER FROM THAT OF THE IDENTITY MATRIX. |
---|
364 | C IF (MLMAS.EQ.N-M1) THIS SUBMATRIX IS SUPPOSED TO BE FULL |
---|
365 | C AM(I,J) = M(I+M1,J+M1) FOR I=1,N-M1 AND J=1,N-M1. |
---|
366 | C ELSE, THE MASS MATRIX IS BANDED |
---|
367 | C AM(I-J+MUMAS+1,J) = M(I+M1,J+M1) |
---|
368 | C - MLMAS: MLMAS=N-M1: IF THE NON-TRIVIAL PART OF M IS FULL |
---|
369 | C 0<=MLMAS<N-M1: LOWER BANDWIDTH OF THE MASS MATRIX |
---|
370 | C - MUMAS: UPPER BANDWIDTH OF THE MASS MATRIX |
---|
371 | C NEED NOT BE DEFINED IF MLMAS=N-M1 |
---|
372 | C |
---|
373 | C IWORK(9) THE VALUE OF M1. DEFAULT M1=0. |
---|
374 | C |
---|
375 | C IWORK(10) THE VALUE OF M2. DEFAULT M2=M1. |
---|
376 | C |
---|
377 | C ---------- |
---|
378 | C |
---|
379 | C WORK(1) UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16. |
---|
380 | C |
---|
381 | C WORK(2) THE SAFETY FACTOR IN STEP SIZE PREDICTION, |
---|
382 | C DEFAULT 0.9D0. |
---|
383 | C |
---|
384 | C WORK(3) DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED; |
---|
385 | C INCREASE WORK(3), TO 0.1 SAY, WHEN JACOBIAN EVALUATIONS |
---|
386 | C ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER |
---|
387 | C (0.001D0, SAY). NEGATIV WORK(3) FORCES THE CODE TO |
---|
388 | C COMPUTE THE JACOBIAN AFTER EVERY ACCEPTED STEP. |
---|
389 | C DEFAULT 0.001D0. |
---|
390 | C |
---|
391 | C WORK(4) STOPPING CRITERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1. |
---|
392 | C SMALLER VALUES OF WORK(4) MAKE THE CODE SLOWER, BUT SAFER. |
---|
393 | C DEFAULT 0.03D0. |
---|
394 | C |
---|
395 | C WORK(5) AND WORK(6) : IF WORK(5) < HNEW/HOLD < WORK(6), THEN THE |
---|
396 | C STEP SIZE IS NOT CHANGED. THIS SAVES, TOGETHER WITH A |
---|
397 | C LARGE WORK(3), LU-DECOMPOSITIONS AND COMPUTING TIME FOR |
---|
398 | C LARGE SYSTEMS. FOR SMALL SYSTEMS ONE MAY HAVE |
---|
399 | C WORK(5)=1.D0, WORK(6)=1.2D0, FOR LARGE FULL SYSTEMS |
---|
400 | C WORK(5)=0.99D0, WORK(6)=2.D0 MIGHT BE GOOD. |
---|
401 | C DEFAULTS WORK(5)=1.D0, WORK(6)=1.2D0 . |
---|
402 | C |
---|
403 | C WORK(7) MAXIMAL STEP SIZE, DEFAULT XEND-X. |
---|
404 | C |
---|
405 | C WORK(8), WORK(9) PARAMETERS FOR STEP SIZE SELECTION |
---|
406 | C THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION |
---|
407 | C WORK(8) <= HNEW/HOLD <= WORK(9) |
---|
408 | C DEFAULT VALUES: WORK(8)=0.2D0, WORK(9)=8.D0 |
---|
409 | C |
---|
410 | C----------------------------------------------------------------------- |
---|
411 | C |
---|
412 | C OUTPUT PARAMETERS |
---|
413 | C ----------------- |
---|
414 | C X X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED |
---|
415 | C (AFTER SUCCESSFUL RETURN X=XEND). |
---|
416 | C |
---|
417 | C Y(N) NUMERICAL SOLUTION AT X |
---|
418 | C |
---|
419 | C H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP |
---|
420 | C |
---|
421 | C IDID REPORTS ON SUCCESSFULNESS UPON RETURN: |
---|
422 | C IDID= 1 COMPUTATION SUCCESSFUL, |
---|
423 | C IDID= 2 COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT) |
---|
424 | C IDID=-1 INPUT IS NOT CONSISTENT, |
---|
425 | C IDID=-2 LARGER NMAX IS NEEDED, |
---|
426 | C IDID=-3 STEP SIZE BECOMES TOO SMALL, |
---|
427 | C IDID=-4 MATRIX IS REPEATEDLY SINGULAR. |
---|
428 | C |
---|
429 | C IWORK(14) NFCN NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL |
---|
430 | C EVALUATION OF THE JACOBIAN ARE NOT COUNTED) |
---|
431 | C IWORK(15) NJAC NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY |
---|
432 | C OR NUMERICALLY) |
---|
433 | C IWORK(16) NSTEP NUMBER OF COMPUTED STEPS |
---|
434 | C IWORK(17) NACCPT NUMBER OF ACCEPTED STEPS |
---|
435 | C IWORK(18) NREJCT NUMBER OF REJECTED STEPS (DUE TO ERROR TEST), |
---|
436 | C (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED) |
---|
437 | C IWORK(19) NDEC NUMBER OF LU-DECOMPOSITIONS OF BOTH MATRICES |
---|
438 | C IWORK(20) NSOL NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS, OF BOTH |
---|
439 | C SYSTEMS; THE NSTEP FORWARD-BACKWARD SUBSTITUTIONS, |
---|
440 | C NEEDED FOR STEP SIZE SELECTION, ARE NOT COUNTED |
---|
441 | C----------------------------------------------------------------------- |
---|
442 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
443 | C DECLARATIONS |
---|
444 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
445 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
446 | DIMENSION Y(N),AbsTol(*),RelTol(*),WORK(LWORK),IWORK(LIWORK) |
---|
447 | DIMENSION RPAR(*),IPAR(*) |
---|
448 | LOGICAL IMPLCT,JBAND,ARRET,STARTN,PRED |
---|
449 | EXTERNAL FCN,JAC,MAS,SOLOUT |
---|
450 | C *** *** *** *** *** *** *** |
---|
451 | C SETTING THE PARAMETERS |
---|
452 | C *** *** *** *** *** *** *** |
---|
453 | NFCN=0 |
---|
454 | NJAC=0 |
---|
455 | NSTEP=0 |
---|
456 | NACCPT=0 |
---|
457 | NREJCT=0 |
---|
458 | NDEC=0 |
---|
459 | NSOL=0 |
---|
460 | ARRET=.FALSE. |
---|
461 | C -------- NMAX , THE MAXIMAL NUMBER OF STEPS ----- |
---|
462 | IF (IWORK(2).EQ.0) THEN |
---|
463 | NMAX=100000 |
---|
464 | ELSE |
---|
465 | NMAX=IWORK(2) |
---|
466 | IF (NMAX.LE.0) THEN |
---|
467 | WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2) |
---|
468 | ARRET=.TRUE. |
---|
469 | END IF |
---|
470 | END IF |
---|
471 | C -------- NIT MAXIMAL NUMBER OF NEWTON ITERATIONS |
---|
472 | IF (IWORK(3).EQ.0) THEN |
---|
473 | NIT=7 |
---|
474 | ELSE |
---|
475 | NIT=IWORK(3) |
---|
476 | IF (NIT.LE.0) THEN |
---|
477 | WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3) |
---|
478 | ARRET=.TRUE. |
---|
479 | END IF |
---|
480 | END IF |
---|
481 | C -------- STARTN SWITCH FOR STARTING VALUES OF NEWTON ITERATIONS |
---|
482 | IF(IWORK(4).EQ.0)THEN |
---|
483 | STARTN=.FALSE. |
---|
484 | ELSE |
---|
485 | STARTN=.TRUE. |
---|
486 | END IF |
---|
487 | C -------- PARAMETER FOR DIFFERENTIAL-ALGEBRAIC COMPONENTS |
---|
488 | NIND1=IWORK(5) |
---|
489 | NIND2=IWORK(6) |
---|
490 | NIND3=IWORK(7) |
---|
491 | IF (NIND1.EQ.0) NIND1=N |
---|
492 | IF (NIND1+NIND2+NIND3.NE.N) THEN |
---|
493 | WRITE(6,*)' CURIOUS INPUT FOR IWORK(5,6,7)=',NIND1,NIND2,NIND3 |
---|
494 | ARRET=.TRUE. |
---|
495 | END IF |
---|
496 | C -------- PRED STEP SIZE CONTROL |
---|
497 | IF(IWORK(8).LE.1)THEN |
---|
498 | PRED=.TRUE. |
---|
499 | ELSE |
---|
500 | PRED=.FALSE. |
---|
501 | END IF |
---|
502 | C -------- PARAMETER FOR SECOND ORDER EQUATIONS |
---|
503 | M1=IWORK(9) |
---|
504 | M2=IWORK(10) |
---|
505 | NM1=N-M1 |
---|
506 | IF (M1.EQ.0) M2=N |
---|
507 | IF (M2.EQ.0) M2=M1 |
---|
508 | IF (M1.LT.0.OR.M2.LT.0.OR.M1+M2.GT.N) THEN |
---|
509 | WRITE(6,*)' CURIOUS INPUT FOR IWORK(9,10)=',M1,M2 |
---|
510 | ARRET=.TRUE. |
---|
511 | END IF |
---|
512 | C -------- UROUND SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0 |
---|
513 | IF (WORK(1).EQ.0.0D0) THEN |
---|
514 | UROUND=1.0D-16 |
---|
515 | ELSE |
---|
516 | UROUND=WORK(1) |
---|
517 | IF (UROUND.LE.1.0D-19.OR.UROUND.GE.1.0D0) THEN |
---|
518 | WRITE(6,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1) |
---|
519 | ARRET=.TRUE. |
---|
520 | END IF |
---|
521 | END IF |
---|
522 | C --------- SAFE SAFETY FACTOR IN STEP SIZE PREDICTION |
---|
523 | IF (WORK(2).EQ.0.0D0) THEN |
---|
524 | SAFE=0.9D0 |
---|
525 | ELSE |
---|
526 | SAFE=WORK(2) |
---|
527 | IF (SAFE.LE.0.001D0.OR.SAFE.GE.1.0D0) THEN |
---|
528 | WRITE(6,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2) |
---|
529 | ARRET=.TRUE. |
---|
530 | END IF |
---|
531 | END IF |
---|
532 | C ------ THET DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED; |
---|
533 | IF (WORK(3).EQ.0.D0) THEN |
---|
534 | THET=0.001D0 |
---|
535 | ELSE |
---|
536 | THET=WORK(3) |
---|
537 | IF (THET.LE.0.0D0.OR.THET.GE.1.0D0) THEN |
---|
538 | WRITE(6,*)' CURIOUS INPUT FOR WORK(3)=',WORK(3) |
---|
539 | ARRET=.TRUE. |
---|
540 | END IF |
---|
541 | END IF |
---|
542 | C --- FNEWT STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1. |
---|
543 | IF (WORK(4).EQ.0.D0) THEN |
---|
544 | FNEWT=0.03D0 |
---|
545 | ELSE |
---|
546 | FNEWT=WORK(4) |
---|
547 | IF (FNEWT.LE.UROUND) THEN |
---|
548 | WRITE(6,*)' CURIOUS INPUT FOR WORK(4)=',WORK(4) |
---|
549 | ARRET=.TRUE. |
---|
550 | END IF |
---|
551 | END IF |
---|
552 | C --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST. |
---|
553 | IF (WORK(5).EQ.0.D0) THEN |
---|
554 | QUOT1=1.D0 |
---|
555 | ELSE |
---|
556 | QUOT1=WORK(5) |
---|
557 | END IF |
---|
558 | IF (WORK(6).EQ.0.D0) THEN |
---|
559 | QUOT2=1.2D0 |
---|
560 | ELSE |
---|
561 | QUOT2=WORK(6) |
---|
562 | END IF |
---|
563 | IF (QUOT1.GT.1.0D0.OR.QUOT2.LT.1.0D0) THEN |
---|
564 | WRITE(6,*)' CURIOUS INPUT FOR WORK(5,6)=',QUOT1,QUOT2 |
---|
565 | ARRET=.TRUE. |
---|
566 | END IF |
---|
567 | C -------- MAXIMAL STEP SIZE |
---|
568 | IF (WORK(7).EQ.0.D0) THEN |
---|
569 | HMAX=XEND-X |
---|
570 | ELSE |
---|
571 | HMAX=WORK(7) |
---|
572 | END IF |
---|
573 | C ------- FACL,FACR PARAMETERS FOR STEP SIZE SELECTION |
---|
574 | IF(WORK(8).EQ.0.D0)THEN |
---|
575 | FACL=5.D0 |
---|
576 | ELSE |
---|
577 | FACL=1.D0/WORK(8) |
---|
578 | END IF |
---|
579 | IF(WORK(9).EQ.0.D0)THEN |
---|
580 | FACR=1.D0/8.0D0 |
---|
581 | ELSE |
---|
582 | FACR=1.D0/WORK(9) |
---|
583 | END IF |
---|
584 | IF (FACL.LT.1.0D0.OR.FACR.GT.1.0D0) THEN |
---|
585 | WRITE(6,*)' CURIOUS INPUT WORK(8,9)=',WORK(8),WORK(9) |
---|
586 | ARRET=.TRUE. |
---|
587 | END IF |
---|
588 | C --------- CHECK IF TOLERANCES ARE O.K. |
---|
589 | IF (ITOL.EQ.0) THEN |
---|
590 | IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN |
---|
591 | WRITE (6,*) ' TOLERANCES ARE TOO SMALL' |
---|
592 | ARRET=.TRUE. |
---|
593 | END IF |
---|
594 | ELSE |
---|
595 | DO I=1,N |
---|
596 | IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN |
---|
597 | WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL' |
---|
598 | ARRET=.TRUE. |
---|
599 | END IF |
---|
600 | END DO |
---|
601 | END IF |
---|
602 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
603 | C COMPUTATION OF ARRAY ENTRIES |
---|
604 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
605 | C ---- IMPLICIT, BANDED OR NOT ? |
---|
606 | IMPLCT=IMAS.NE.0 |
---|
607 | JBAND=MLJAC.LT.NM1 |
---|
608 | C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS --- |
---|
609 | C -- JACOBIAN AND MATRICES E1, E2 |
---|
610 | IF (JBAND) THEN |
---|
611 | LDJAC=MLJAC+MUJAC+1 |
---|
612 | LDE1=MLJAC+LDJAC |
---|
613 | ELSE |
---|
614 | MLJAC=NM1 |
---|
615 | MUJAC=NM1 |
---|
616 | LDJAC=NM1 |
---|
617 | LDE1=NM1 |
---|
618 | END IF |
---|
619 | C -- MASS MATRIX |
---|
620 | IF (IMPLCT) THEN |
---|
621 | IF (MLMAS.NE.NM1) THEN |
---|
622 | LDMAS=MLMAS+MUMAS+1 |
---|
623 | IF (JBAND) THEN |
---|
624 | IJOB=4 |
---|
625 | ELSE |
---|
626 | IJOB=3 |
---|
627 | END IF |
---|
628 | ELSE |
---|
629 | LDMAS=NM1 |
---|
630 | IJOB=5 |
---|
631 | END IF |
---|
632 | C ------ BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF "JAC" |
---|
633 | IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN |
---|
634 | WRITE (6,*) 'BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF |
---|
635 | & "JAC"' |
---|
636 | ARRET=.TRUE. |
---|
637 | END IF |
---|
638 | ELSE |
---|
639 | LDMAS=0 |
---|
640 | IF (JBAND) THEN |
---|
641 | IJOB=2 |
---|
642 | ELSE |
---|
643 | IJOB=1 |
---|
644 | IF (N.GT.2.AND.IWORK(1).NE.0) IJOB=7 |
---|
645 | END IF |
---|
646 | END IF |
---|
647 | LDMAS2=MAX(1,LDMAS) |
---|
648 | C ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN |
---|
649 | IF ((IMPLCT.OR.JBAND).AND.IJOB.EQ.7) THEN |
---|
650 | WRITE(6,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH |
---|
651 | &FULL JACOBIAN' |
---|
652 | ARRET=.TRUE. |
---|
653 | END IF |
---|
654 | C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- |
---|
655 | IEZ1=21 |
---|
656 | IEZ2=IEZ1+N |
---|
657 | IEZ3=IEZ2+N |
---|
658 | IEY0=IEZ3+N |
---|
659 | IESCAL=IEY0+N |
---|
660 | IEF1=IESCAL+N |
---|
661 | IEF2=IEF1+N |
---|
662 | IEF3=IEF2+N |
---|
663 | IECON=IEF3+N |
---|
664 | IEJAC=IECON+4*N |
---|
665 | IEMAS=IEJAC+N*LDJAC |
---|
666 | IEE1=IEMAS+NM1*LDMAS |
---|
667 | IEE2R=IEE1+NM1*LDE1 |
---|
668 | IEE2I=IEE2R+NM1*LDE1 |
---|
669 | C ------ TOTAL STORAGE REQUIREMENT ----------- |
---|
670 | ISTORE=IEE2I+NM1*LDE1-1 |
---|
671 | IF(ISTORE.GT.LWORK)THEN |
---|
672 | WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE |
---|
673 | ARRET=.TRUE. |
---|
674 | END IF |
---|
675 | C ------- ENTRY POINTS FOR INTEGER WORKSPACE ----- |
---|
676 | IEIP1=21 |
---|
677 | IEIP2=IEIP1+NM1 |
---|
678 | IEIPH=IEIP2+NM1 |
---|
679 | C --------- TOTAL REQUIREMENT --------------- |
---|
680 | ISTORE=IEIPH+NM1-1 |
---|
681 | IF (ISTORE.GT.LIWORK) THEN |
---|
682 | WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE |
---|
683 | ARRET=.TRUE. |
---|
684 | END IF |
---|
685 | C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 |
---|
686 | IF (ARRET) THEN |
---|
687 | IDID=-1 |
---|
688 | RETURN |
---|
689 | END IF |
---|
690 | C -------- CALL TO CORE INTEGRATOR ------------ |
---|
691 | CALL RADCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL, |
---|
692 | & JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID, |
---|
693 | & NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,IJOB,STARTN, |
---|
694 | & NIND1,NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1, |
---|
695 | & IMPLCT,JBAND,LDJAC,LDE1,LDMAS2,WORK(IEZ1),WORK(IEZ2), |
---|
696 | & WORK(IEZ3),WORK(IEY0),WORK(IESCAL),WORK(IEF1),WORK(IEF2), |
---|
697 | & WORK(IEF3),WORK(IEJAC),WORK(IEE1),WORK(IEE2R),WORK(IEE2I), |
---|
698 | & WORK(IEMAS),IWORK(IEIP1),IWORK(IEIP2),IWORK(IEIPH), |
---|
699 | & WORK(IECON),NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,RPAR,IPAR) |
---|
700 | IWORK(14)=NFCN |
---|
701 | IWORK(15)=NJAC |
---|
702 | IWORK(16)=NSTEP |
---|
703 | IWORK(17)=NACCPT |
---|
704 | IWORK(18)=NREJCT |
---|
705 | IWORK(19)=NDEC |
---|
706 | IWORK(20)=NSOL |
---|
707 | C ----------- RETURN ----------- |
---|
708 | RETURN |
---|
709 | END |
---|
710 | C |
---|
711 | C END OF SUBROUTINE RADAU5 |
---|
712 | C |
---|
713 | C *********************************************************** |
---|
714 | C |
---|
715 | SUBROUTINE RADCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL, |
---|
716 | & JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID, |
---|
717 | & NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,IJOB,STARTN, |
---|
718 | & NIND1,NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1, |
---|
719 | & IMPLCT,BANDED,LDJAC,LDE1,LDMAS,Z1,Z2,Z3, |
---|
720 | & Y0,SCAL,F1,F2,F3,FJAC,E1,E2R,E2I,FMAS,IP1,IP2,IPHES, |
---|
721 | & CONT,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,RPAR,IPAR) |
---|
722 | C ---------------------------------------------------------- |
---|
723 | C CORE INTEGRATOR FOR RADAU5 |
---|
724 | C PARAMETERS SAME AS IN RADAU5 WITH WORKSPACE ADDED |
---|
725 | C ---------------------------------------------------------- |
---|
726 | C DECLARATIONS |
---|
727 | C ---------------------------------------------------------- |
---|
728 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
729 | DIMENSION Y(N),Z1(N),Z2(N),Z3(N),Y0(N),SCAL(N),F1(N),F2(N),F3(N) |
---|
730 | DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),CONT(4*N) |
---|
731 | DIMENSION E1(LDE1,NM1),E2R(LDE1,NM1),E2I(LDE1,NM1) |
---|
732 | DIMENSION AbsTol(*),RelTol(*),RPAR(*),IPAR(*) |
---|
733 | INTEGER IP1(NM1),IP2(NM1),IPHES(NM1) |
---|
734 | COMMON /CONRA5/NN,NN2,NN3,NN4,XSOL,HSOL,C2M1,C1M1 |
---|
735 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
736 | LOGICAL REJECT,FIRST,IMPLCT,BANDED,CALJAC,STARTN,CALHES |
---|
737 | LOGICAL INDEX1,INDEX2,INDEX3,LAST,PRED |
---|
738 | EXTERNAL FCN, JAC, MAS, SOLOUT |
---|
739 | C *** *** *** *** *** *** *** |
---|
740 | C INITIALISATIONS |
---|
741 | C *** *** *** *** *** *** *** |
---|
742 | C --------- DUPLIFY N FOR COMMON BLOCK CONT ----- |
---|
743 | NN=N |
---|
744 | NN2=2*N |
---|
745 | NN3=3*N |
---|
746 | LRC=4*N |
---|
747 | C -------- CHECK THE INDEX OF THE PROBLEM ----- |
---|
748 | INDEX1=NIND1.NE.0 |
---|
749 | INDEX2=NIND2.NE.0 |
---|
750 | INDEX3=NIND3.NE.0 |
---|
751 | C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ---------- |
---|
752 | IF (IMPLCT) CALL MAS(NM1,FMAS,LDMAS,RPAR,IPAR) |
---|
753 | C ---------- CONSTANTS --------- |
---|
754 | SQ6=DSQRT(6.D0) |
---|
755 | C1=(4.D0-SQ6)/10.D0 |
---|
756 | C2=(4.D0+SQ6)/10.D0 |
---|
757 | C1M1=C1-1.D0 |
---|
758 | C2M1=C2-1.D0 |
---|
759 | C1MC2=C1-C2 |
---|
760 | DD1=-(13.D0+7.D0*SQ6)/3.D0 |
---|
761 | DD2=(-13.D0+7.D0*SQ6)/3.D0 |
---|
762 | DD3=-1.D0/3.D0 |
---|
763 | U1=(6.D0+81.D0**(1.D0/3.D0)-9.D0**(1.D0/3.D0))/30.D0 |
---|
764 | ALPH=(12.D0-81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))/60.D0 |
---|
765 | BETA=(81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))*DSQRT(3.D0)/60.D0 |
---|
766 | CNO=ALPH**2+BETA**2 |
---|
767 | U1=1.0D0/U1 |
---|
768 | ALPH=ALPH/CNO |
---|
769 | BETA=BETA/CNO |
---|
770 | T11=9.1232394870892942792D-02 |
---|
771 | T12=-0.14125529502095420843D0 |
---|
772 | T13=-3.0029194105147424492D-02 |
---|
773 | T21=0.24171793270710701896D0 |
---|
774 | T22=0.20412935229379993199D0 |
---|
775 | T23=0.38294211275726193779D0 |
---|
776 | T31=0.96604818261509293619D0 |
---|
777 | TI11=4.3255798900631553510D0 |
---|
778 | TI12=0.33919925181580986954D0 |
---|
779 | TI13=0.54177053993587487119D0 |
---|
780 | TI21=-4.1787185915519047273D0 |
---|
781 | TI22=-0.32768282076106238708D0 |
---|
782 | TI23=0.47662355450055045196D0 |
---|
783 | TI31=-0.50287263494578687595D0 |
---|
784 | TI32=2.5719269498556054292D0 |
---|
785 | TI33=-0.59603920482822492497D0 |
---|
786 | IF (M1.GT.0) IJOB=IJOB+10 |
---|
787 | POSNEG=SIGN(1.D0,XEND-X) |
---|
788 | HMAXN=MIN(ABS(HMAX),ABS(XEND-X)) |
---|
789 | IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6 |
---|
790 | H=MIN(ABS(H),HMAXN) |
---|
791 | H=SIGN(H,POSNEG) |
---|
792 | HOLD=H |
---|
793 | REJECT=.FALSE. |
---|
794 | FIRST=.TRUE. |
---|
795 | LAST=.FALSE. |
---|
796 | IF ((X+H*1.0001D0-XEND)*POSNEG.GE.0.D0) THEN |
---|
797 | H=XEND-X |
---|
798 | LAST=.TRUE. |
---|
799 | END IF |
---|
800 | FACCON=1.D0 |
---|
801 | CFAC=SAFE*(1+2*NIT) |
---|
802 | NSING=0 |
---|
803 | XOLD=X |
---|
804 | IF (IOUT.NE.0) THEN |
---|
805 | IRTRN=1 |
---|
806 | NRSOL=1 |
---|
807 | XOSOL=XOLD |
---|
808 | XSOL=X |
---|
809 | DO I=1,N |
---|
810 | CONT(I)=Y(I) |
---|
811 | END DO |
---|
812 | NSOLU=N |
---|
813 | HSOL=HOLD |
---|
814 | CALL SOLOUT(NRSOL,XOSOL,XSOL,Y,CONT,LRC,NSOLU, |
---|
815 | & RPAR,IPAR,IRTRN) |
---|
816 | IF (IRTRN.LT.0) GOTO 179 |
---|
817 | END IF |
---|
818 | MLE=MLJAC |
---|
819 | MUE=MUJAC |
---|
820 | MBJAC=MLJAC+MUJAC+1 |
---|
821 | MBB=MLMAS+MUMAS+1 |
---|
822 | MDIAG=MLE+MUE+1 |
---|
823 | MDIFF=MLE+MUE-MUMAS |
---|
824 | MBDIAG=MUMAS+1 |
---|
825 | N2=2*N |
---|
826 | N3=3*N |
---|
827 | IF (ITOL.EQ.0) THEN |
---|
828 | DO I=1,N |
---|
829 | SCAL(I)=AbsTol(1)+RelTol(1)*ABS(Y(I)) |
---|
830 | END DO |
---|
831 | ELSE |
---|
832 | DO I=1,N |
---|
833 | SCAL(I)=AbsTol(I)+RelTol(I)*ABS(Y(I)) |
---|
834 | END DO |
---|
835 | END IF |
---|
836 | HHFAC=H |
---|
837 | CALL FCN(N,X,Y,Y0,RPAR,IPAR) |
---|
838 | NFCN=NFCN+1 |
---|
839 | C --- BASIC INTEGRATION STEP |
---|
840 | 10 CONTINUE |
---|
841 | C *** *** *** *** *** *** *** |
---|
842 | C COMPUTATION OF THE JACOBIAN |
---|
843 | C *** *** *** *** *** *** *** |
---|
844 | NJAC=NJAC+1 |
---|
845 | IF (IJAC.EQ.0) THEN |
---|
846 | C --- COMPUTE JACOBIAN MATRIX NUMERICALLY |
---|
847 | IF (BANDED) THEN |
---|
848 | C --- JACOBIAN IS BANDED |
---|
849 | MUJACP=MUJAC+1 |
---|
850 | MD=MIN(MBJAC,M2) |
---|
851 | DO MM=1,M1/M2+1 |
---|
852 | DO K=1,MD |
---|
853 | J=K+(MM-1)*M2 |
---|
854 | 12 F1(J)=Y(J) |
---|
855 | F2(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J)))) |
---|
856 | Y(J)=Y(J)+F2(J) |
---|
857 | J=J+MD |
---|
858 | IF (J.LE.MM*M2) GOTO 12 |
---|
859 | CALL FCN(N,X,Y,CONT,RPAR,IPAR) |
---|
860 | J=K+(MM-1)*M2 |
---|
861 | J1=K |
---|
862 | LBEG=MAX(1,J1-MUJAC)+M1 |
---|
863 | 14 LEND=MIN(M2,J1+MLJAC)+M1 |
---|
864 | Y(J)=F1(J) |
---|
865 | MUJACJ=MUJACP-J1-M1 |
---|
866 | DO L=LBEG,LEND |
---|
867 | FJAC(L+MUJACJ,J)=(CONT(L)-Y0(L))/F2(J) |
---|
868 | END DO |
---|
869 | J=J+MD |
---|
870 | J1=J1+MD |
---|
871 | LBEG=LEND+1 |
---|
872 | IF (J.LE.MM*M2) GOTO 14 |
---|
873 | END DO |
---|
874 | END DO |
---|
875 | ELSE |
---|
876 | C --- JACOBIAN IS FULL |
---|
877 | DO I=1,N |
---|
878 | YSAFE=Y(I) |
---|
879 | DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE))) |
---|
880 | Y(I)=YSAFE+DELT |
---|
881 | CALL FCN(N,X,Y,CONT,RPAR,IPAR) |
---|
882 | DO J=M1+1,N |
---|
883 | FJAC(J-M1,I)=(CONT(J)-Y0(J))/DELT |
---|
884 | END DO |
---|
885 | Y(I)=YSAFE |
---|
886 | END DO |
---|
887 | END IF |
---|
888 | ELSE |
---|
889 | C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY |
---|
890 | CALL JAC_CHEM(N,X,Y,FJAC,LDJAC,RPAR,IPAR) |
---|
891 | END IF |
---|
892 | CALJAC=.TRUE. |
---|
893 | CALHES=.TRUE. |
---|
894 | 20 CONTINUE |
---|
895 | C --- COMPUTE THE MATRICES E1 AND E2 AND THEIR DECOMPOSITIONS |
---|
896 | FAC1=U1/H |
---|
897 | ALPHN=ALPH/H |
---|
898 | BETAN=BETA/H |
---|
899 | CALL DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
900 | & M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES) |
---|
901 | IF (IER.NE.0) GOTO 78 |
---|
902 | CALL DECOMC(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
903 | & M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,IP2,IER,IJOB) |
---|
904 | IF (IER.NE.0) GOTO 78 |
---|
905 | NDEC=NDEC+1 |
---|
906 | 30 CONTINUE |
---|
907 | NSTEP=NSTEP+1 |
---|
908 | IF (NSTEP.GT.NMAX) GOTO 178 |
---|
909 | IF (0.1D0*ABS(H).LE.ABS(X)*UROUND) GOTO 177 |
---|
910 | IF (INDEX2) THEN |
---|
911 | DO I=NIND1+1,NIND1+NIND2 |
---|
912 | SCAL(I)=SCAL(I)/HHFAC |
---|
913 | END DO |
---|
914 | END IF |
---|
915 | IF (INDEX3) THEN |
---|
916 | DO I=NIND1+NIND2+1,NIND1+NIND2+NIND3 |
---|
917 | SCAL(I)=SCAL(I)/(HHFAC*HHFAC) |
---|
918 | END DO |
---|
919 | END IF |
---|
920 | XPH=X+H |
---|
921 | C *** *** *** *** *** *** *** |
---|
922 | C STARTING VALUES FOR NEWTON ITERATION |
---|
923 | C *** *** *** *** *** *** *** |
---|
924 | IF (FIRST.OR.STARTN) THEN |
---|
925 | DO I=1,N |
---|
926 | Z1(I)=0.D0 |
---|
927 | Z2(I)=0.D0 |
---|
928 | Z3(I)=0.D0 |
---|
929 | F1(I)=0.D0 |
---|
930 | F2(I)=0.D0 |
---|
931 | F3(I)=0.D0 |
---|
932 | END DO |
---|
933 | ELSE |
---|
934 | C3Q=H/HOLD |
---|
935 | C1Q=C1*C3Q |
---|
936 | C2Q=C2*C3Q |
---|
937 | DO I=1,N |
---|
938 | AK1=CONT(I+N) |
---|
939 | AK2=CONT(I+N2) |
---|
940 | AK3=CONT(I+N3) |
---|
941 | Z1I=C1Q*(AK1+(C1Q-C2M1)*(AK2+(C1Q-C1M1)*AK3)) |
---|
942 | Z2I=C2Q*(AK1+(C2Q-C2M1)*(AK2+(C2Q-C1M1)*AK3)) |
---|
943 | Z3I=C3Q*(AK1+(C3Q-C2M1)*(AK2+(C3Q-C1M1)*AK3)) |
---|
944 | Z1(I)=Z1I |
---|
945 | Z2(I)=Z2I |
---|
946 | Z3(I)=Z3I |
---|
947 | F1(I)=TI11*Z1I+TI12*Z2I+TI13*Z3I |
---|
948 | F2(I)=TI21*Z1I+TI22*Z2I+TI23*Z3I |
---|
949 | F3(I)=TI31*Z1I+TI32*Z2I+TI33*Z3I |
---|
950 | END DO |
---|
951 | END IF |
---|
952 | C *** *** *** *** *** *** *** |
---|
953 | C LOOP FOR THE SIMPLIFIED NEWTON ITERATION |
---|
954 | C *** *** *** *** *** *** *** |
---|
955 | NEWT=0 |
---|
956 | FACCON=MAX(FACCON,UROUND)**0.8D0 |
---|
957 | THETA=ABS(THET) |
---|
958 | 40 CONTINUE |
---|
959 | IF (NEWT.GE.NIT) GOTO 78 |
---|
960 | C --- COMPUTE THE RIGHT-HAND SIDE |
---|
961 | DO I=1,N |
---|
962 | CONT(I)=Y(I)+Z1(I) |
---|
963 | END DO |
---|
964 | CALL FCN(N,X+C1*H,CONT,Z1,RPAR,IPAR) |
---|
965 | DO I=1,N |
---|
966 | CONT(I)=Y(I)+Z2(I) |
---|
967 | END DO |
---|
968 | CALL FCN(N,X+C2*H,CONT,Z2,RPAR,IPAR) |
---|
969 | DO I=1,N |
---|
970 | CONT(I)=Y(I)+Z3(I) |
---|
971 | END DO |
---|
972 | CALL FCN(N,XPH,CONT,Z3,RPAR,IPAR) |
---|
973 | NFCN=NFCN+3 |
---|
974 | C --- KppSolve THE LINEAR SYSTEMS |
---|
975 | DO I=1,N |
---|
976 | A1=Z1(I) |
---|
977 | A2=Z2(I) |
---|
978 | A3=Z3(I) |
---|
979 | Z1(I)=TI11*A1+TI12*A2+TI13*A3 |
---|
980 | Z2(I)=TI21*A1+TI22*A2+TI23*A3 |
---|
981 | Z3(I)=TI31*A1+TI32*A2+TI33*A3 |
---|
982 | END DO |
---|
983 | CALL SLVRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
984 | & M1,M2,NM1,FAC1,ALPHN,BETAN,E1,E2R,E2I,LDE1,Z1,Z2,Z3, |
---|
985 | & F1,F2,F3,CONT,IP1,IP2,IPHES,IER,IJOB) |
---|
986 | NSOL=NSOL+1 |
---|
987 | NEWT=NEWT+1 |
---|
988 | DYNO=0.D0 |
---|
989 | DO I=1,N |
---|
990 | DENOM=SCAL(I) |
---|
991 | DYNO=DYNO+(Z1(I)/DENOM)**2+(Z2(I)/DENOM)**2 |
---|
992 | & +(Z3(I)/DENOM)**2 |
---|
993 | END DO |
---|
994 | DYNO=DSQRT(DYNO/N3) |
---|
995 | C --- BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE |
---|
996 | IF (NEWT.GT.1.AND.NEWT.LT.NIT) THEN |
---|
997 | THQ=DYNO/DYNOLD |
---|
998 | IF (NEWT.EQ.2) THEN |
---|
999 | THETA=THQ |
---|
1000 | ELSE |
---|
1001 | THETA=SQRT(THQ*THQOLD) |
---|
1002 | END IF |
---|
1003 | THQOLD=THQ |
---|
1004 | IF (THETA.LT.0.99D0) THEN |
---|
1005 | FACCON=THETA/(1.0D0-THETA) |
---|
1006 | DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)/FNEWT |
---|
1007 | IF (DYTH.GE.1.0D0) THEN |
---|
1008 | QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH)) |
---|
1009 | HHFAC=.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT)) |
---|
1010 | H=HHFAC*H |
---|
1011 | REJECT=.TRUE. |
---|
1012 | LAST=.FALSE. |
---|
1013 | IF (CALJAC) GOTO 20 |
---|
1014 | GOTO 10 |
---|
1015 | END IF |
---|
1016 | ELSE |
---|
1017 | GOTO 78 |
---|
1018 | END IF |
---|
1019 | END IF |
---|
1020 | DYNOLD=MAX(DYNO,UROUND) |
---|
1021 | DO I=1,N |
---|
1022 | F1I=F1(I)+Z1(I) |
---|
1023 | F2I=F2(I)+Z2(I) |
---|
1024 | F3I=F3(I)+Z3(I) |
---|
1025 | F1(I)=F1I |
---|
1026 | F2(I)=F2I |
---|
1027 | F3(I)=F3I |
---|
1028 | Z1(I)=T11*F1I+T12*F2I+T13*F3I |
---|
1029 | Z2(I)=T21*F1I+T22*F2I+T23*F3I |
---|
1030 | Z3(I)=T31*F1I+ F2I |
---|
1031 | END DO |
---|
1032 | IF (FACCON*DYNO.GT.FNEWT) GOTO 40 |
---|
1033 | C --- ERROR ESTIMATION |
---|
1034 | CALL ESTRAD (N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
1035 | & H,DD1,DD2,DD3,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1, |
---|
1036 | & E1,LDE1,Z1,Z2,Z3,CONT,F1,F2,IP1,IPHES,SCAL,ERR, |
---|
1037 | & FIRST,REJECT,FAC1,RPAR,IPAR) |
---|
1038 | C --- COMPUTATION OF HNEW |
---|
1039 | C --- WE REQUIRE .2<=HNEW/H<=8. |
---|
1040 | FAC=MIN(SAFE,CFAC/(NEWT+2*NIT)) |
---|
1041 | QUOT=MAX(FACR,MIN(FACL,ERR**.25D0/FAC)) |
---|
1042 | HNEW=H/QUOT |
---|
1043 | C *** *** *** *** *** *** *** |
---|
1044 | C IS THE ERROR SMALL ENOUGH ? |
---|
1045 | C *** *** *** *** *** *** *** |
---|
1046 | IF (ERR.LT.1.D0) THEN |
---|
1047 | C --- STEP IS ACCEPTED |
---|
1048 | FIRST=.FALSE. |
---|
1049 | NACCPT=NACCPT+1 |
---|
1050 | IF (PRED) THEN |
---|
1051 | C --- PREDICTIVE CONTROLLER OF GUSTAFSSON |
---|
1052 | IF (NACCPT.GT.1) THEN |
---|
1053 | FACGUS=(HACC/H)*(ERR**2/ERRACC)**0.25D0/SAFE |
---|
1054 | FACGUS=MAX(FACR,MIN(FACL,FACGUS)) |
---|
1055 | QUOT=MAX(QUOT,FACGUS) |
---|
1056 | HNEW=H/QUOT |
---|
1057 | END IF |
---|
1058 | HACC=H |
---|
1059 | ERRACC=MAX(1.0D-2,ERR) |
---|
1060 | END IF |
---|
1061 | XOLD=X |
---|
1062 | HOLD=H |
---|
1063 | X=XPH |
---|
1064 | DO I=1,N |
---|
1065 | Y(I)=Y(I)+Z3(I) |
---|
1066 | Z2I=Z2(I) |
---|
1067 | Z1I=Z1(I) |
---|
1068 | CONT(I+N)=(Z2I-Z3(I))/C2M1 |
---|
1069 | AK=(Z1I-Z2I)/C1MC2 |
---|
1070 | ACONT3=Z1I/C1 |
---|
1071 | ACONT3=(AK-ACONT3)/C2 |
---|
1072 | CONT(I+N2)=(AK-CONT(I+N))/C1M1 |
---|
1073 | CONT(I+N3)=CONT(I+N2)-ACONT3 |
---|
1074 | END DO |
---|
1075 | IF (ITOL.EQ.0) THEN |
---|
1076 | DO I=1,N |
---|
1077 | SCAL(I)=AbsTol(1)+RelTol(1)*ABS(Y(I)) |
---|
1078 | END DO |
---|
1079 | ELSE |
---|
1080 | DO I=1,N |
---|
1081 | SCAL(I)=AbsTol(I)+RelTol(I)*ABS(Y(I)) |
---|
1082 | END DO |
---|
1083 | END IF |
---|
1084 | IF (IOUT.NE.0) THEN |
---|
1085 | NRSOL=NACCPT+1 |
---|
1086 | XSOL=X |
---|
1087 | XOSOL=XOLD |
---|
1088 | DO I=1,N |
---|
1089 | CONT(I)=Y(I) |
---|
1090 | END DO |
---|
1091 | NSOLU=N |
---|
1092 | HSOL=HOLD |
---|
1093 | CALL SOLOUT(NRSOL,XOSOL,XSOL,Y,CONT,LRC,NSOLU, |
---|
1094 | & RPAR,IPAR,IRTRN) |
---|
1095 | IF (IRTRN.LT.0) GOTO 179 |
---|
1096 | END IF |
---|
1097 | CALJAC=.FALSE. |
---|
1098 | IF (LAST) THEN |
---|
1099 | H=HOPT |
---|
1100 | IDID=1 |
---|
1101 | RETURN |
---|
1102 | END IF |
---|
1103 | CALL FCN(N,X,Y,Y0,RPAR,IPAR) |
---|
1104 | NFCN=NFCN+1 |
---|
1105 | HNEW=POSNEG*MIN(ABS(HNEW),HMAXN) |
---|
1106 | HOPT=HNEW |
---|
1107 | HOPT=MIN(H,HNEW) |
---|
1108 | IF (REJECT) HNEW=POSNEG*MIN(ABS(HNEW),ABS(H)) |
---|
1109 | REJECT=.FALSE. |
---|
1110 | IF ((X+HNEW/QUOT1-XEND)*POSNEG.GE.0.D0) THEN |
---|
1111 | H=XEND-X |
---|
1112 | LAST=.TRUE. |
---|
1113 | ELSE |
---|
1114 | QT=HNEW/H |
---|
1115 | HHFAC=H |
---|
1116 | IF (THETA.LE.THET.AND.QT.GE.QUOT1.AND.QT.LE.QUOT2) GOTO 30 |
---|
1117 | H=HNEW |
---|
1118 | END IF |
---|
1119 | HHFAC=H |
---|
1120 | IF (THETA.LE.THET) GOTO 20 |
---|
1121 | GOTO 10 |
---|
1122 | ELSE |
---|
1123 | C --- STEP IS REJECTED |
---|
1124 | REJECT=.TRUE. |
---|
1125 | LAST=.FALSE. |
---|
1126 | IF (FIRST) THEN |
---|
1127 | H=H*0.1D0 |
---|
1128 | HHFAC=0.1D0 |
---|
1129 | ELSE |
---|
1130 | HHFAC=HNEW/H |
---|
1131 | H=HNEW |
---|
1132 | END IF |
---|
1133 | IF (NACCPT.GE.1) NREJCT=NREJCT+1 |
---|
1134 | IF (CALJAC) GOTO 20 |
---|
1135 | GOTO 10 |
---|
1136 | END IF |
---|
1137 | C --- UNEXPECTED STEP-REJECTION |
---|
1138 | 78 CONTINUE |
---|
1139 | IF (IER.NE.0) THEN |
---|
1140 | NSING=NSING+1 |
---|
1141 | IF (NSING.GE.5) GOTO 176 |
---|
1142 | END IF |
---|
1143 | H=H*0.5D0 |
---|
1144 | HHFAC=0.5D0 |
---|
1145 | REJECT=.TRUE. |
---|
1146 | LAST=.FALSE. |
---|
1147 | IF (CALJAC) GOTO 20 |
---|
1148 | GOTO 10 |
---|
1149 | C --- FAIL EXIT |
---|
1150 | 176 CONTINUE |
---|
1151 | WRITE(6,979)X |
---|
1152 | WRITE(6,*) ' MATRIX IS REPEATEDLY SINGULAR, IER=',IER |
---|
1153 | IDID=-4 |
---|
1154 | RETURN |
---|
1155 | 177 CONTINUE |
---|
1156 | WRITE(6,979)X |
---|
1157 | WRITE(6,*) ' STEP SIZE T0O SMALL, H=',H |
---|
1158 | IDID=-3 |
---|
1159 | RETURN |
---|
1160 | 178 CONTINUE |
---|
1161 | WRITE(6,979)X |
---|
1162 | WRITE(6,*) ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED' |
---|
1163 | IDID=-2 |
---|
1164 | RETURN |
---|
1165 | C --- EXIT CAUSED BY SOLOUT |
---|
1166 | 179 CONTINUE |
---|
1167 | WRITE(6,979)X |
---|
1168 | 979 FORMAT(' EXIT OF RADAU5 AT X=',E18.4) |
---|
1169 | IDID=2 |
---|
1170 | RETURN |
---|
1171 | END |
---|
1172 | C |
---|
1173 | C END OF SUBROUTINE RADCOR |
---|
1174 | C |
---|
1175 | C *********************************************************** |
---|
1176 | C |
---|
1177 | KPP_REAL FUNCTION CONTR5(I,X,CONT,LRC) |
---|
1178 | C ---------------------------------------------------------- |
---|
1179 | C THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT. IT PROVIDES AN |
---|
1180 | C APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT X. |
---|
1181 | C IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR |
---|
1182 | C THE LAST SUCCESSFULLY COMPUTED STEP (BY RADAU5). |
---|
1183 | C ---------------------------------------------------------- |
---|
1184 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
1185 | DIMENSION CONT(LRC) |
---|
1186 | COMMON /CONRA5/NN,NN2,NN3,NN4,XSOL,HSOL,C2M1,C1M1 |
---|
1187 | S=(X-XSOL)/HSOL |
---|
1188 | CONTR5=CONT(I)+S*(CONT(I+NN)+(S-C2M1)*(CONT(I+NN2) |
---|
1189 | & +(S-C1M1)*CONT(I+NN3))) |
---|
1190 | RETURN |
---|
1191 | END |
---|
1192 | C |
---|
1193 | C END OF FUNCTION CONTR5 |
---|
1194 | C |
---|
1195 | C *********************************************************** |
---|
1196 | C ****************************************** |
---|
1197 | C VERSION OF SEPTEMBER 18, 1995 |
---|
1198 | C ****************************************** |
---|
1199 | C |
---|
1200 | SUBROUTINE DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
1201 | & M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES) |
---|
1202 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
1203 | DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1), |
---|
1204 | & IP1(NM1),IPHES(N) |
---|
1205 | LOGICAL CALHES |
---|
1206 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
1207 | C |
---|
1208 | GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB |
---|
1209 | C |
---|
1210 | C ----------------------------------------------------------- |
---|
1211 | C |
---|
1212 | 1 CONTINUE |
---|
1213 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
1214 | DO J=1,N |
---|
1215 | DO I=1,N |
---|
1216 | E1(I,J)=-FJAC(I,J) |
---|
1217 | END DO |
---|
1218 | E1(J,J)=E1(J,J)+FAC1 |
---|
1219 | END DO |
---|
1220 | CALL DEC (N,LDE1,E1,IP1,IER) |
---|
1221 | RETURN |
---|
1222 | C |
---|
1223 | C ----------------------------------------------------------- |
---|
1224 | C |
---|
1225 | 11 CONTINUE |
---|
1226 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1227 | DO J=1,NM1 |
---|
1228 | JM1=J+M1 |
---|
1229 | DO I=1,NM1 |
---|
1230 | E1(I,J)=-FJAC(I,JM1) |
---|
1231 | END DO |
---|
1232 | E1(J,J)=E1(J,J)+FAC1 |
---|
1233 | END DO |
---|
1234 | 45 MM=M1/M2 |
---|
1235 | DO J=1,M2 |
---|
1236 | DO I=1,NM1 |
---|
1237 | SUM=0.D0 |
---|
1238 | DO K=0,MM-1 |
---|
1239 | SUM=(SUM+FJAC(I,J+K*M2))/FAC1 |
---|
1240 | END DO |
---|
1241 | E1(I,J)=E1(I,J)-SUM |
---|
1242 | END DO |
---|
1243 | END DO |
---|
1244 | CALL DEC (NM1,LDE1,E1,IP1,IER) |
---|
1245 | RETURN |
---|
1246 | C |
---|
1247 | C ----------------------------------------------------------- |
---|
1248 | C |
---|
1249 | 2 CONTINUE |
---|
1250 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX |
---|
1251 | DO J=1,N |
---|
1252 | DO I=1,MBJAC |
---|
1253 | E1(I+MLE,J)=-FJAC(I,J) |
---|
1254 | END DO |
---|
1255 | E1(MDIAG,J)=E1(MDIAG,J)+FAC1 |
---|
1256 | END DO |
---|
1257 | CALL DECB (N,LDE1,E1,MLE,MUE,IP1,IER) |
---|
1258 | RETURN |
---|
1259 | C |
---|
1260 | C ----------------------------------------------------------- |
---|
1261 | C |
---|
1262 | 12 CONTINUE |
---|
1263 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
1264 | DO J=1,NM1 |
---|
1265 | JM1=J+M1 |
---|
1266 | DO I=1,MBJAC |
---|
1267 | E1(I+MLE,J)=-FJAC(I,JM1) |
---|
1268 | END DO |
---|
1269 | E1(MDIAG,J)=E1(MDIAG,J)+FAC1 |
---|
1270 | END DO |
---|
1271 | 46 MM=M1/M2 |
---|
1272 | DO J=1,M2 |
---|
1273 | DO I=1,MBJAC |
---|
1274 | SUM=0.D0 |
---|
1275 | DO K=0,MM-1 |
---|
1276 | SUM=(SUM+FJAC(I,J+K*M2))/FAC1 |
---|
1277 | END DO |
---|
1278 | E1(I+MLE,J)=E1(I+MLE,J)-SUM |
---|
1279 | END DO |
---|
1280 | END DO |
---|
1281 | CALL DECB (NM1,LDE1,E1,MLE,MUE,IP1,IER) |
---|
1282 | RETURN |
---|
1283 | C |
---|
1284 | C ----------------------------------------------------------- |
---|
1285 | C |
---|
1286 | 3 CONTINUE |
---|
1287 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX |
---|
1288 | DO J=1,N |
---|
1289 | DO I=1,N |
---|
1290 | E1(I,J)=-FJAC(I,J) |
---|
1291 | END DO |
---|
1292 | DO I=MAX(1,J-MUMAS),MIN(N,J+MLMAS) |
---|
1293 | E1(I,J)=E1(I,J)+FAC1*FMAS(I-J+MBDIAG,J) |
---|
1294 | END DO |
---|
1295 | END DO |
---|
1296 | CALL DEC (N,LDE1,E1,IP1,IER) |
---|
1297 | RETURN |
---|
1298 | C |
---|
1299 | C ----------------------------------------------------------- |
---|
1300 | C |
---|
1301 | 13 CONTINUE |
---|
1302 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1303 | DO J=1,NM1 |
---|
1304 | JM1=J+M1 |
---|
1305 | DO I=1,NM1 |
---|
1306 | E1(I,J)=-FJAC(I,JM1) |
---|
1307 | END DO |
---|
1308 | DO I=MAX(1,J-MUMAS),MIN(NM1,J+MLMAS) |
---|
1309 | E1(I,J)=E1(I,J)+FAC1*FMAS(I-J+MBDIAG,J) |
---|
1310 | END DO |
---|
1311 | END DO |
---|
1312 | GOTO 45 |
---|
1313 | C |
---|
1314 | C ----------------------------------------------------------- |
---|
1315 | C |
---|
1316 | 4 CONTINUE |
---|
1317 | C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX |
---|
1318 | DO J=1,N |
---|
1319 | DO I=1,MBJAC |
---|
1320 | E1(I+MLE,J)=-FJAC(I,J) |
---|
1321 | END DO |
---|
1322 | DO I=1,MBB |
---|
1323 | IB=I+MDIFF |
---|
1324 | E1(IB,J)=E1(IB,J)+FAC1*FMAS(I,J) |
---|
1325 | END DO |
---|
1326 | END DO |
---|
1327 | CALL DECB (N,LDE1,E1,MLE,MUE,IP1,IER) |
---|
1328 | RETURN |
---|
1329 | C |
---|
1330 | C ----------------------------------------------------------- |
---|
1331 | C |
---|
1332 | 14 CONTINUE |
---|
1333 | C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
1334 | DO J=1,NM1 |
---|
1335 | JM1=J+M1 |
---|
1336 | DO I=1,MBJAC |
---|
1337 | E1(I+MLE,J)=-FJAC(I,JM1) |
---|
1338 | END DO |
---|
1339 | DO I=1,MBB |
---|
1340 | IB=I+MDIFF |
---|
1341 | E1(IB,J)=E1(IB,J)+FAC1*FMAS(I,J) |
---|
1342 | END DO |
---|
1343 | END DO |
---|
1344 | GOTO 46 |
---|
1345 | C |
---|
1346 | C ----------------------------------------------------------- |
---|
1347 | C |
---|
1348 | 5 CONTINUE |
---|
1349 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX |
---|
1350 | DO J=1,N |
---|
1351 | DO I=1,N |
---|
1352 | E1(I,J)=FMAS(I,J)*FAC1-FJAC(I,J) |
---|
1353 | END DO |
---|
1354 | END DO |
---|
1355 | CALL DEC (N,LDE1,E1,IP1,IER) |
---|
1356 | RETURN |
---|
1357 | C |
---|
1358 | C ----------------------------------------------------------- |
---|
1359 | C |
---|
1360 | 15 CONTINUE |
---|
1361 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1362 | DO J=1,NM1 |
---|
1363 | JM1=J+M1 |
---|
1364 | DO I=1,NM1 |
---|
1365 | E1(I,J)=FMAS(I,J)*FAC1-FJAC(I,JM1) |
---|
1366 | END DO |
---|
1367 | END DO |
---|
1368 | GOTO 45 |
---|
1369 | C |
---|
1370 | C ----------------------------------------------------------- |
---|
1371 | C |
---|
1372 | 6 CONTINUE |
---|
1373 | C --- B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX |
---|
1374 | C --- THIS OPTION IS NOT PROVIDED |
---|
1375 | RETURN |
---|
1376 | C |
---|
1377 | C ----------------------------------------------------------- |
---|
1378 | C |
---|
1379 | 7 CONTINUE |
---|
1380 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION |
---|
1381 | IF (CALHES) CALL Elmhes (LDJAC,N,1,N,FJAC,IPHES) |
---|
1382 | CALHES=.FALSE. |
---|
1383 | DO J=1,N-1 |
---|
1384 | J1=J+1 |
---|
1385 | E1(J1,J)=-FJAC(J1,J) |
---|
1386 | END DO |
---|
1387 | DO J=1,N |
---|
1388 | DO I=1,J |
---|
1389 | E1(I,J)=-FJAC(I,J) |
---|
1390 | END DO |
---|
1391 | E1(J,J)=E1(J,J)+FAC1 |
---|
1392 | END DO |
---|
1393 | CALL DECH(N,LDE1,E1,1,IP1,IER) |
---|
1394 | RETURN |
---|
1395 | C |
---|
1396 | C ----------------------------------------------------------- |
---|
1397 | C |
---|
1398 | 55 CONTINUE |
---|
1399 | RETURN |
---|
1400 | END |
---|
1401 | C |
---|
1402 | C END OF SUBROUTINE DECOMR |
---|
1403 | C |
---|
1404 | C *********************************************************** |
---|
1405 | C |
---|
1406 | SUBROUTINE DECOMC(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
1407 | & M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,IP2,IER,IJOB) |
---|
1408 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
1409 | DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1), |
---|
1410 | & E2R(LDE1,NM1),E2I(LDE1,NM1),IP2(NM1) |
---|
1411 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
1412 | C |
---|
1413 | GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB |
---|
1414 | C |
---|
1415 | C ----------------------------------------------------------- |
---|
1416 | C |
---|
1417 | 1 CONTINUE |
---|
1418 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
1419 | DO J=1,N |
---|
1420 | DO I=1,N |
---|
1421 | E2R(I,J)=-FJAC(I,J) |
---|
1422 | E2I(I,J)=0.D0 |
---|
1423 | END DO |
---|
1424 | E2R(J,J)=E2R(J,J)+ALPHN |
---|
1425 | E2I(J,J)=BETAN |
---|
1426 | END DO |
---|
1427 | CALL DECC (N,LDE1,E2R,E2I,IP2,IER) |
---|
1428 | RETURN |
---|
1429 | C |
---|
1430 | C ----------------------------------------------------------- |
---|
1431 | C |
---|
1432 | 11 CONTINUE |
---|
1433 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1434 | DO J=1,NM1 |
---|
1435 | JM1=J+M1 |
---|
1436 | DO I=1,NM1 |
---|
1437 | E2R(I,J)=-FJAC(I,JM1) |
---|
1438 | E2I(I,J)=0.D0 |
---|
1439 | END DO |
---|
1440 | E2R(J,J)=E2R(J,J)+ALPHN |
---|
1441 | E2I(J,J)=BETAN |
---|
1442 | END DO |
---|
1443 | 45 MM=M1/M2 |
---|
1444 | ABNO=ALPHN**2+BETAN**2 |
---|
1445 | ALP=ALPHN/ABNO |
---|
1446 | BET=BETAN/ABNO |
---|
1447 | DO J=1,M2 |
---|
1448 | DO I=1,NM1 |
---|
1449 | SUMR=0.D0 |
---|
1450 | SUMI=0.D0 |
---|
1451 | DO K=0,MM-1 |
---|
1452 | SUMS=SUMR+FJAC(I,J+K*M2) |
---|
1453 | SUMR=SUMS*ALP+SUMI*BET |
---|
1454 | SUMI=SUMI*ALP-SUMS*BET |
---|
1455 | END DO |
---|
1456 | E2R(I,J)=E2R(I,J)-SUMR |
---|
1457 | E2I(I,J)=E2I(I,J)-SUMI |
---|
1458 | END DO |
---|
1459 | END DO |
---|
1460 | CALL DECC (NM1,LDE1,E2R,E2I,IP2,IER) |
---|
1461 | RETURN |
---|
1462 | C |
---|
1463 | C ----------------------------------------------------------- |
---|
1464 | C |
---|
1465 | 2 CONTINUE |
---|
1466 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX |
---|
1467 | DO J=1,N |
---|
1468 | DO I=1,MBJAC |
---|
1469 | IMLE=I+MLE |
---|
1470 | E2R(IMLE,J)=-FJAC(I,J) |
---|
1471 | E2I(IMLE,J)=0.D0 |
---|
1472 | END DO |
---|
1473 | E2R(MDIAG,J)=E2R(MDIAG,J)+ALPHN |
---|
1474 | E2I(MDIAG,J)=BETAN |
---|
1475 | END DO |
---|
1476 | CALL DECBC (N,LDE1,E2R,E2I,MLE,MUE,IP2,IER) |
---|
1477 | RETURN |
---|
1478 | C |
---|
1479 | C ----------------------------------------------------------- |
---|
1480 | C |
---|
1481 | 12 CONTINUE |
---|
1482 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
1483 | DO J=1,NM1 |
---|
1484 | JM1=J+M1 |
---|
1485 | DO I=1,MBJAC |
---|
1486 | E2R(I+MLE,J)=-FJAC(I,JM1) |
---|
1487 | E2I(I+MLE,J)=0.D0 |
---|
1488 | END DO |
---|
1489 | E2R(MDIAG,J)=E2R(MDIAG,J)+ALPHN |
---|
1490 | E2I(MDIAG,J)=E2I(MDIAG,J)+BETAN |
---|
1491 | END DO |
---|
1492 | 46 MM=M1/M2 |
---|
1493 | ABNO=ALPHN**2+BETAN**2 |
---|
1494 | ALP=ALPHN/ABNO |
---|
1495 | BET=BETAN/ABNO |
---|
1496 | DO J=1,M2 |
---|
1497 | DO I=1,MBJAC |
---|
1498 | SUMR=0.D0 |
---|
1499 | SUMI=0.D0 |
---|
1500 | DO K=0,MM-1 |
---|
1501 | SUMS=SUMR+FJAC(I,J+K*M2) |
---|
1502 | SUMR=SUMS*ALP+SUMI*BET |
---|
1503 | SUMI=SUMI*ALP-SUMS*BET |
---|
1504 | END DO |
---|
1505 | IMLE=I+MLE |
---|
1506 | E2R(IMLE,J)=E2R(IMLE,J)-SUMR |
---|
1507 | E2I(IMLE,J)=E2I(IMLE,J)-SUMI |
---|
1508 | END DO |
---|
1509 | END DO |
---|
1510 | CALL DECBC (NM1,LDE1,E2R,E2I,MLE,MUE,IP2,IER) |
---|
1511 | RETURN |
---|
1512 | C |
---|
1513 | C ----------------------------------------------------------- |
---|
1514 | C |
---|
1515 | 3 CONTINUE |
---|
1516 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX |
---|
1517 | DO J=1,N |
---|
1518 | DO I=1,N |
---|
1519 | E2R(I,J)=-FJAC(I,J) |
---|
1520 | E2I(I,J)=0.D0 |
---|
1521 | END DO |
---|
1522 | END DO |
---|
1523 | DO J=1,N |
---|
1524 | DO I=MAX(1,J-MUMAS),MIN(N,J+MLMAS) |
---|
1525 | BB=FMAS(I-J+MBDIAG,J) |
---|
1526 | E2R(I,J)=E2R(I,J)+ALPHN*BB |
---|
1527 | E2I(I,J)=BETAN*BB |
---|
1528 | END DO |
---|
1529 | END DO |
---|
1530 | CALL DECC(N,LDE1,E2R,E2I,IP2,IER) |
---|
1531 | RETURN |
---|
1532 | C |
---|
1533 | C ----------------------------------------------------------- |
---|
1534 | C |
---|
1535 | 13 CONTINUE |
---|
1536 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1537 | DO J=1,NM1 |
---|
1538 | JM1=J+M1 |
---|
1539 | DO I=1,NM1 |
---|
1540 | E2R(I,J)=-FJAC(I,JM1) |
---|
1541 | E2I(I,J)=0.D0 |
---|
1542 | END DO |
---|
1543 | DO I=MAX(1,J-MUMAS),MIN(NM1,J+MLMAS) |
---|
1544 | FFMA=FMAS(I-J+MBDIAG,J) |
---|
1545 | E2R(I,J)=E2R(I,J)+ALPHN*FFMA |
---|
1546 | E2I(I,J)=E2I(I,J)+BETAN*FFMA |
---|
1547 | END DO |
---|
1548 | END DO |
---|
1549 | GOTO 45 |
---|
1550 | C |
---|
1551 | C ----------------------------------------------------------- |
---|
1552 | C |
---|
1553 | 4 CONTINUE |
---|
1554 | C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX |
---|
1555 | DO J=1,N |
---|
1556 | DO I=1,MBJAC |
---|
1557 | IMLE=I+MLE |
---|
1558 | E2R(IMLE,J)=-FJAC(I,J) |
---|
1559 | E2I(IMLE,J)=0.D0 |
---|
1560 | END DO |
---|
1561 | DO I=MAX(1,MUMAS+2-J),MIN(MBB,MUMAS+1-J+N) |
---|
1562 | IB=I+MDIFF |
---|
1563 | BB=FMAS(I,J) |
---|
1564 | E2R(IB,J)=E2R(IB,J)+ALPHN*BB |
---|
1565 | E2I(IB,J)=BETAN*BB |
---|
1566 | END DO |
---|
1567 | END DO |
---|
1568 | CALL DECBC (N,LDE1,E2R,E2I,MLE,MUE,IP2,IER) |
---|
1569 | RETURN |
---|
1570 | C |
---|
1571 | C ----------------------------------------------------------- |
---|
1572 | C |
---|
1573 | 14 CONTINUE |
---|
1574 | C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
1575 | DO J=1,NM1 |
---|
1576 | JM1=J+M1 |
---|
1577 | DO I=1,MBJAC |
---|
1578 | E2R(I+MLE,J)=-FJAC(I,JM1) |
---|
1579 | E2I(I+MLE,J)=0.D0 |
---|
1580 | END DO |
---|
1581 | DO I=1,MBB |
---|
1582 | IB=I+MDIFF |
---|
1583 | FFMA=FMAS(I,J) |
---|
1584 | E2R(IB,J)=E2R(IB,J)+ALPHN*FFMA |
---|
1585 | E2I(IB,J)=E2I(IB,J)+BETAN*FFMA |
---|
1586 | END DO |
---|
1587 | END DO |
---|
1588 | GOTO 46 |
---|
1589 | C |
---|
1590 | C ----------------------------------------------------------- |
---|
1591 | C |
---|
1592 | 5 CONTINUE |
---|
1593 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX |
---|
1594 | DO J=1,N |
---|
1595 | DO I=1,N |
---|
1596 | BB=FMAS(I,J) |
---|
1597 | E2R(I,J)=BB*ALPHN-FJAC(I,J) |
---|
1598 | E2I(I,J)=BB*BETAN |
---|
1599 | END DO |
---|
1600 | END DO |
---|
1601 | CALL DECC(N,LDE1,E2R,E2I,IP2,IER) |
---|
1602 | RETURN |
---|
1603 | C |
---|
1604 | C ----------------------------------------------------------- |
---|
1605 | C |
---|
1606 | 15 CONTINUE |
---|
1607 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1608 | DO J=1,NM1 |
---|
1609 | JM1=J+M1 |
---|
1610 | DO I=1,NM1 |
---|
1611 | E2R(I,J)=ALPHN*FMAS(I,J)-FJAC(I,JM1) |
---|
1612 | E2I(I,J)=BETAN*FMAS(I,J) |
---|
1613 | END DO |
---|
1614 | END DO |
---|
1615 | GOTO 45 |
---|
1616 | C |
---|
1617 | C ----------------------------------------------------------- |
---|
1618 | C |
---|
1619 | 6 CONTINUE |
---|
1620 | C --- B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX |
---|
1621 | C --- THIS OPTION IS NOT PROVIDED |
---|
1622 | RETURN |
---|
1623 | C |
---|
1624 | C ----------------------------------------------------------- |
---|
1625 | C |
---|
1626 | 7 CONTINUE |
---|
1627 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION |
---|
1628 | DO J=1,N-1 |
---|
1629 | J1=J+1 |
---|
1630 | E2R(J1,J)=-FJAC(J1,J) |
---|
1631 | E2I(J1,J)=0.D0 |
---|
1632 | END DO |
---|
1633 | DO J=1,N |
---|
1634 | DO I=1,J |
---|
1635 | E2I(I,J)=0.D0 |
---|
1636 | E2R(I,J)=-FJAC(I,J) |
---|
1637 | END DO |
---|
1638 | E2R(J,J)=E2R(J,J)+ALPHN |
---|
1639 | E2I(J,J)=BETAN |
---|
1640 | END DO |
---|
1641 | CALL DECHC(N,LDE1,E2R,E2I,1,IP2,IER) |
---|
1642 | RETURN |
---|
1643 | C |
---|
1644 | C ----------------------------------------------------------- |
---|
1645 | C |
---|
1646 | 55 CONTINUE |
---|
1647 | RETURN |
---|
1648 | END |
---|
1649 | C |
---|
1650 | C END OF SUBROUTINE DECOMC |
---|
1651 | C |
---|
1652 | C *********************************************************** |
---|
1653 | C |
---|
1654 | SUBROUTINE SLVRAR(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
1655 | & M1,M2,NM1,FAC1,E1,LDE1,Z1,F1,IP1,IPHES,IER,IJOB) |
---|
1656 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
1657 | DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1), |
---|
1658 | & IP1(NM1),IPHES(N),Z1(N),F1(N) |
---|
1659 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
1660 | C |
---|
1661 | GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,13,15), IJOB |
---|
1662 | C |
---|
1663 | C ----------------------------------------------------------- |
---|
1664 | C |
---|
1665 | 1 CONTINUE |
---|
1666 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
1667 | DO I=1,N |
---|
1668 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
1669 | END DO |
---|
1670 | CALL SOL (N,LDE1,E1,Z1,IP1) |
---|
1671 | RETURN |
---|
1672 | C |
---|
1673 | C ----------------------------------------------------------- |
---|
1674 | C |
---|
1675 | 11 CONTINUE |
---|
1676 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1677 | DO I=1,N |
---|
1678 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
1679 | END DO |
---|
1680 | 48 CONTINUE |
---|
1681 | MM=M1/M2 |
---|
1682 | DO J=1,M2 |
---|
1683 | SUM1=0.D0 |
---|
1684 | DO K=MM-1,0,-1 |
---|
1685 | JKM=J+K*M2 |
---|
1686 | SUM1=(Z1(JKM)+SUM1)/FAC1 |
---|
1687 | DO I=1,NM1 |
---|
1688 | IM1=I+M1 |
---|
1689 | Z1(IM1)=Z1(IM1)+FJAC(I,JKM)*SUM1 |
---|
1690 | END DO |
---|
1691 | END DO |
---|
1692 | END DO |
---|
1693 | CALL SOL (NM1,LDE1,E1,Z1(M1+1),IP1) |
---|
1694 | 49 CONTINUE |
---|
1695 | DO I=M1,1,-1 |
---|
1696 | Z1(I)=(Z1(I)+Z1(M2+I))/FAC1 |
---|
1697 | END DO |
---|
1698 | RETURN |
---|
1699 | C |
---|
1700 | C ----------------------------------------------------------- |
---|
1701 | C |
---|
1702 | 2 CONTINUE |
---|
1703 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX |
---|
1704 | DO I=1,N |
---|
1705 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
1706 | END DO |
---|
1707 | CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1) |
---|
1708 | RETURN |
---|
1709 | C |
---|
1710 | C ----------------------------------------------------------- |
---|
1711 | C |
---|
1712 | 12 CONTINUE |
---|
1713 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
1714 | DO I=1,N |
---|
1715 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
1716 | END DO |
---|
1717 | 45 CONTINUE |
---|
1718 | MM=M1/M2 |
---|
1719 | DO J=1,M2 |
---|
1720 | SUM1=0.D0 |
---|
1721 | DO K=MM-1,0,-1 |
---|
1722 | JKM=J+K*M2 |
---|
1723 | SUM1=(Z1(JKM)+SUM1)/FAC1 |
---|
1724 | DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC) |
---|
1725 | IM1=I+M1 |
---|
1726 | Z1(IM1)=Z1(IM1)+FJAC(I+MUJAC+1-J,JKM)*SUM1 |
---|
1727 | END DO |
---|
1728 | END DO |
---|
1729 | END DO |
---|
1730 | CALL SOLB (NM1,LDE1,E1,MLE,MUE,Z1(M1+1),IP1) |
---|
1731 | GOTO 49 |
---|
1732 | C |
---|
1733 | C ----------------------------------------------------------- |
---|
1734 | C |
---|
1735 | 3 CONTINUE |
---|
1736 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX |
---|
1737 | DO I=1,N |
---|
1738 | S1=0.0D0 |
---|
1739 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
1740 | S1=S1-FMAS(I-J+MBDIAG,J)*F1(J) |
---|
1741 | END DO |
---|
1742 | Z1(I)=Z1(I)+S1*FAC1 |
---|
1743 | END DO |
---|
1744 | CALL SOL (N,LDE1,E1,Z1,IP1) |
---|
1745 | RETURN |
---|
1746 | C |
---|
1747 | C ----------------------------------------------------------- |
---|
1748 | C |
---|
1749 | 13 CONTINUE |
---|
1750 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1751 | DO I=1,M1 |
---|
1752 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
1753 | END DO |
---|
1754 | DO I=1,NM1 |
---|
1755 | IM1=I+M1 |
---|
1756 | S1=0.0D0 |
---|
1757 | DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS) |
---|
1758 | S1=S1-FMAS(I-J+MBDIAG,J)*F1(J+M1) |
---|
1759 | END DO |
---|
1760 | Z1(IM1)=Z1(IM1)+S1*FAC1 |
---|
1761 | END DO |
---|
1762 | IF (IJOB.EQ.14) GOTO 45 |
---|
1763 | GOTO 48 |
---|
1764 | C |
---|
1765 | C ----------------------------------------------------------- |
---|
1766 | C |
---|
1767 | 4 CONTINUE |
---|
1768 | C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX |
---|
1769 | DO I=1,N |
---|
1770 | S1=0.0D0 |
---|
1771 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
1772 | S1=S1-FMAS(I-J+MBDIAG,J)*F1(J) |
---|
1773 | END DO |
---|
1774 | Z1(I)=Z1(I)+S1*FAC1 |
---|
1775 | END DO |
---|
1776 | CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1) |
---|
1777 | RETURN |
---|
1778 | C |
---|
1779 | C ----------------------------------------------------------- |
---|
1780 | C |
---|
1781 | 5 CONTINUE |
---|
1782 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX |
---|
1783 | DO I=1,N |
---|
1784 | S1=0.0D0 |
---|
1785 | DO J=1,N |
---|
1786 | S1=S1-FMAS(I,J)*F1(J) |
---|
1787 | END DO |
---|
1788 | Z1(I)=Z1(I)+S1*FAC1 |
---|
1789 | END DO |
---|
1790 | CALL SOL (N,LDE1,E1,Z1,IP1) |
---|
1791 | RETURN |
---|
1792 | C |
---|
1793 | C ----------------------------------------------------------- |
---|
1794 | C |
---|
1795 | 15 CONTINUE |
---|
1796 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1797 | DO I=1,M1 |
---|
1798 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
1799 | END DO |
---|
1800 | DO I=1,NM1 |
---|
1801 | IM1=I+M1 |
---|
1802 | S1=0.0D0 |
---|
1803 | DO J=1,NM1 |
---|
1804 | S1=S1-FMAS(I,J)*F1(J+M1) |
---|
1805 | END DO |
---|
1806 | Z1(IM1)=Z1(IM1)+S1*FAC1 |
---|
1807 | END DO |
---|
1808 | GOTO 48 |
---|
1809 | C |
---|
1810 | C ----------------------------------------------------------- |
---|
1811 | C |
---|
1812 | 6 CONTINUE |
---|
1813 | C --- B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX |
---|
1814 | C --- THIS OPTION IS NOT PROVIDED |
---|
1815 | RETURN |
---|
1816 | C |
---|
1817 | C ----------------------------------------------------------- |
---|
1818 | C |
---|
1819 | 7 CONTINUE |
---|
1820 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION |
---|
1821 | DO I=1,N |
---|
1822 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
1823 | END DO |
---|
1824 | DO MM=N-2,1,-1 |
---|
1825 | MP=N-MM |
---|
1826 | MP1=MP-1 |
---|
1827 | I=IPHES(MP) |
---|
1828 | IF (I.EQ.MP) GOTO 746 |
---|
1829 | ZSAFE=Z1(MP) |
---|
1830 | Z1(MP)=Z1(I) |
---|
1831 | Z1(I)=ZSAFE |
---|
1832 | 746 CONTINUE |
---|
1833 | DO I=MP+1,N |
---|
1834 | Z1(I)=Z1(I)-FJAC(I,MP1)*Z1(MP) |
---|
1835 | END DO |
---|
1836 | END DO |
---|
1837 | CALL SOLH(N,LDE1,E1,1,Z1,IP1) |
---|
1838 | DO MM=1,N-2 |
---|
1839 | MP=N-MM |
---|
1840 | MP1=MP-1 |
---|
1841 | DO I=MP+1,N |
---|
1842 | Z1(I)=Z1(I)+FJAC(I,MP1)*Z1(MP) |
---|
1843 | END DO |
---|
1844 | I=IPHES(MP) |
---|
1845 | IF (I.EQ.MP) GOTO 750 |
---|
1846 | ZSAFE=Z1(MP) |
---|
1847 | Z1(MP)=Z1(I) |
---|
1848 | Z1(I)=ZSAFE |
---|
1849 | 750 CONTINUE |
---|
1850 | END DO |
---|
1851 | RETURN |
---|
1852 | C |
---|
1853 | C ----------------------------------------------------------- |
---|
1854 | C |
---|
1855 | 55 CONTINUE |
---|
1856 | RETURN |
---|
1857 | END |
---|
1858 | C |
---|
1859 | C END OF SUBROUTINE SLVRAR |
---|
1860 | C |
---|
1861 | C *********************************************************** |
---|
1862 | C |
---|
1863 | SUBROUTINE SLVRAI(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
1864 | & M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,Z2,Z3, |
---|
1865 | & F2,F3,CONT,IP2,IPHES,IER,IJOB) |
---|
1866 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
1867 | DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1), |
---|
1868 | & IP2(NM1),IPHES(N),Z2(N),Z3(N),F2(N),F3(N) |
---|
1869 | DIMENSION E2R(LDE1,NM1),E2I(LDE1,NM1) |
---|
1870 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
1871 | C |
---|
1872 | GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,13,15), IJOB |
---|
1873 | C |
---|
1874 | C ----------------------------------------------------------- |
---|
1875 | C |
---|
1876 | 1 CONTINUE |
---|
1877 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
1878 | DO I=1,N |
---|
1879 | S2=-F2(I) |
---|
1880 | S3=-F3(I) |
---|
1881 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
1882 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
1883 | END DO |
---|
1884 | CALL SOLC (N,LDE1,E2R,E2I,Z2,Z3,IP2) |
---|
1885 | RETURN |
---|
1886 | C |
---|
1887 | C ----------------------------------------------------------- |
---|
1888 | C |
---|
1889 | 11 CONTINUE |
---|
1890 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1891 | DO I=1,N |
---|
1892 | S2=-F2(I) |
---|
1893 | S3=-F3(I) |
---|
1894 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
1895 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
1896 | END DO |
---|
1897 | 48 ABNO=ALPHN**2+BETAN**2 |
---|
1898 | MM=M1/M2 |
---|
1899 | DO J=1,M2 |
---|
1900 | SUM2=0.D0 |
---|
1901 | SUM3=0.D0 |
---|
1902 | DO K=MM-1,0,-1 |
---|
1903 | JKM=J+K*M2 |
---|
1904 | SUMH=(Z2(JKM)+SUM2)/ABNO |
---|
1905 | SUM3=(Z3(JKM)+SUM3)/ABNO |
---|
1906 | SUM2=SUMH*ALPHN+SUM3*BETAN |
---|
1907 | SUM3=SUM3*ALPHN-SUMH*BETAN |
---|
1908 | DO I=1,NM1 |
---|
1909 | IM1=I+M1 |
---|
1910 | Z2(IM1)=Z2(IM1)+FJAC(I,JKM)*SUM2 |
---|
1911 | Z3(IM1)=Z3(IM1)+FJAC(I,JKM)*SUM3 |
---|
1912 | END DO |
---|
1913 | END DO |
---|
1914 | END DO |
---|
1915 | CALL SOLC (NM1,LDE1,E2R,E2I,Z2(M1+1),Z3(M1+1),IP2) |
---|
1916 | 49 CONTINUE |
---|
1917 | DO I=M1,1,-1 |
---|
1918 | MPI=M2+I |
---|
1919 | Z2I=Z2(I)+Z2(MPI) |
---|
1920 | Z3I=Z3(I)+Z3(MPI) |
---|
1921 | Z3(I)=(Z3I*ALPHN-Z2I*BETAN)/ABNO |
---|
1922 | Z2(I)=(Z2I*ALPHN+Z3I*BETAN)/ABNO |
---|
1923 | END DO |
---|
1924 | RETURN |
---|
1925 | C |
---|
1926 | C ----------------------------------------------------------- |
---|
1927 | C |
---|
1928 | 2 CONTINUE |
---|
1929 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX |
---|
1930 | DO I=1,N |
---|
1931 | S2=-F2(I) |
---|
1932 | S3=-F3(I) |
---|
1933 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
1934 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
1935 | END DO |
---|
1936 | CALL SOLBC (N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2) |
---|
1937 | RETURN |
---|
1938 | C |
---|
1939 | C ----------------------------------------------------------- |
---|
1940 | C |
---|
1941 | 12 CONTINUE |
---|
1942 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
1943 | DO I=1,N |
---|
1944 | S2=-F2(I) |
---|
1945 | S3=-F3(I) |
---|
1946 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
1947 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
1948 | END DO |
---|
1949 | 45 ABNO=ALPHN**2+BETAN**2 |
---|
1950 | MM=M1/M2 |
---|
1951 | DO J=1,M2 |
---|
1952 | SUM2=0.D0 |
---|
1953 | SUM3=0.D0 |
---|
1954 | DO K=MM-1,0,-1 |
---|
1955 | JKM=J+K*M2 |
---|
1956 | SUMH=(Z2(JKM)+SUM2)/ABNO |
---|
1957 | SUM3=(Z3(JKM)+SUM3)/ABNO |
---|
1958 | SUM2=SUMH*ALPHN+SUM3*BETAN |
---|
1959 | SUM3=SUM3*ALPHN-SUMH*BETAN |
---|
1960 | DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC) |
---|
1961 | IM1=I+M1 |
---|
1962 | IIMU=I+MUJAC+1-J |
---|
1963 | Z2(IM1)=Z2(IM1)+FJAC(IIMU,JKM)*SUM2 |
---|
1964 | Z3(IM1)=Z3(IM1)+FJAC(IIMU,JKM)*SUM3 |
---|
1965 | END DO |
---|
1966 | END DO |
---|
1967 | END DO |
---|
1968 | CALL SOLBC (NM1,LDE1,E2R,E2I,MLE,MUE,Z2(M1+1),Z3(M1+1),IP2) |
---|
1969 | GOTO 49 |
---|
1970 | C |
---|
1971 | C ----------------------------------------------------------- |
---|
1972 | C |
---|
1973 | 3 CONTINUE |
---|
1974 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX |
---|
1975 | DO I=1,N |
---|
1976 | S2=0.0D0 |
---|
1977 | S3=0.0D0 |
---|
1978 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
1979 | BB=FMAS(I-J+MBDIAG,J) |
---|
1980 | S2=S2-BB*F2(J) |
---|
1981 | S3=S3-BB*F3(J) |
---|
1982 | END DO |
---|
1983 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
1984 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
1985 | END DO |
---|
1986 | CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2) |
---|
1987 | RETURN |
---|
1988 | C |
---|
1989 | C ----------------------------------------------------------- |
---|
1990 | C |
---|
1991 | 13 CONTINUE |
---|
1992 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
1993 | DO I=1,M1 |
---|
1994 | S2=-F2(I) |
---|
1995 | S3=-F3(I) |
---|
1996 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
1997 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
1998 | END DO |
---|
1999 | DO I=1,NM1 |
---|
2000 | IM1=I+M1 |
---|
2001 | S2=0.0D0 |
---|
2002 | S3=0.0D0 |
---|
2003 | DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS) |
---|
2004 | JM1=J+M1 |
---|
2005 | BB=FMAS(I-J+MBDIAG,J) |
---|
2006 | S2=S2-BB*F2(JM1) |
---|
2007 | S3=S3-BB*F3(JM1) |
---|
2008 | END DO |
---|
2009 | Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN |
---|
2010 | Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN |
---|
2011 | END DO |
---|
2012 | IF (IJOB.EQ.14) GOTO 45 |
---|
2013 | GOTO 48 |
---|
2014 | C |
---|
2015 | C ----------------------------------------------------------- |
---|
2016 | C |
---|
2017 | 4 CONTINUE |
---|
2018 | C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX |
---|
2019 | DO I=1,N |
---|
2020 | S2=0.0D0 |
---|
2021 | S3=0.0D0 |
---|
2022 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
2023 | BB=FMAS(I-J+MBDIAG,J) |
---|
2024 | S2=S2-BB*F2(J) |
---|
2025 | S3=S3-BB*F3(J) |
---|
2026 | END DO |
---|
2027 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2028 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2029 | END DO |
---|
2030 | CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2) |
---|
2031 | RETURN |
---|
2032 | C |
---|
2033 | C ----------------------------------------------------------- |
---|
2034 | C |
---|
2035 | 5 CONTINUE |
---|
2036 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX |
---|
2037 | DO I=1,N |
---|
2038 | S2=0.0D0 |
---|
2039 | S3=0.0D0 |
---|
2040 | DO J=1,N |
---|
2041 | BB=FMAS(I,J) |
---|
2042 | S2=S2-BB*F2(J) |
---|
2043 | S3=S3-BB*F3(J) |
---|
2044 | END DO |
---|
2045 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2046 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2047 | END DO |
---|
2048 | CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2) |
---|
2049 | RETURN |
---|
2050 | C |
---|
2051 | C ----------------------------------------------------------- |
---|
2052 | C |
---|
2053 | 15 CONTINUE |
---|
2054 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
2055 | DO I=1,M1 |
---|
2056 | S2=-F2(I) |
---|
2057 | S3=-F3(I) |
---|
2058 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2059 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2060 | END DO |
---|
2061 | DO I=1,NM1 |
---|
2062 | IM1=I+M1 |
---|
2063 | S2=0.0D0 |
---|
2064 | S3=0.0D0 |
---|
2065 | DO J=1,NM1 |
---|
2066 | JM1=J+M1 |
---|
2067 | BB=FMAS(I,J) |
---|
2068 | S2=S2-BB*F2(JM1) |
---|
2069 | S3=S3-BB*F3(JM1) |
---|
2070 | END DO |
---|
2071 | Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN |
---|
2072 | Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN |
---|
2073 | END DO |
---|
2074 | GOTO 48 |
---|
2075 | C |
---|
2076 | C ----------------------------------------------------------- |
---|
2077 | C |
---|
2078 | 6 CONTINUE |
---|
2079 | C --- B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX |
---|
2080 | C --- THIS OPTION IS NOT PROVIDED |
---|
2081 | RETURN |
---|
2082 | C |
---|
2083 | C ----------------------------------------------------------- |
---|
2084 | C |
---|
2085 | 7 CONTINUE |
---|
2086 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION |
---|
2087 | DO I=1,N |
---|
2088 | S2=-F2(I) |
---|
2089 | S3=-F3(I) |
---|
2090 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2091 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2092 | END DO |
---|
2093 | DO MM=N-2,1,-1 |
---|
2094 | MP=N-MM |
---|
2095 | MP1=MP-1 |
---|
2096 | I=IPHES(MP) |
---|
2097 | IF (I.EQ.MP) GOTO 746 |
---|
2098 | ZSAFE=Z2(MP) |
---|
2099 | Z2(MP)=Z2(I) |
---|
2100 | Z2(I)=ZSAFE |
---|
2101 | ZSAFE=Z3(MP) |
---|
2102 | Z3(MP)=Z3(I) |
---|
2103 | Z3(I)=ZSAFE |
---|
2104 | 746 CONTINUE |
---|
2105 | DO I=MP+1,N |
---|
2106 | E1IMP=FJAC(I,MP1) |
---|
2107 | Z2(I)=Z2(I)-E1IMP*Z2(MP) |
---|
2108 | Z3(I)=Z3(I)-E1IMP*Z3(MP) |
---|
2109 | END DO |
---|
2110 | END DO |
---|
2111 | CALL SOLHC(N,LDE1,E2R,E2I,1,Z2,Z3,IP2) |
---|
2112 | DO MM=1,N-2 |
---|
2113 | MP=N-MM |
---|
2114 | MP1=MP-1 |
---|
2115 | DO I=MP+1,N |
---|
2116 | E1IMP=FJAC(I,MP1) |
---|
2117 | Z2(I)=Z2(I)+E1IMP*Z2(MP) |
---|
2118 | Z3(I)=Z3(I)+E1IMP*Z3(MP) |
---|
2119 | END DO |
---|
2120 | I=IPHES(MP) |
---|
2121 | IF (I.EQ.MP) GOTO 750 |
---|
2122 | ZSAFE=Z2(MP) |
---|
2123 | Z2(MP)=Z2(I) |
---|
2124 | Z2(I)=ZSAFE |
---|
2125 | ZSAFE=Z3(MP) |
---|
2126 | Z3(MP)=Z3(I) |
---|
2127 | Z3(I)=ZSAFE |
---|
2128 | 750 CONTINUE |
---|
2129 | END DO |
---|
2130 | RETURN |
---|
2131 | C |
---|
2132 | C ----------------------------------------------------------- |
---|
2133 | C |
---|
2134 | 55 CONTINUE |
---|
2135 | RETURN |
---|
2136 | END |
---|
2137 | C |
---|
2138 | C END OF SUBROUTINE SLVRAI |
---|
2139 | C |
---|
2140 | C *********************************************************** |
---|
2141 | C |
---|
2142 | SUBROUTINE SLVRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
2143 | & M1,M2,NM1,FAC1,ALPHN,BETAN,E1,E2R,E2I,LDE1,Z1,Z2,Z3, |
---|
2144 | & F1,F2,F3,CONT,IP1,IP2,IPHES,IER,IJOB) |
---|
2145 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
2146 | DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1), |
---|
2147 | & E2R(LDE1,NM1),E2I(LDE1,NM1),IP1(NM1),IP2(NM1), |
---|
2148 | & IPHES(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N),F3(N) |
---|
2149 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
2150 | C |
---|
2151 | GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,13,15), IJOB |
---|
2152 | C |
---|
2153 | C ----------------------------------------------------------- |
---|
2154 | C |
---|
2155 | 1 CONTINUE |
---|
2156 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
2157 | DO I=1,N |
---|
2158 | S2=-F2(I) |
---|
2159 | S3=-F3(I) |
---|
2160 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
2161 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2162 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2163 | END DO |
---|
2164 | CALL SOL (N,LDE1,E1,Z1,IP1) |
---|
2165 | CALL SOLC (N,LDE1,E2R,E2I,Z2,Z3,IP2) |
---|
2166 | RETURN |
---|
2167 | C |
---|
2168 | C ----------------------------------------------------------- |
---|
2169 | C |
---|
2170 | 11 CONTINUE |
---|
2171 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
2172 | DO I=1,N |
---|
2173 | S2=-F2(I) |
---|
2174 | S3=-F3(I) |
---|
2175 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
2176 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2177 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2178 | END DO |
---|
2179 | 48 ABNO=ALPHN**2+BETAN**2 |
---|
2180 | MM=M1/M2 |
---|
2181 | DO J=1,M2 |
---|
2182 | SUM1=0.D0 |
---|
2183 | SUM2=0.D0 |
---|
2184 | SUM3=0.D0 |
---|
2185 | DO K=MM-1,0,-1 |
---|
2186 | JKM=J+K*M2 |
---|
2187 | SUM1=(Z1(JKM)+SUM1)/FAC1 |
---|
2188 | SUMH=(Z2(JKM)+SUM2)/ABNO |
---|
2189 | SUM3=(Z3(JKM)+SUM3)/ABNO |
---|
2190 | SUM2=SUMH*ALPHN+SUM3*BETAN |
---|
2191 | SUM3=SUM3*ALPHN-SUMH*BETAN |
---|
2192 | DO I=1,NM1 |
---|
2193 | IM1=I+M1 |
---|
2194 | Z1(IM1)=Z1(IM1)+FJAC(I,JKM)*SUM1 |
---|
2195 | Z2(IM1)=Z2(IM1)+FJAC(I,JKM)*SUM2 |
---|
2196 | Z3(IM1)=Z3(IM1)+FJAC(I,JKM)*SUM3 |
---|
2197 | END DO |
---|
2198 | END DO |
---|
2199 | END DO |
---|
2200 | CALL SOL (NM1,LDE1,E1,Z1(M1+1),IP1) |
---|
2201 | CALL SOLC (NM1,LDE1,E2R,E2I,Z2(M1+1),Z3(M1+1),IP2) |
---|
2202 | 49 CONTINUE |
---|
2203 | DO I=M1,1,-1 |
---|
2204 | MPI=M2+I |
---|
2205 | Z1(I)=(Z1(I)+Z1(MPI))/FAC1 |
---|
2206 | Z2I=Z2(I)+Z2(MPI) |
---|
2207 | Z3I=Z3(I)+Z3(MPI) |
---|
2208 | Z3(I)=(Z3I*ALPHN-Z2I*BETAN)/ABNO |
---|
2209 | Z2(I)=(Z2I*ALPHN+Z3I*BETAN)/ABNO |
---|
2210 | END DO |
---|
2211 | RETURN |
---|
2212 | C |
---|
2213 | C ----------------------------------------------------------- |
---|
2214 | C |
---|
2215 | 2 CONTINUE |
---|
2216 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX |
---|
2217 | DO I=1,N |
---|
2218 | S2=-F2(I) |
---|
2219 | S3=-F3(I) |
---|
2220 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
2221 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2222 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2223 | END DO |
---|
2224 | CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1) |
---|
2225 | CALL SOLBC (N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2) |
---|
2226 | RETURN |
---|
2227 | C |
---|
2228 | C ----------------------------------------------------------- |
---|
2229 | C |
---|
2230 | 12 CONTINUE |
---|
2231 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
2232 | DO I=1,N |
---|
2233 | S2=-F2(I) |
---|
2234 | S3=-F3(I) |
---|
2235 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
2236 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2237 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2238 | END DO |
---|
2239 | 45 ABNO=ALPHN**2+BETAN**2 |
---|
2240 | MM=M1/M2 |
---|
2241 | DO J=1,M2 |
---|
2242 | SUM1=0.D0 |
---|
2243 | SUM2=0.D0 |
---|
2244 | SUM3=0.D0 |
---|
2245 | DO K=MM-1,0,-1 |
---|
2246 | JKM=J+K*M2 |
---|
2247 | SUM1=(Z1(JKM)+SUM1)/FAC1 |
---|
2248 | SUMH=(Z2(JKM)+SUM2)/ABNO |
---|
2249 | SUM3=(Z3(JKM)+SUM3)/ABNO |
---|
2250 | SUM2=SUMH*ALPHN+SUM3*BETAN |
---|
2251 | SUM3=SUM3*ALPHN-SUMH*BETAN |
---|
2252 | DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC) |
---|
2253 | IM1=I+M1 |
---|
2254 | FFJA=FJAC(I+MUJAC+1-J,JKM) |
---|
2255 | Z1(IM1)=Z1(IM1)+FFJA*SUM1 |
---|
2256 | Z2(IM1)=Z2(IM1)+FFJA*SUM2 |
---|
2257 | Z3(IM1)=Z3(IM1)+FFJA*SUM3 |
---|
2258 | END DO |
---|
2259 | END DO |
---|
2260 | END DO |
---|
2261 | CALL SOLB (NM1,LDE1,E1,MLE,MUE,Z1(M1+1),IP1) |
---|
2262 | CALL SOLBC (NM1,LDE1,E2R,E2I,MLE,MUE,Z2(M1+1),Z3(M1+1),IP2) |
---|
2263 | GOTO 49 |
---|
2264 | C |
---|
2265 | C ----------------------------------------------------------- |
---|
2266 | C |
---|
2267 | 3 CONTINUE |
---|
2268 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX |
---|
2269 | DO I=1,N |
---|
2270 | S1=0.0D0 |
---|
2271 | S2=0.0D0 |
---|
2272 | S3=0.0D0 |
---|
2273 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
2274 | BB=FMAS(I-J+MBDIAG,J) |
---|
2275 | S1=S1-BB*F1(J) |
---|
2276 | S2=S2-BB*F2(J) |
---|
2277 | S3=S3-BB*F3(J) |
---|
2278 | END DO |
---|
2279 | Z1(I)=Z1(I)+S1*FAC1 |
---|
2280 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2281 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2282 | END DO |
---|
2283 | CALL SOL (N,LDE1,E1,Z1,IP1) |
---|
2284 | CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2) |
---|
2285 | RETURN |
---|
2286 | C |
---|
2287 | C ----------------------------------------------------------- |
---|
2288 | C |
---|
2289 | 13 CONTINUE |
---|
2290 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
2291 | DO I=1,M1 |
---|
2292 | S2=-F2(I) |
---|
2293 | S3=-F3(I) |
---|
2294 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
2295 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2296 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2297 | END DO |
---|
2298 | DO I=1,NM1 |
---|
2299 | IM1=I+M1 |
---|
2300 | S1=0.0D0 |
---|
2301 | S2=0.0D0 |
---|
2302 | S3=0.0D0 |
---|
2303 | J1B=MAX(1,I-MLMAS) |
---|
2304 | J2B=MIN(NM1,I+MUMAS) |
---|
2305 | DO J=J1B,J2B |
---|
2306 | JM1=J+M1 |
---|
2307 | BB=FMAS(I-J+MBDIAG,J) |
---|
2308 | S1=S1-BB*F1(JM1) |
---|
2309 | S2=S2-BB*F2(JM1) |
---|
2310 | S3=S3-BB*F3(JM1) |
---|
2311 | END DO |
---|
2312 | Z1(IM1)=Z1(IM1)+S1*FAC1 |
---|
2313 | Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN |
---|
2314 | Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN |
---|
2315 | END DO |
---|
2316 | IF (IJOB.EQ.14) GOTO 45 |
---|
2317 | GOTO 48 |
---|
2318 | C |
---|
2319 | C ----------------------------------------------------------- |
---|
2320 | C |
---|
2321 | 4 CONTINUE |
---|
2322 | C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX |
---|
2323 | DO I=1,N |
---|
2324 | S1=0.0D0 |
---|
2325 | S2=0.0D0 |
---|
2326 | S3=0.0D0 |
---|
2327 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
2328 | BB=FMAS(I-J+MBDIAG,J) |
---|
2329 | S1=S1-BB*F1(J) |
---|
2330 | S2=S2-BB*F2(J) |
---|
2331 | S3=S3-BB*F3(J) |
---|
2332 | END DO |
---|
2333 | Z1(I)=Z1(I)+S1*FAC1 |
---|
2334 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2335 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2336 | END DO |
---|
2337 | CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1) |
---|
2338 | CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2) |
---|
2339 | RETURN |
---|
2340 | C |
---|
2341 | C ----------------------------------------------------------- |
---|
2342 | C |
---|
2343 | 5 CONTINUE |
---|
2344 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX |
---|
2345 | DO I=1,N |
---|
2346 | S1=0.0D0 |
---|
2347 | S2=0.0D0 |
---|
2348 | S3=0.0D0 |
---|
2349 | DO J=1,N |
---|
2350 | BB=FMAS(I,J) |
---|
2351 | S1=S1-BB*F1(J) |
---|
2352 | S2=S2-BB*F2(J) |
---|
2353 | S3=S3-BB*F3(J) |
---|
2354 | END DO |
---|
2355 | Z1(I)=Z1(I)+S1*FAC1 |
---|
2356 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2357 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2358 | END DO |
---|
2359 | CALL SOL (N,LDE1,E1,Z1,IP1) |
---|
2360 | CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2) |
---|
2361 | RETURN |
---|
2362 | C |
---|
2363 | C ----------------------------------------------------------- |
---|
2364 | C |
---|
2365 | 15 CONTINUE |
---|
2366 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
2367 | DO I=1,M1 |
---|
2368 | S2=-F2(I) |
---|
2369 | S3=-F3(I) |
---|
2370 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
2371 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2372 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2373 | END DO |
---|
2374 | DO I=1,NM1 |
---|
2375 | IM1=I+M1 |
---|
2376 | S1=0.0D0 |
---|
2377 | S2=0.0D0 |
---|
2378 | S3=0.0D0 |
---|
2379 | DO J=1,NM1 |
---|
2380 | JM1=J+M1 |
---|
2381 | BB=FMAS(I,J) |
---|
2382 | S1=S1-BB*F1(JM1) |
---|
2383 | S2=S2-BB*F2(JM1) |
---|
2384 | S3=S3-BB*F3(JM1) |
---|
2385 | END DO |
---|
2386 | Z1(IM1)=Z1(IM1)+S1*FAC1 |
---|
2387 | Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN |
---|
2388 | Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN |
---|
2389 | END DO |
---|
2390 | GOTO 48 |
---|
2391 | C |
---|
2392 | C ----------------------------------------------------------- |
---|
2393 | C |
---|
2394 | 6 CONTINUE |
---|
2395 | C --- B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX |
---|
2396 | C --- THIS OPTION IS NOT PROVIDED |
---|
2397 | RETURN |
---|
2398 | C |
---|
2399 | C ----------------------------------------------------------- |
---|
2400 | C |
---|
2401 | 7 CONTINUE |
---|
2402 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION |
---|
2403 | DO I=1,N |
---|
2404 | S2=-F2(I) |
---|
2405 | S3=-F3(I) |
---|
2406 | Z1(I)=Z1(I)-F1(I)*FAC1 |
---|
2407 | Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN |
---|
2408 | Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN |
---|
2409 | END DO |
---|
2410 | DO MM=N-2,1,-1 |
---|
2411 | MP=N-MM |
---|
2412 | MP1=MP-1 |
---|
2413 | I=IPHES(MP) |
---|
2414 | IF (I.EQ.MP) GOTO 746 |
---|
2415 | ZSAFE=Z1(MP) |
---|
2416 | Z1(MP)=Z1(I) |
---|
2417 | Z1(I)=ZSAFE |
---|
2418 | ZSAFE=Z2(MP) |
---|
2419 | Z2(MP)=Z2(I) |
---|
2420 | Z2(I)=ZSAFE |
---|
2421 | ZSAFE=Z3(MP) |
---|
2422 | Z3(MP)=Z3(I) |
---|
2423 | Z3(I)=ZSAFE |
---|
2424 | 746 CONTINUE |
---|
2425 | DO I=MP+1,N |
---|
2426 | E1IMP=FJAC(I,MP1) |
---|
2427 | Z1(I)=Z1(I)-E1IMP*Z1(MP) |
---|
2428 | Z2(I)=Z2(I)-E1IMP*Z2(MP) |
---|
2429 | Z3(I)=Z3(I)-E1IMP*Z3(MP) |
---|
2430 | END DO |
---|
2431 | END DO |
---|
2432 | CALL SOLH(N,LDE1,E1,1,Z1,IP1) |
---|
2433 | CALL SOLHC(N,LDE1,E2R,E2I,1,Z2,Z3,IP2) |
---|
2434 | DO MM=1,N-2 |
---|
2435 | MP=N-MM |
---|
2436 | MP1=MP-1 |
---|
2437 | DO I=MP+1,N |
---|
2438 | E1IMP=FJAC(I,MP1) |
---|
2439 | Z1(I)=Z1(I)+E1IMP*Z1(MP) |
---|
2440 | Z2(I)=Z2(I)+E1IMP*Z2(MP) |
---|
2441 | Z3(I)=Z3(I)+E1IMP*Z3(MP) |
---|
2442 | END DO |
---|
2443 | I=IPHES(MP) |
---|
2444 | IF (I.EQ.MP) GOTO 750 |
---|
2445 | ZSAFE=Z1(MP) |
---|
2446 | Z1(MP)=Z1(I) |
---|
2447 | Z1(I)=ZSAFE |
---|
2448 | ZSAFE=Z2(MP) |
---|
2449 | Z2(MP)=Z2(I) |
---|
2450 | Z2(I)=ZSAFE |
---|
2451 | ZSAFE=Z3(MP) |
---|
2452 | Z3(MP)=Z3(I) |
---|
2453 | Z3(I)=ZSAFE |
---|
2454 | 750 CONTINUE |
---|
2455 | END DO |
---|
2456 | RETURN |
---|
2457 | C |
---|
2458 | C ----------------------------------------------------------- |
---|
2459 | C |
---|
2460 | 55 CONTINUE |
---|
2461 | RETURN |
---|
2462 | END |
---|
2463 | C |
---|
2464 | C END OF SUBROUTINE SLVRAD |
---|
2465 | C |
---|
2466 | C *********************************************************** |
---|
2467 | C |
---|
2468 | SUBROUTINE ESTRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
2469 | & H,DD1,DD2,DD3,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1, |
---|
2470 | & E1,LDE1,Z1,Z2,Z3,CONT,F1,F2,IP1,IPHES,SCAL,ERR, |
---|
2471 | & FIRST,REJECT,FAC1,RPAR,IPAR) |
---|
2472 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
2473 | DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),IP1(NM1), |
---|
2474 | & SCAL(N),IPHES(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N),Y0(N),Y(N) |
---|
2475 | DIMENSION CONT(N),RPAR(1),IPAR(1) |
---|
2476 | LOGICAL FIRST,REJECT |
---|
2477 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
2478 | EXTERNAL FCN |
---|
2479 | HEE1=DD1/H |
---|
2480 | HEE2=DD2/H |
---|
2481 | HEE3=DD3/H |
---|
2482 | GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB |
---|
2483 | C |
---|
2484 | 1 CONTINUE |
---|
2485 | C ------ B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
2486 | DO I=1,N |
---|
2487 | F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2488 | CONT(I)=F2(I)+Y0(I) |
---|
2489 | END DO |
---|
2490 | CALL SOL (N,LDE1,E1,CONT,IP1) |
---|
2491 | GOTO 77 |
---|
2492 | C |
---|
2493 | 11 CONTINUE |
---|
2494 | C ------ B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
2495 | DO I=1,N |
---|
2496 | F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2497 | CONT(I)=F2(I)+Y0(I) |
---|
2498 | END DO |
---|
2499 | 48 MM=M1/M2 |
---|
2500 | DO J=1,M2 |
---|
2501 | SUM1=0.D0 |
---|
2502 | DO K=MM-1,0,-1 |
---|
2503 | SUM1=(CONT(J+K*M2)+SUM1)/FAC1 |
---|
2504 | DO I=1,NM1 |
---|
2505 | IM1=I+M1 |
---|
2506 | CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1 |
---|
2507 | END DO |
---|
2508 | END DO |
---|
2509 | END DO |
---|
2510 | CALL SOL (NM1,LDE1,E1,CONT(M1+1),IP1) |
---|
2511 | DO I=M1,1,-1 |
---|
2512 | CONT(I)=(CONT(I)+CONT(M2+I))/FAC1 |
---|
2513 | END DO |
---|
2514 | GOTO 77 |
---|
2515 | C |
---|
2516 | 2 CONTINUE |
---|
2517 | C ------ B=IDENTITY, JACOBIAN A BANDED MATRIX |
---|
2518 | DO I=1,N |
---|
2519 | F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2520 | CONT(I)=F2(I)+Y0(I) |
---|
2521 | END DO |
---|
2522 | CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1) |
---|
2523 | GOTO 77 |
---|
2524 | C |
---|
2525 | 12 CONTINUE |
---|
2526 | C ------ B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
2527 | DO I=1,N |
---|
2528 | F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2529 | CONT(I)=F2(I)+Y0(I) |
---|
2530 | END DO |
---|
2531 | 45 MM=M1/M2 |
---|
2532 | DO J=1,M2 |
---|
2533 | SUM1=0.D0 |
---|
2534 | DO K=MM-1,0,-1 |
---|
2535 | SUM1=(CONT(J+K*M2)+SUM1)/FAC1 |
---|
2536 | DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC) |
---|
2537 | IM1=I+M1 |
---|
2538 | CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1 |
---|
2539 | END DO |
---|
2540 | END DO |
---|
2541 | END DO |
---|
2542 | CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1) |
---|
2543 | DO I=M1,1,-1 |
---|
2544 | CONT(I)=(CONT(I)+CONT(M2+I))/FAC1 |
---|
2545 | END DO |
---|
2546 | GOTO 77 |
---|
2547 | C |
---|
2548 | 3 CONTINUE |
---|
2549 | C ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX |
---|
2550 | DO I=1,N |
---|
2551 | F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2552 | END DO |
---|
2553 | DO I=1,N |
---|
2554 | SUM=0.D0 |
---|
2555 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
2556 | SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J) |
---|
2557 | END DO |
---|
2558 | F2(I)=SUM |
---|
2559 | CONT(I)=SUM+Y0(I) |
---|
2560 | END DO |
---|
2561 | CALL SOL (N,LDE1,E1,CONT,IP1) |
---|
2562 | GOTO 77 |
---|
2563 | C |
---|
2564 | 13 CONTINUE |
---|
2565 | C ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
2566 | DO I=1,M1 |
---|
2567 | F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2568 | CONT(I)=F2(I)+Y0(I) |
---|
2569 | END DO |
---|
2570 | DO I=M1+1,N |
---|
2571 | F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2572 | END DO |
---|
2573 | DO I=1,NM1 |
---|
2574 | SUM=0.D0 |
---|
2575 | DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS) |
---|
2576 | SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J+M1) |
---|
2577 | END DO |
---|
2578 | IM1=I+M1 |
---|
2579 | F2(IM1)=SUM |
---|
2580 | CONT(IM1)=SUM+Y0(IM1) |
---|
2581 | END DO |
---|
2582 | GOTO 48 |
---|
2583 | C |
---|
2584 | 4 CONTINUE |
---|
2585 | C ------ B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX |
---|
2586 | DO I=1,N |
---|
2587 | F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2588 | END DO |
---|
2589 | DO I=1,N |
---|
2590 | SUM=0.D0 |
---|
2591 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
2592 | SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J) |
---|
2593 | END DO |
---|
2594 | F2(I)=SUM |
---|
2595 | CONT(I)=SUM+Y0(I) |
---|
2596 | END DO |
---|
2597 | CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1) |
---|
2598 | GOTO 77 |
---|
2599 | C |
---|
2600 | 14 CONTINUE |
---|
2601 | C ------ B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
2602 | DO I=1,M1 |
---|
2603 | F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2604 | CONT(I)=F2(I)+Y0(I) |
---|
2605 | END DO |
---|
2606 | DO I=M1+1,N |
---|
2607 | F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2608 | END DO |
---|
2609 | DO I=1,NM1 |
---|
2610 | SUM=0.D0 |
---|
2611 | DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS) |
---|
2612 | SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J+M1) |
---|
2613 | END DO |
---|
2614 | IM1=I+M1 |
---|
2615 | F2(IM1)=SUM |
---|
2616 | CONT(IM1)=SUM+Y0(IM1) |
---|
2617 | END DO |
---|
2618 | GOTO 45 |
---|
2619 | C |
---|
2620 | 5 CONTINUE |
---|
2621 | C ------ B IS A FULL MATRIX, JACOBIAN A FULL MATRIX |
---|
2622 | DO I=1,N |
---|
2623 | F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2624 | END DO |
---|
2625 | DO I=1,N |
---|
2626 | SUM=0.D0 |
---|
2627 | DO J=1,N |
---|
2628 | SUM=SUM+FMAS(I,J)*F1(J) |
---|
2629 | END DO |
---|
2630 | F2(I)=SUM |
---|
2631 | CONT(I)=SUM+Y0(I) |
---|
2632 | END DO |
---|
2633 | CALL SOL (N,LDE1,E1,CONT,IP1) |
---|
2634 | GOTO 77 |
---|
2635 | C |
---|
2636 | 15 CONTINUE |
---|
2637 | C ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
2638 | DO I=1,M1 |
---|
2639 | F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2640 | CONT(I)=F2(I)+Y0(I) |
---|
2641 | END DO |
---|
2642 | DO I=M1+1,N |
---|
2643 | F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2644 | END DO |
---|
2645 | DO I=1,NM1 |
---|
2646 | SUM=0.D0 |
---|
2647 | DO J=1,NM1 |
---|
2648 | SUM=SUM+FMAS(I,J)*F1(J+M1) |
---|
2649 | END DO |
---|
2650 | IM1=I+M1 |
---|
2651 | F2(IM1)=SUM |
---|
2652 | CONT(IM1)=SUM+Y0(IM1) |
---|
2653 | END DO |
---|
2654 | GOTO 48 |
---|
2655 | C |
---|
2656 | 6 CONTINUE |
---|
2657 | C ------ B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX |
---|
2658 | C ------ THIS OPTION IS NOT PROVIDED |
---|
2659 | RETURN |
---|
2660 | C |
---|
2661 | 7 CONTINUE |
---|
2662 | C ------ B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION |
---|
2663 | DO I=1,N |
---|
2664 | F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) |
---|
2665 | CONT(I)=F2(I)+Y0(I) |
---|
2666 | END DO |
---|
2667 | DO MM=N-2,1,-1 |
---|
2668 | MP=N-MM |
---|
2669 | I=IPHES(MP) |
---|
2670 | IF (I.EQ.MP) GOTO 310 |
---|
2671 | ZSAFE=CONT(MP) |
---|
2672 | CONT(MP)=CONT(I) |
---|
2673 | CONT(I)=ZSAFE |
---|
2674 | 310 CONTINUE |
---|
2675 | DO I=MP+1,N |
---|
2676 | CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP) |
---|
2677 | END DO |
---|
2678 | END DO |
---|
2679 | CALL SOLH(N,LDE1,E1,1,CONT,IP1) |
---|
2680 | DO MM=1,N-2 |
---|
2681 | MP=N-MM |
---|
2682 | DO I=MP+1,N |
---|
2683 | CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP) |
---|
2684 | END DO |
---|
2685 | I=IPHES(MP) |
---|
2686 | IF (I.EQ.MP) GOTO 440 |
---|
2687 | ZSAFE=CONT(MP) |
---|
2688 | CONT(MP)=CONT(I) |
---|
2689 | CONT(I)=ZSAFE |
---|
2690 | 440 CONTINUE |
---|
2691 | END DO |
---|
2692 | C |
---|
2693 | C -------------------------------------- |
---|
2694 | C |
---|
2695 | 77 CONTINUE |
---|
2696 | ERR=0.D0 |
---|
2697 | DO I=1,N |
---|
2698 | ERR=ERR+(CONT(I)/SCAL(I))**2 |
---|
2699 | END DO |
---|
2700 | ERR=MAX(SQRT(ERR/N),1.D-10) |
---|
2701 | C |
---|
2702 | IF (ERR.LT.1.D0) RETURN |
---|
2703 | IF (FIRST.OR.REJECT) THEN |
---|
2704 | DO I=1,N |
---|
2705 | CONT(I)=Y(I)+CONT(I) |
---|
2706 | END DO |
---|
2707 | CALL FCN(N,X,CONT,F1,RPAR,IPAR) |
---|
2708 | NFCN=NFCN+1 |
---|
2709 | DO I=1,N |
---|
2710 | CONT(I)=F1(I)+F2(I) |
---|
2711 | END DO |
---|
2712 | GOTO (31,32,31,32,31,32,33,55,55,55,41,42,41,42,41), IJOB |
---|
2713 | C ------ FULL MATRIX OPTION |
---|
2714 | 31 CONTINUE |
---|
2715 | CALL SOL(N,LDE1,E1,CONT,IP1) |
---|
2716 | GOTO 88 |
---|
2717 | C ------ FULL MATRIX OPTION, SECOND ORDER |
---|
2718 | 41 CONTINUE |
---|
2719 | DO J=1,M2 |
---|
2720 | SUM1=0.D0 |
---|
2721 | DO K=MM-1,0,-1 |
---|
2722 | SUM1=(CONT(J+K*M2)+SUM1)/FAC1 |
---|
2723 | DO I=1,NM1 |
---|
2724 | IM1=I+M1 |
---|
2725 | CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1 |
---|
2726 | END DO |
---|
2727 | END DO |
---|
2728 | END DO |
---|
2729 | CALL SOL(NM1,LDE1,E1,CONT(M1+1),IP1) |
---|
2730 | DO I=M1,1,-1 |
---|
2731 | CONT(I)=(CONT(I)+CONT(M2+I))/FAC1 |
---|
2732 | END DO |
---|
2733 | GOTO 88 |
---|
2734 | C ------ BANDED MATRIX OPTION |
---|
2735 | 32 CONTINUE |
---|
2736 | CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1) |
---|
2737 | GOTO 88 |
---|
2738 | C ------ BANDED MATRIX OPTION, SECOND ORDER |
---|
2739 | 42 CONTINUE |
---|
2740 | DO J=1,M2 |
---|
2741 | SUM1=0.D0 |
---|
2742 | DO K=MM-1,0,-1 |
---|
2743 | SUM1=(CONT(J+K*M2)+SUM1)/FAC1 |
---|
2744 | DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC) |
---|
2745 | IM1=I+M1 |
---|
2746 | CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1 |
---|
2747 | END DO |
---|
2748 | END DO |
---|
2749 | END DO |
---|
2750 | CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1) |
---|
2751 | DO I=M1,1,-1 |
---|
2752 | CONT(I)=(CONT(I)+CONT(M2+I))/FAC1 |
---|
2753 | END DO |
---|
2754 | GOTO 88 |
---|
2755 | C ------ HESSENBERG MATRIX OPTION |
---|
2756 | 33 CONTINUE |
---|
2757 | DO MM=N-2,1,-1 |
---|
2758 | MP=N-MM |
---|
2759 | I=IPHES(MP) |
---|
2760 | IF (I.EQ.MP) GOTO 510 |
---|
2761 | ZSAFE=CONT(MP) |
---|
2762 | CONT(MP)=CONT(I) |
---|
2763 | CONT(I)=ZSAFE |
---|
2764 | 510 CONTINUE |
---|
2765 | DO I=MP+1,N |
---|
2766 | CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP) |
---|
2767 | END DO |
---|
2768 | END DO |
---|
2769 | CALL SOLH(N,LDE1,E1,1,CONT,IP1) |
---|
2770 | DO MM=1,N-2 |
---|
2771 | MP=N-MM |
---|
2772 | DO I=MP+1,N |
---|
2773 | CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP) |
---|
2774 | END DO |
---|
2775 | I=IPHES(MP) |
---|
2776 | IF (I.EQ.MP) GOTO 640 |
---|
2777 | ZSAFE=CONT(MP) |
---|
2778 | CONT(MP)=CONT(I) |
---|
2779 | CONT(I)=ZSAFE |
---|
2780 | 640 CONTINUE |
---|
2781 | END DO |
---|
2782 | C ----------------------------------- |
---|
2783 | 88 CONTINUE |
---|
2784 | ERR=0.D0 |
---|
2785 | DO I=1,N |
---|
2786 | ERR=ERR+(CONT(I)/SCAL(I))**2 |
---|
2787 | END DO |
---|
2788 | ERR=MAX(SQRT(ERR/N),1.D-10) |
---|
2789 | END IF |
---|
2790 | RETURN |
---|
2791 | C ----------------------------------------------------------- |
---|
2792 | 55 CONTINUE |
---|
2793 | RETURN |
---|
2794 | END |
---|
2795 | C |
---|
2796 | C END OF SUBROUTINE ESTRAD |
---|
2797 | C |
---|
2798 | C *********************************************************** |
---|
2799 | C |
---|
2800 | SUBROUTINE ESTRAV(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
2801 | & H,DD,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,NS,NNS, |
---|
2802 | & E1,LDE1,ZZ,CONT,FF,IP1,IPHES,SCAL,ERR, |
---|
2803 | & FIRST,REJECT,FAC1,RPAR,IPAR) |
---|
2804 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
2805 | DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),IP1(NM1), |
---|
2806 | & SCAL(N),IPHES(N),ZZ(NNS),FF(NNS),Y0(N),Y(N) |
---|
2807 | DIMENSION DD(NS),CONT(N),RPAR(1),IPAR(1) |
---|
2808 | LOGICAL FIRST,REJECT |
---|
2809 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
2810 | EXTERNAL FCN |
---|
2811 | |
---|
2812 | GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB |
---|
2813 | C |
---|
2814 | 1 CONTINUE |
---|
2815 | C ------ B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
2816 | DO I=1,N |
---|
2817 | SUM=0.D0 |
---|
2818 | DO K=1,NS |
---|
2819 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2820 | END DO |
---|
2821 | FF(I+N)=SUM/H |
---|
2822 | CONT(I)=FF(I+N)+Y0(I) |
---|
2823 | END DO |
---|
2824 | CALL SOL (N,LDE1,E1,CONT,IP1) |
---|
2825 | GOTO 77 |
---|
2826 | C |
---|
2827 | 11 CONTINUE |
---|
2828 | C ------ B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
2829 | DO I=1,N |
---|
2830 | SUM=0.D0 |
---|
2831 | DO K=1,NS |
---|
2832 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2833 | END DO |
---|
2834 | FF(I+N)=SUM/H |
---|
2835 | CONT(I)=FF(I+N)+Y0(I) |
---|
2836 | END DO |
---|
2837 | 48 MM=M1/M2 |
---|
2838 | DO J=1,M2 |
---|
2839 | SUM1=0.D0 |
---|
2840 | DO K=MM-1,0,-1 |
---|
2841 | SUM1=(CONT(J+K*M2)+SUM1)/FAC1 |
---|
2842 | DO I=1,NM1 |
---|
2843 | IM1=I+M1 |
---|
2844 | CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1 |
---|
2845 | END DO |
---|
2846 | END DO |
---|
2847 | END DO |
---|
2848 | CALL SOL (NM1,LDE1,E1,CONT(M1+1),IP1) |
---|
2849 | DO I=M1,1,-1 |
---|
2850 | CONT(I)=(CONT(I)+CONT(M2+I))/FAC1 |
---|
2851 | END DO |
---|
2852 | GOTO 77 |
---|
2853 | C |
---|
2854 | 2 CONTINUE |
---|
2855 | C ------ B=IDENTITY, JACOBIAN A BANDED MATRIX |
---|
2856 | DO I=1,N |
---|
2857 | SUM=0.D0 |
---|
2858 | DO K=1,NS |
---|
2859 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2860 | END DO |
---|
2861 | FF(I+N)=SUM/H |
---|
2862 | CONT(I)=FF(I+N)+Y0(I) |
---|
2863 | END DO |
---|
2864 | CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1) |
---|
2865 | GOTO 77 |
---|
2866 | C |
---|
2867 | 12 CONTINUE |
---|
2868 | C ------ B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
2869 | DO I=1,N |
---|
2870 | SUM=0.D0 |
---|
2871 | DO K=1,NS |
---|
2872 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2873 | END DO |
---|
2874 | FF(I+N)=SUM/H |
---|
2875 | CONT(I)=FF(I+N)+Y0(I) |
---|
2876 | END DO |
---|
2877 | 45 MM=M1/M2 |
---|
2878 | DO J=1,M2 |
---|
2879 | SUM1=0.D0 |
---|
2880 | DO K=MM-1,0,-1 |
---|
2881 | SUM1=(CONT(J+K*M2)+SUM1)/FAC1 |
---|
2882 | DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC) |
---|
2883 | IM1=I+M1 |
---|
2884 | CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1 |
---|
2885 | END DO |
---|
2886 | END DO |
---|
2887 | END DO |
---|
2888 | CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1) |
---|
2889 | DO I=M1,1,-1 |
---|
2890 | CONT(I)=(CONT(I)+CONT(M2+I))/FAC1 |
---|
2891 | END DO |
---|
2892 | GOTO 77 |
---|
2893 | C |
---|
2894 | 3 CONTINUE |
---|
2895 | C ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX |
---|
2896 | DO I=1,N |
---|
2897 | SUM=0.D0 |
---|
2898 | DO K=1,NS |
---|
2899 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2900 | END DO |
---|
2901 | FF(I)=SUM/H |
---|
2902 | END DO |
---|
2903 | DO I=1,N |
---|
2904 | SUM=0.D0 |
---|
2905 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
2906 | SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J) |
---|
2907 | END DO |
---|
2908 | FF(I+N)=SUM |
---|
2909 | CONT(I)=SUM+Y0(I) |
---|
2910 | END DO |
---|
2911 | CALL SOL (N,LDE1,E1,CONT,IP1) |
---|
2912 | GOTO 77 |
---|
2913 | C |
---|
2914 | 13 CONTINUE |
---|
2915 | C ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
2916 | DO I=1,M1 |
---|
2917 | SUM=0.D0 |
---|
2918 | DO K=1,NS |
---|
2919 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2920 | END DO |
---|
2921 | FF(I+N)=SUM/H |
---|
2922 | CONT(I)=FF(I+N)+Y0(I) |
---|
2923 | END DO |
---|
2924 | DO I=M1+1,N |
---|
2925 | SUM=0.D0 |
---|
2926 | DO K=1,NS |
---|
2927 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2928 | END DO |
---|
2929 | FF(I)=SUM/H |
---|
2930 | END DO |
---|
2931 | DO I=1,NM1 |
---|
2932 | SUM=0.D0 |
---|
2933 | DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS) |
---|
2934 | SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J+M1) |
---|
2935 | END DO |
---|
2936 | IM1=I+M1 |
---|
2937 | FF(IM1+N)=SUM |
---|
2938 | CONT(IM1)=SUM+Y0(IM1) |
---|
2939 | END DO |
---|
2940 | GOTO 48 |
---|
2941 | C |
---|
2942 | 4 CONTINUE |
---|
2943 | C ------ B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX |
---|
2944 | DO I=1,N |
---|
2945 | SUM=0.D0 |
---|
2946 | DO K=1,NS |
---|
2947 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2948 | END DO |
---|
2949 | FF(I)=SUM/H |
---|
2950 | END DO |
---|
2951 | DO I=1,N |
---|
2952 | SUM=0.D0 |
---|
2953 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
2954 | SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J) |
---|
2955 | END DO |
---|
2956 | FF(I+N)=SUM |
---|
2957 | CONT(I)=SUM+Y0(I) |
---|
2958 | END DO |
---|
2959 | CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1) |
---|
2960 | GOTO 77 |
---|
2961 | C |
---|
2962 | 14 CONTINUE |
---|
2963 | C ------ B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
2964 | DO I=1,M1 |
---|
2965 | SUM=0.D0 |
---|
2966 | DO K=1,NS |
---|
2967 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2968 | END DO |
---|
2969 | FF(I+N)=SUM/H |
---|
2970 | CONT(I)=FF(I+N)+Y0(I) |
---|
2971 | END DO |
---|
2972 | DO I=M1+1,N |
---|
2973 | SUM=0.D0 |
---|
2974 | DO K=1,NS |
---|
2975 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2976 | END DO |
---|
2977 | FF(I)=SUM/H |
---|
2978 | END DO |
---|
2979 | DO I=1,NM1 |
---|
2980 | SUM=0.D0 |
---|
2981 | DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS) |
---|
2982 | SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J+M1) |
---|
2983 | END DO |
---|
2984 | IM1=I+M1 |
---|
2985 | FF(IM1+N)=SUM |
---|
2986 | CONT(IM1)=SUM+Y0(IM1) |
---|
2987 | END DO |
---|
2988 | GOTO 45 |
---|
2989 | C |
---|
2990 | 5 CONTINUE |
---|
2991 | C ------ B IS A FULL MATRIX, JACOBIAN A FULL MATRIX |
---|
2992 | DO I=1,N |
---|
2993 | SUM=0.D0 |
---|
2994 | DO K=1,NS |
---|
2995 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
2996 | END DO |
---|
2997 | FF(I)=SUM/H |
---|
2998 | END DO |
---|
2999 | DO I=1,N |
---|
3000 | SUM=0.D0 |
---|
3001 | DO J=1,N |
---|
3002 | SUM=SUM+FMAS(I,J)*FF(J) |
---|
3003 | END DO |
---|
3004 | FF(I+N)=SUM |
---|
3005 | CONT(I)=SUM+Y0(I) |
---|
3006 | END DO |
---|
3007 | CALL SOL (N,LDE1,E1,CONT,IP1) |
---|
3008 | GOTO 77 |
---|
3009 | C |
---|
3010 | 15 CONTINUE |
---|
3011 | C ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
3012 | DO I=1,M1 |
---|
3013 | SUM=0.D0 |
---|
3014 | DO K=1,NS |
---|
3015 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
3016 | END DO |
---|
3017 | FF(I+N)=SUM/H |
---|
3018 | CONT(I)=FF(I+N)+Y0(I) |
---|
3019 | END DO |
---|
3020 | DO I=M1+1,N |
---|
3021 | SUM=0.D0 |
---|
3022 | DO K=1,NS |
---|
3023 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
3024 | END DO |
---|
3025 | FF(I)=SUM/H |
---|
3026 | END DO |
---|
3027 | DO I=1,NM1 |
---|
3028 | SUM=0.D0 |
---|
3029 | DO J=1,NM1 |
---|
3030 | SUM=SUM+FMAS(I,J)*FF(J+M1) |
---|
3031 | END DO |
---|
3032 | IM1=I+M1 |
---|
3033 | FF(IM1+N)=SUM |
---|
3034 | CONT(IM1)=SUM+Y0(IM1) |
---|
3035 | END DO |
---|
3036 | GOTO 48 |
---|
3037 | C |
---|
3038 | 6 CONTINUE |
---|
3039 | C ------ B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX |
---|
3040 | C ------ THIS OPTION IS NOT PROVIDED |
---|
3041 | RETURN |
---|
3042 | C |
---|
3043 | 7 CONTINUE |
---|
3044 | C ------ B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION |
---|
3045 | DO I=1,N |
---|
3046 | SUM=0.D0 |
---|
3047 | DO K=1,NS |
---|
3048 | SUM=SUM+DD(K)*ZZ(I+(K-1)*N) |
---|
3049 | END DO |
---|
3050 | FF(I+N)=SUM/H |
---|
3051 | CONT(I)=FF(I+N)+Y0(I) |
---|
3052 | END DO |
---|
3053 | DO MM=N-2,1,-1 |
---|
3054 | MP=N-MM |
---|
3055 | I=IPHES(MP) |
---|
3056 | IF (I.EQ.MP) GOTO 310 |
---|
3057 | ZSAFE=CONT(MP) |
---|
3058 | CONT(MP)=CONT(I) |
---|
3059 | CONT(I)=ZSAFE |
---|
3060 | 310 CONTINUE |
---|
3061 | DO I=MP+1,N |
---|
3062 | CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP) |
---|
3063 | END DO |
---|
3064 | END DO |
---|
3065 | CALL SOLH(N,LDE1,E1,1,CONT,IP1) |
---|
3066 | DO MM=1,N-2 |
---|
3067 | MP=N-MM |
---|
3068 | DO I=MP+1,N |
---|
3069 | CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP) |
---|
3070 | END DO |
---|
3071 | I=IPHES(MP) |
---|
3072 | IF (I.EQ.MP) GOTO 440 |
---|
3073 | ZSAFE=CONT(MP) |
---|
3074 | CONT(MP)=CONT(I) |
---|
3075 | CONT(I)=ZSAFE |
---|
3076 | 440 CONTINUE |
---|
3077 | END DO |
---|
3078 | C |
---|
3079 | C -------------------------------------- |
---|
3080 | C |
---|
3081 | 77 CONTINUE |
---|
3082 | ERR=0.D0 |
---|
3083 | DO I=1,N |
---|
3084 | ERR=ERR+(CONT(I)/SCAL(I))**2 |
---|
3085 | END DO |
---|
3086 | ERR=MAX(SQRT(ERR/N),1.D-10) |
---|
3087 | C |
---|
3088 | IF (ERR.LT.1.D0) RETURN |
---|
3089 | IF (FIRST.OR.REJECT) THEN |
---|
3090 | DO I=1,N |
---|
3091 | CONT(I)=Y(I)+CONT(I) |
---|
3092 | END DO |
---|
3093 | CALL FCN(N,X,CONT,FF,RPAR,IPAR) |
---|
3094 | NFCN=NFCN+1 |
---|
3095 | DO I=1,N |
---|
3096 | CONT(I)=FF(I)+FF(I+N) |
---|
3097 | END DO |
---|
3098 | GOTO (31,32,31,32,31,32,33,55,55,55,41,42,41,42,41), IJOB |
---|
3099 | C ------ FULL MATRIX OPTION |
---|
3100 | 31 CONTINUE |
---|
3101 | CALL SOL (N,LDE1,E1,CONT,IP1) |
---|
3102 | GOTO 88 |
---|
3103 | C ------ FULL MATRIX OPTION, SECOND ORDER |
---|
3104 | 41 CONTINUE |
---|
3105 | DO J=1,M2 |
---|
3106 | SUM1=0.D0 |
---|
3107 | DO K=MM-1,0,-1 |
---|
3108 | SUM1=(CONT(J+K*M2)+SUM1)/FAC1 |
---|
3109 | DO I=1,NM1 |
---|
3110 | IM1=I+M1 |
---|
3111 | CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1 |
---|
3112 | END DO |
---|
3113 | END DO |
---|
3114 | END DO |
---|
3115 | CALL SOL (NM1,LDE1,E1,CONT(M1+1),IP1) |
---|
3116 | DO I=M1,1,-1 |
---|
3117 | CONT(I)=(CONT(I)+CONT(M2+I))/FAC1 |
---|
3118 | END DO |
---|
3119 | GOTO 88 |
---|
3120 | C ------ BANDED MATRIX OPTION |
---|
3121 | 32 CONTINUE |
---|
3122 | CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1) |
---|
3123 | GOTO 88 |
---|
3124 | C ------ BANDED MATRIX OPTION, SECOND ORDER |
---|
3125 | 42 CONTINUE |
---|
3126 | DO J=1,M2 |
---|
3127 | SUM1=0.D0 |
---|
3128 | DO K=MM-1,0,-1 |
---|
3129 | SUM1=(CONT(J+K*M2)+SUM1)/FAC1 |
---|
3130 | DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC) |
---|
3131 | IM1=I+M1 |
---|
3132 | CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1 |
---|
3133 | END DO |
---|
3134 | END DO |
---|
3135 | END DO |
---|
3136 | CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1) |
---|
3137 | DO I=M1,1,-1 |
---|
3138 | CONT(I)=(CONT(I)+CONT(M2+I))/FAC1 |
---|
3139 | END DO |
---|
3140 | GOTO 88 |
---|
3141 | C ------ HESSENBERG MATRIX OPTION |
---|
3142 | 33 CONTINUE |
---|
3143 | DO MM=N-2,1,-1 |
---|
3144 | MP=N-MM |
---|
3145 | I=IPHES(MP) |
---|
3146 | IF (I.EQ.MP) GOTO 510 |
---|
3147 | ZSAFE=CONT(MP) |
---|
3148 | CONT(MP)=CONT(I) |
---|
3149 | CONT(I)=ZSAFE |
---|
3150 | 510 CONTINUE |
---|
3151 | DO I=MP+1,N |
---|
3152 | CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP) |
---|
3153 | END DO |
---|
3154 | END DO |
---|
3155 | CALL SOLH(N,LDE1,E1,1,CONT,IP1) |
---|
3156 | DO MM=1,N-2 |
---|
3157 | MP=N-MM |
---|
3158 | DO I=MP+1,N |
---|
3159 | CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP) |
---|
3160 | END DO |
---|
3161 | I=IPHES(MP) |
---|
3162 | IF (I.EQ.MP) GOTO 640 |
---|
3163 | ZSAFE=CONT(MP) |
---|
3164 | CONT(MP)=CONT(I) |
---|
3165 | CONT(I)=ZSAFE |
---|
3166 | 640 CONTINUE |
---|
3167 | END DO |
---|
3168 | C ----------------------------------- |
---|
3169 | 88 CONTINUE |
---|
3170 | ERR=0.D0 |
---|
3171 | DO I=1,N |
---|
3172 | ERR=ERR+(CONT(I)/SCAL(I))**2 |
---|
3173 | END DO |
---|
3174 | ERR=MAX(SQRT(ERR/N),1.D-10) |
---|
3175 | END IF |
---|
3176 | RETURN |
---|
3177 | C |
---|
3178 | C ----------------------------------------------------------- |
---|
3179 | C |
---|
3180 | 55 CONTINUE |
---|
3181 | RETURN |
---|
3182 | END |
---|
3183 | C |
---|
3184 | C END OF SUBROUTINE ESTRAV |
---|
3185 | C |
---|
3186 | C *********************************************************** |
---|
3187 | C |
---|
3188 | SUBROUTINE SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
3189 | & M1,M2,NM1,FAC1,E,LDE,IP,DY,AK,FX,YNEW,HD,IJOB,STAGE1) |
---|
3190 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
3191 | DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E(LDE,NM1), |
---|
3192 | & IP(NM1),DY(N),AK(N),FX(N),YNEW(N) |
---|
3193 | LOGICAL STAGE1 |
---|
3194 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
3195 | C |
---|
3196 | IF (HD.EQ.0.D0) THEN |
---|
3197 | DO I=1,N |
---|
3198 | AK(I)=DY(I) |
---|
3199 | END DO |
---|
3200 | ELSE |
---|
3201 | DO I=1,N |
---|
3202 | AK(I)=DY(I)+HD*FX(I) |
---|
3203 | END DO |
---|
3204 | END IF |
---|
3205 | C |
---|
3206 | GOTO (1,2,3,4,5,6,55,55,55,55,11,12,13,13,15), IJOB |
---|
3207 | C |
---|
3208 | C ----------------------------------------------------------- |
---|
3209 | C |
---|
3210 | 1 CONTINUE |
---|
3211 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
3212 | IF (STAGE1) THEN |
---|
3213 | DO I=1,N |
---|
3214 | AK(I)=AK(I)+YNEW(I) |
---|
3215 | END DO |
---|
3216 | END IF |
---|
3217 | CALL SOL (N,LDE,E,AK,IP) |
---|
3218 | RETURN |
---|
3219 | C |
---|
3220 | C ----------------------------------------------------------- |
---|
3221 | C |
---|
3222 | 11 CONTINUE |
---|
3223 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
3224 | IF (STAGE1) THEN |
---|
3225 | DO I=1,N |
---|
3226 | AK(I)=AK(I)+YNEW(I) |
---|
3227 | END DO |
---|
3228 | END IF |
---|
3229 | 48 MM=M1/M2 |
---|
3230 | DO J=1,M2 |
---|
3231 | SUM=0.D0 |
---|
3232 | DO K=MM-1,0,-1 |
---|
3233 | JKM=J+K*M2 |
---|
3234 | SUM=(AK(JKM)+SUM)/FAC1 |
---|
3235 | DO I=1,NM1 |
---|
3236 | IM1=I+M1 |
---|
3237 | AK(IM1)=AK(IM1)+FJAC(I,JKM)*SUM |
---|
3238 | END DO |
---|
3239 | END DO |
---|
3240 | END DO |
---|
3241 | CALL SOL (NM1,LDE,E,AK(M1+1),IP) |
---|
3242 | DO I=M1,1,-1 |
---|
3243 | AK(I)=(AK(I)+AK(M2+I))/FAC1 |
---|
3244 | END DO |
---|
3245 | RETURN |
---|
3246 | C |
---|
3247 | C ----------------------------------------------------------- |
---|
3248 | C |
---|
3249 | 2 CONTINUE |
---|
3250 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX |
---|
3251 | IF (STAGE1) THEN |
---|
3252 | DO I=1,N |
---|
3253 | AK(I)=AK(I)+YNEW(I) |
---|
3254 | END DO |
---|
3255 | END IF |
---|
3256 | CALL SOLB (N,LDE,E,MLE,MUE,AK,IP) |
---|
3257 | RETURN |
---|
3258 | C |
---|
3259 | C ----------------------------------------------------------- |
---|
3260 | C |
---|
3261 | 12 CONTINUE |
---|
3262 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
3263 | IF (STAGE1) THEN |
---|
3264 | DO I=1,N |
---|
3265 | AK(I)=AK(I)+YNEW(I) |
---|
3266 | END DO |
---|
3267 | END IF |
---|
3268 | 45 MM=M1/M2 |
---|
3269 | DO J=1,M2 |
---|
3270 | SUM=0.D0 |
---|
3271 | DO K=MM-1,0,-1 |
---|
3272 | JKM=J+K*M2 |
---|
3273 | SUM=(AK(JKM)+SUM)/FAC1 |
---|
3274 | DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC) |
---|
3275 | IM1=I+M1 |
---|
3276 | AK(IM1)=AK(IM1)+FJAC(I+MUJAC+1-J,JKM)*SUM |
---|
3277 | END DO |
---|
3278 | END DO |
---|
3279 | END DO |
---|
3280 | CALL SOLB (NM1,LDE,E,MLE,MUE,AK(M1+1),IP) |
---|
3281 | DO I=M1,1,-1 |
---|
3282 | AK(I)=(AK(I)+AK(M2+I))/FAC1 |
---|
3283 | END DO |
---|
3284 | RETURN |
---|
3285 | C |
---|
3286 | C ----------------------------------------------------------- |
---|
3287 | C |
---|
3288 | 3 CONTINUE |
---|
3289 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX |
---|
3290 | IF (STAGE1) THEN |
---|
3291 | DO I=1,N |
---|
3292 | SUM=0.D0 |
---|
3293 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
3294 | SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J) |
---|
3295 | END DO |
---|
3296 | AK(I)=AK(I)+SUM |
---|
3297 | END DO |
---|
3298 | END IF |
---|
3299 | CALL SOL (N,LDE,E,AK,IP) |
---|
3300 | RETURN |
---|
3301 | C |
---|
3302 | C ----------------------------------------------------------- |
---|
3303 | C |
---|
3304 | 13 CONTINUE |
---|
3305 | C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
3306 | IF (STAGE1) THEN |
---|
3307 | DO I=1,M1 |
---|
3308 | AK(I)=AK(I)+YNEW(I) |
---|
3309 | END DO |
---|
3310 | DO I=1,NM1 |
---|
3311 | SUM=0.D0 |
---|
3312 | DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS) |
---|
3313 | SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J+M1) |
---|
3314 | END DO |
---|
3315 | IM1=I+M1 |
---|
3316 | AK(IM1)=AK(IM1)+SUM |
---|
3317 | END DO |
---|
3318 | END IF |
---|
3319 | IF (IJOB.EQ.14) GOTO 45 |
---|
3320 | GOTO 48 |
---|
3321 | C |
---|
3322 | C ----------------------------------------------------------- |
---|
3323 | C |
---|
3324 | 4 CONTINUE |
---|
3325 | C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX |
---|
3326 | IF (STAGE1) THEN |
---|
3327 | DO I=1,N |
---|
3328 | SUM=0.D0 |
---|
3329 | DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) |
---|
3330 | SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J) |
---|
3331 | END DO |
---|
3332 | AK(I)=AK(I)+SUM |
---|
3333 | END DO |
---|
3334 | END IF |
---|
3335 | CALL SOLB (N,LDE,E,MLE,MUE,AK,IP) |
---|
3336 | RETURN |
---|
3337 | C |
---|
3338 | C ----------------------------------------------------------- |
---|
3339 | C |
---|
3340 | 5 CONTINUE |
---|
3341 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX |
---|
3342 | IF (STAGE1) THEN |
---|
3343 | DO I=1,N |
---|
3344 | SUM=0.D0 |
---|
3345 | DO J=1,N |
---|
3346 | SUM=SUM+FMAS(I,J)*YNEW(J) |
---|
3347 | END DO |
---|
3348 | AK(I)=AK(I)+SUM |
---|
3349 | END DO |
---|
3350 | END IF |
---|
3351 | CALL SOL (N,LDE,E,AK,IP) |
---|
3352 | RETURN |
---|
3353 | C |
---|
3354 | C ----------------------------------------------------------- |
---|
3355 | C |
---|
3356 | 15 CONTINUE |
---|
3357 | C --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
3358 | IF (STAGE1) THEN |
---|
3359 | DO I=1,M1 |
---|
3360 | AK(I)=AK(I)+YNEW(I) |
---|
3361 | END DO |
---|
3362 | DO I=1,NM1 |
---|
3363 | SUM=0.D0 |
---|
3364 | DO J=1,NM1 |
---|
3365 | SUM=SUM+FMAS(I,J)*YNEW(J+M1) |
---|
3366 | END DO |
---|
3367 | IM1=I+M1 |
---|
3368 | AK(IM1)=AK(IM1)+SUM |
---|
3369 | END DO |
---|
3370 | END IF |
---|
3371 | GOTO 48 |
---|
3372 | C |
---|
3373 | C ----------------------------------------------------------- |
---|
3374 | C |
---|
3375 | 6 CONTINUE |
---|
3376 | C --- B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX |
---|
3377 | C --- THIS OPTION IS NOT PROVIDED |
---|
3378 | IF (STAGE1) THEN |
---|
3379 | DO 624 I=1,N |
---|
3380 | SUM=0.D0 |
---|
3381 | DO 623 J=1,N |
---|
3382 | 623 SUM=SUM+FMAS(I,J)*YNEW(J) |
---|
3383 | 624 AK(I)=AK(I)+SUM |
---|
3384 | CALL SOLB (N,LDE,E,MLE,MUE,AK,IP) |
---|
3385 | END IF |
---|
3386 | RETURN |
---|
3387 | C |
---|
3388 | C ----------------------------------------------------------- |
---|
3389 | C |
---|
3390 | 55 CONTINUE |
---|
3391 | RETURN |
---|
3392 | END |
---|
3393 | C |
---|
3394 | C END OF SUBROUTINE SLVROD |
---|
3395 | C |
---|
3396 | C |
---|
3397 | C *********************************************************** |
---|
3398 | C |
---|
3399 | SUBROUTINE SLVSEU(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
3400 | & M1,M2,NM1,FAC1,E,LDE,IP,IPHES,DEL,IJOB) |
---|
3401 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
3402 | DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E(LDE,NM1),DEL(N) |
---|
3403 | DIMENSION IP(NM1),IPHES(N) |
---|
3404 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
3405 | C |
---|
3406 | GOTO (1,2,1,2,1,55,7,55,55,55,11,12,11,12,11), IJOB |
---|
3407 | C |
---|
3408 | C ----------------------------------------------------------- |
---|
3409 | C |
---|
3410 | 1 CONTINUE |
---|
3411 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
3412 | CALL SOL (N,LDE,E,DEL,IP) |
---|
3413 | RETURN |
---|
3414 | C |
---|
3415 | C ----------------------------------------------------------- |
---|
3416 | C |
---|
3417 | 11 CONTINUE |
---|
3418 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER |
---|
3419 | MM=M1/M2 |
---|
3420 | DO J=1,M2 |
---|
3421 | SUM=0.D0 |
---|
3422 | DO K=MM-1,0,-1 |
---|
3423 | JKM=J+K*M2 |
---|
3424 | SUM=(DEL(JKM)+SUM)/FAC1 |
---|
3425 | DO I=1,NM1 |
---|
3426 | IM1=I+M1 |
---|
3427 | DEL(IM1)=DEL(IM1)+FJAC(I,JKM)*SUM |
---|
3428 | END DO |
---|
3429 | END DO |
---|
3430 | END DO |
---|
3431 | CALL SOL (NM1,LDE,E,DEL(M1+1),IP) |
---|
3432 | DO I=M1,1,-1 |
---|
3433 | DEL(I)=(DEL(I)+DEL(M2+I))/FAC1 |
---|
3434 | END DO |
---|
3435 | RETURN |
---|
3436 | C |
---|
3437 | C ----------------------------------------------------------- |
---|
3438 | C |
---|
3439 | 2 CONTINUE |
---|
3440 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX |
---|
3441 | CALL SOLB (N,LDE,E,MLE,MUE,DEL,IP) |
---|
3442 | RETURN |
---|
3443 | C |
---|
3444 | C ----------------------------------------------------------- |
---|
3445 | C |
---|
3446 | 12 CONTINUE |
---|
3447 | C --- B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER |
---|
3448 | MM=M1/M2 |
---|
3449 | DO J=1,M2 |
---|
3450 | SUM=0.D0 |
---|
3451 | DO K=MM-1,0,-1 |
---|
3452 | JKM=J+K*M2 |
---|
3453 | SUM=(DEL(JKM)+SUM)/FAC1 |
---|
3454 | DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC) |
---|
3455 | IM1=I+M1 |
---|
3456 | DEL(IM1)=DEL(IM1)+FJAC(I+MUJAC+1-J,JKM)*SUM |
---|
3457 | END DO |
---|
3458 | END DO |
---|
3459 | END DO |
---|
3460 | CALL SOLB (NM1,LDE,E,MLE,MUE,DEL(M1+1),IP) |
---|
3461 | DO I=M1,1,-1 |
---|
3462 | DEL(I)=(DEL(I)+DEL(M2+I))/FAC1 |
---|
3463 | END DO |
---|
3464 | RETURN |
---|
3465 | C |
---|
3466 | C ----------------------------------------------------------- |
---|
3467 | C |
---|
3468 | 7 CONTINUE |
---|
3469 | C --- HESSENBERG OPTION |
---|
3470 | DO MMM=N-2,1,-1 |
---|
3471 | MP=N-MMM |
---|
3472 | MP1=MP-1 |
---|
3473 | I=IPHES(MP) |
---|
3474 | IF (I.EQ.MP) GOTO 110 |
---|
3475 | ZSAFE=DEL(MP) |
---|
3476 | DEL(MP)=DEL(I) |
---|
3477 | DEL(I)=ZSAFE |
---|
3478 | 110 CONTINUE |
---|
3479 | DO I=MP+1,N |
---|
3480 | DEL(I)=DEL(I)-FJAC(I,MP1)*DEL(MP) |
---|
3481 | END DO |
---|
3482 | END DO |
---|
3483 | CALL SOLH(N,LDE,E,1,DEL,IP) |
---|
3484 | DO MMM=1,N-2 |
---|
3485 | MP=N-MMM |
---|
3486 | MP1=MP-1 |
---|
3487 | DO I=MP+1,N |
---|
3488 | DEL(I)=DEL(I)+FJAC(I,MP1)*DEL(MP) |
---|
3489 | END DO |
---|
3490 | I=IPHES(MP) |
---|
3491 | IF (I.EQ.MP) GOTO 240 |
---|
3492 | ZSAFE=DEL(MP) |
---|
3493 | DEL(MP)=DEL(I) |
---|
3494 | DEL(I)=ZSAFE |
---|
3495 | 240 CONTINUE |
---|
3496 | END DO |
---|
3497 | RETURN |
---|
3498 | C |
---|
3499 | C ----------------------------------------------------------- |
---|
3500 | C |
---|
3501 | 55 CONTINUE |
---|
3502 | RETURN |
---|
3503 | END |
---|
3504 | C |
---|
3505 | C END OF SUBROUTINE SLVSEU |
---|
3506 | C |
---|
3507 | SUBROUTINE DEC (N, NDIM, A, IP, IER) |
---|
3508 | C VERSION REAL KPP_REAL |
---|
3509 | INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J |
---|
3510 | KPP_REAL A,T |
---|
3511 | DIMENSION A(NDIM,N), IP(N) |
---|
3512 | C----------------------------------------------------------------------- |
---|
3513 | C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION. |
---|
3514 | C INPUT.. |
---|
3515 | C N = ORDER OF MATRIX. |
---|
3516 | C NDIM = DECLARED DIMENSION OF ARRAY A . |
---|
3517 | C A = MATRIX TO BE TRIANGULARIZED. |
---|
3518 | C OUTPUT.. |
---|
3519 | C A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U . |
---|
3520 | C A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. |
---|
3521 | C IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. |
---|
3522 | C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O . |
---|
3523 | C IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE |
---|
3524 | C SINGULAR AT STAGE K. |
---|
3525 | C USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. |
---|
3526 | C DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N). |
---|
3527 | C IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. |
---|
3528 | C |
---|
3529 | C REFERENCE.. |
---|
3530 | C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR, |
---|
3531 | C C.A.C.M. 15 (1972), P. 274. |
---|
3532 | C----------------------------------------------------------------------- |
---|
3533 | IER = 0 |
---|
3534 | IP(N) = 1 |
---|
3535 | IF (N .EQ. 1) GO TO 70 |
---|
3536 | NM1 = N - 1 |
---|
3537 | DO 60 K = 1,NM1 |
---|
3538 | KP1 = K + 1 |
---|
3539 | M = K |
---|
3540 | DO 10 I = KP1,N |
---|
3541 | IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I |
---|
3542 | 10 CONTINUE |
---|
3543 | IP(K) = M |
---|
3544 | T = A(M,K) |
---|
3545 | IF (M .EQ. K) GO TO 20 |
---|
3546 | IP(N) = -IP(N) |
---|
3547 | A(M,K) = A(K,K) |
---|
3548 | A(K,K) = T |
---|
3549 | 20 CONTINUE |
---|
3550 | IF (T .EQ. 0.D0) GO TO 80 |
---|
3551 | T = 1.D0/T |
---|
3552 | DO 30 I = KP1,N |
---|
3553 | 30 A(I,K) = -A(I,K)*T |
---|
3554 | DO 50 J = KP1,N |
---|
3555 | T = A(M,J) |
---|
3556 | A(M,J) = A(K,J) |
---|
3557 | A(K,J) = T |
---|
3558 | IF (T .EQ. 0.D0) GO TO 45 |
---|
3559 | DO 40 I = KP1,N |
---|
3560 | 40 A(I,J) = A(I,J) + A(I,K)*T |
---|
3561 | 45 CONTINUE |
---|
3562 | 50 CONTINUE |
---|
3563 | 60 CONTINUE |
---|
3564 | 70 K = N |
---|
3565 | IF (A(N,N) .EQ. 0.D0) GO TO 80 |
---|
3566 | RETURN |
---|
3567 | 80 IER = K |
---|
3568 | IP(N) = 0 |
---|
3569 | RETURN |
---|
3570 | C----------------------- END OF SUBROUTINE DEC ------------------------- |
---|
3571 | END |
---|
3572 | C |
---|
3573 | C |
---|
3574 | SUBROUTINE SOL (N, NDIM, A, B, IP) |
---|
3575 | C VERSION REAL KPP_REAL |
---|
3576 | INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1 |
---|
3577 | KPP_REAL A,B,T |
---|
3578 | DIMENSION A(NDIM,N), B(N), IP(N) |
---|
3579 | C----------------------------------------------------------------------- |
---|
3580 | C SOLUTION OF LINEAR SYSTEM, A*X = B . |
---|
3581 | C INPUT.. |
---|
3582 | C N = ORDER OF MATRIX. |
---|
3583 | C NDIM = DECLARED DIMENSION OF ARRAY A . |
---|
3584 | C A = TRIANGULARIZED MATRIX OBTAINED FROM DEC. |
---|
3585 | C B = RIGHT HAND SIDE VECTOR. |
---|
3586 | C IP = PIVOT VECTOR OBTAINED FROM DEC. |
---|
3587 | C DO NOT USE IF DEC HAS SET IER .NE. 0. |
---|
3588 | C OUTPUT.. |
---|
3589 | C B = SOLUTION VECTOR, X . |
---|
3590 | C----------------------------------------------------------------------- |
---|
3591 | IF (N .EQ. 1) GO TO 50 |
---|
3592 | NM1 = N - 1 |
---|
3593 | DO 20 K = 1,NM1 |
---|
3594 | KP1 = K + 1 |
---|
3595 | M = IP(K) |
---|
3596 | T = B(M) |
---|
3597 | B(M) = B(K) |
---|
3598 | B(K) = T |
---|
3599 | DO 10 I = KP1,N |
---|
3600 | 10 B(I) = B(I) + A(I,K)*T |
---|
3601 | 20 CONTINUE |
---|
3602 | DO 40 KB = 1,NM1 |
---|
3603 | KM1 = N - KB |
---|
3604 | K = KM1 + 1 |
---|
3605 | B(K) = B(K)/A(K,K) |
---|
3606 | T = -B(K) |
---|
3607 | DO 30 I = 1,KM1 |
---|
3608 | 30 B(I) = B(I) + A(I,K)*T |
---|
3609 | 40 CONTINUE |
---|
3610 | 50 B(1) = B(1)/A(1,1) |
---|
3611 | RETURN |
---|
3612 | C----------------------- END OF SUBROUTINE SOL ------------------------- |
---|
3613 | END |
---|
3614 | c |
---|
3615 | c |
---|
3616 | SUBROUTINE DECH (N, NDIM, A, LB, IP, IER) |
---|
3617 | C VERSION REAL KPP_REAL |
---|
3618 | INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J,LB,NA |
---|
3619 | KPP_REAL A,T |
---|
3620 | DIMENSION A(NDIM,N), IP(N) |
---|
3621 | C----------------------------------------------------------------------- |
---|
3622 | C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A HESSENBERG |
---|
3623 | C MATRIX WITH LOWER BANDWIDTH LB |
---|
3624 | C INPUT.. |
---|
3625 | C N = ORDER OF MATRIX A. |
---|
3626 | C NDIM = DECLARED DIMENSION OF ARRAY A . |
---|
3627 | C A = MATRIX TO BE TRIANGULARIZED. |
---|
3628 | C LB = LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED, LB.GE.1). |
---|
3629 | C OUTPUT.. |
---|
3630 | C A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U . |
---|
3631 | C A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. |
---|
3632 | C IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. |
---|
3633 | C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O . |
---|
3634 | C IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE |
---|
3635 | C SINGULAR AT STAGE K. |
---|
3636 | C USE SOLH TO OBTAIN SOLUTION OF LINEAR SYSTEM. |
---|
3637 | C DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N). |
---|
3638 | C IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. |
---|
3639 | C |
---|
3640 | C REFERENCE.. |
---|
3641 | C THIS IS A SLIGHT MODIFICATION OF |
---|
3642 | C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR, |
---|
3643 | C C.A.C.M. 15 (1972), P. 274. |
---|
3644 | C----------------------------------------------------------------------- |
---|
3645 | IER = 0 |
---|
3646 | IP(N) = 1 |
---|
3647 | IF (N .EQ. 1) GO TO 70 |
---|
3648 | NM1 = N - 1 |
---|
3649 | DO 60 K = 1,NM1 |
---|
3650 | KP1 = K + 1 |
---|
3651 | M = K |
---|
3652 | NA = MIN0(N,LB+K) |
---|
3653 | DO 10 I = KP1,NA |
---|
3654 | IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I |
---|
3655 | 10 CONTINUE |
---|
3656 | IP(K) = M |
---|
3657 | T = A(M,K) |
---|
3658 | IF (M .EQ. K) GO TO 20 |
---|
3659 | IP(N) = -IP(N) |
---|
3660 | A(M,K) = A(K,K) |
---|
3661 | A(K,K) = T |
---|
3662 | 20 CONTINUE |
---|
3663 | IF (T .EQ. 0.D0) GO TO 80 |
---|
3664 | T = 1.D0/T |
---|
3665 | DO 30 I = KP1,NA |
---|
3666 | 30 A(I,K) = -A(I,K)*T |
---|
3667 | DO 50 J = KP1,N |
---|
3668 | T = A(M,J) |
---|
3669 | A(M,J) = A(K,J) |
---|
3670 | A(K,J) = T |
---|
3671 | IF (T .EQ. 0.D0) GO TO 45 |
---|
3672 | DO 40 I = KP1,NA |
---|
3673 | 40 A(I,J) = A(I,J) + A(I,K)*T |
---|
3674 | 45 CONTINUE |
---|
3675 | 50 CONTINUE |
---|
3676 | 60 CONTINUE |
---|
3677 | 70 K = N |
---|
3678 | IF (A(N,N) .EQ. 0.D0) GO TO 80 |
---|
3679 | RETURN |
---|
3680 | 80 IER = K |
---|
3681 | IP(N) = 0 |
---|
3682 | RETURN |
---|
3683 | C----------------------- END OF SUBROUTINE DECH ------------------------ |
---|
3684 | END |
---|
3685 | C |
---|
3686 | C |
---|
3687 | SUBROUTINE SOLH (N, NDIM, A, LB, B, IP) |
---|
3688 | C VERSION REAL KPP_REAL |
---|
3689 | INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1,LB,NA |
---|
3690 | KPP_REAL A,B,T |
---|
3691 | DIMENSION A(NDIM,N), B(N), IP(N) |
---|
3692 | C----------------------------------------------------------------------- |
---|
3693 | C SOLUTION OF LINEAR SYSTEM, A*X = B . |
---|
3694 | C INPUT.. |
---|
3695 | C N = ORDER OF MATRIX A. |
---|
3696 | C NDIM = DECLARED DIMENSION OF ARRAY A . |
---|
3697 | C A = TRIANGULARIZED MATRIX OBTAINED FROM DECH. |
---|
3698 | C LB = LOWER BANDWIDTH OF A. |
---|
3699 | C B = RIGHT HAND SIDE VECTOR. |
---|
3700 | C IP = PIVOT VECTOR OBTAINED FROM DEC. |
---|
3701 | C DO NOT USE IF DECH HAS SET IER .NE. 0. |
---|
3702 | C OUTPUT.. |
---|
3703 | C B = SOLUTION VECTOR, X . |
---|
3704 | C----------------------------------------------------------------------- |
---|
3705 | IF (N .EQ. 1) GO TO 50 |
---|
3706 | NM1 = N - 1 |
---|
3707 | DO 20 K = 1,NM1 |
---|
3708 | KP1 = K + 1 |
---|
3709 | M = IP(K) |
---|
3710 | T = B(M) |
---|
3711 | B(M) = B(K) |
---|
3712 | B(K) = T |
---|
3713 | NA = MIN0(N,LB+K) |
---|
3714 | DO 10 I = KP1,NA |
---|
3715 | 10 B(I) = B(I) + A(I,K)*T |
---|
3716 | 20 CONTINUE |
---|
3717 | DO 40 KB = 1,NM1 |
---|
3718 | KM1 = N - KB |
---|
3719 | K = KM1 + 1 |
---|
3720 | B(K) = B(K)/A(K,K) |
---|
3721 | T = -B(K) |
---|
3722 | DO 30 I = 1,KM1 |
---|
3723 | 30 B(I) = B(I) + A(I,K)*T |
---|
3724 | 40 CONTINUE |
---|
3725 | 50 B(1) = B(1)/A(1,1) |
---|
3726 | RETURN |
---|
3727 | C----------------------- END OF SUBROUTINE SOLH ------------------------ |
---|
3728 | END |
---|
3729 | C |
---|
3730 | SUBROUTINE DECC (N, NDIM, AR, AI, IP, IER) |
---|
3731 | C VERSION COMPLEX KPP_REAL |
---|
3732 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
3733 | INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J |
---|
3734 | DIMENSION AR(NDIM,N), AI(NDIM,N), IP(N) |
---|
3735 | C----------------------------------------------------------------------- |
---|
3736 | C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION |
---|
3737 | C ------ MODIFICATION FOR COMPLEX MATRICES -------- |
---|
3738 | C INPUT.. |
---|
3739 | C N = ORDER OF MATRIX. |
---|
3740 | C NDIM = DECLARED DIMENSION OF ARRAYS AR AND AI . |
---|
3741 | C (AR, AI) = MATRIX TO BE TRIANGULARIZED. |
---|
3742 | C OUTPUT.. |
---|
3743 | C AR(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; REAL PART. |
---|
3744 | C AI(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; IMAGINARY PART. |
---|
3745 | C AR(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. |
---|
3746 | C REAL PART. |
---|
3747 | C AI(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. |
---|
3748 | C IMAGINARY PART. |
---|
3749 | C IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. |
---|
3750 | C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O . |
---|
3751 | C IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE |
---|
3752 | C SINGULAR AT STAGE K. |
---|
3753 | C USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. |
---|
3754 | C IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. |
---|
3755 | C |
---|
3756 | C REFERENCE.. |
---|
3757 | C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR, |
---|
3758 | C C.A.C.M. 15 (1972), P. 274. |
---|
3759 | C----------------------------------------------------------------------- |
---|
3760 | IER = 0 |
---|
3761 | IP(N) = 1 |
---|
3762 | IF (N .EQ. 1) GO TO 70 |
---|
3763 | NM1 = N - 1 |
---|
3764 | DO 60 K = 1,NM1 |
---|
3765 | KP1 = K + 1 |
---|
3766 | M = K |
---|
3767 | DO 10 I = KP1,N |
---|
3768 | IF (DABS(AR(I,K))+DABS(AI(I,K)) .GT. |
---|
3769 | & DABS(AR(M,K))+DABS(AI(M,K))) M = I |
---|
3770 | 10 CONTINUE |
---|
3771 | IP(K) = M |
---|
3772 | TR = AR(M,K) |
---|
3773 | TI = AI(M,K) |
---|
3774 | IF (M .EQ. K) GO TO 20 |
---|
3775 | IP(N) = -IP(N) |
---|
3776 | AR(M,K) = AR(K,K) |
---|
3777 | AI(M,K) = AI(K,K) |
---|
3778 | AR(K,K) = TR |
---|
3779 | AI(K,K) = TI |
---|
3780 | 20 CONTINUE |
---|
3781 | IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 80 |
---|
3782 | DEN=TR*TR+TI*TI |
---|
3783 | TR=TR/DEN |
---|
3784 | TI=-TI/DEN |
---|
3785 | DO 30 I = KP1,N |
---|
3786 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
3787 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
3788 | AR(I,K)=-PRODR |
---|
3789 | AI(I,K)=-PRODI |
---|
3790 | 30 CONTINUE |
---|
3791 | DO 50 J = KP1,N |
---|
3792 | TR = AR(M,J) |
---|
3793 | TI = AI(M,J) |
---|
3794 | AR(M,J) = AR(K,J) |
---|
3795 | AI(M,J) = AI(K,J) |
---|
3796 | AR(K,J) = TR |
---|
3797 | AI(K,J) = TI |
---|
3798 | IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 48 |
---|
3799 | IF (TI .EQ. 0.D0) THEN |
---|
3800 | DO 40 I = KP1,N |
---|
3801 | PRODR=AR(I,K)*TR |
---|
3802 | PRODI=AI(I,K)*TR |
---|
3803 | AR(I,J) = AR(I,J) + PRODR |
---|
3804 | AI(I,J) = AI(I,J) + PRODI |
---|
3805 | 40 CONTINUE |
---|
3806 | GO TO 48 |
---|
3807 | END IF |
---|
3808 | IF (TR .EQ. 0.D0) THEN |
---|
3809 | DO 45 I = KP1,N |
---|
3810 | PRODR=-AI(I,K)*TI |
---|
3811 | PRODI=AR(I,K)*TI |
---|
3812 | AR(I,J) = AR(I,J) + PRODR |
---|
3813 | AI(I,J) = AI(I,J) + PRODI |
---|
3814 | 45 CONTINUE |
---|
3815 | GO TO 48 |
---|
3816 | END IF |
---|
3817 | DO 47 I = KP1,N |
---|
3818 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
3819 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
3820 | AR(I,J) = AR(I,J) + PRODR |
---|
3821 | AI(I,J) = AI(I,J) + PRODI |
---|
3822 | 47 CONTINUE |
---|
3823 | 48 CONTINUE |
---|
3824 | 50 CONTINUE |
---|
3825 | 60 CONTINUE |
---|
3826 | 70 K = N |
---|
3827 | IF (DABS(AR(N,N))+DABS(AI(N,N)) .EQ. 0.D0) GO TO 80 |
---|
3828 | RETURN |
---|
3829 | 80 IER = K |
---|
3830 | IP(N) = 0 |
---|
3831 | RETURN |
---|
3832 | C----------------------- END OF SUBROUTINE DECC ------------------------ |
---|
3833 | END |
---|
3834 | C |
---|
3835 | C |
---|
3836 | SUBROUTINE SOLC (N, NDIM, AR, AI, BR, BI, IP) |
---|
3837 | C VERSION COMPLEX KPP_REAL |
---|
3838 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
3839 | INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1 |
---|
3840 | DIMENSION AR(NDIM,N), AI(NDIM,N), BR(N), BI(N), IP(N) |
---|
3841 | C----------------------------------------------------------------------- |
---|
3842 | C SOLUTION OF LINEAR SYSTEM, A*X = B . |
---|
3843 | C INPUT.. |
---|
3844 | C N = ORDER OF MATRIX. |
---|
3845 | C NDIM = DECLARED DIMENSION OF ARRAYS AR AND AI. |
---|
3846 | C (AR,AI) = TRIANGULARIZED MATRIX OBTAINED FROM DEC. |
---|
3847 | C (BR,BI) = RIGHT HAND SIDE VECTOR. |
---|
3848 | C IP = PIVOT VECTOR OBTAINED FROM DEC. |
---|
3849 | C DO NOT USE IF DEC HAS SET IER .NE. 0. |
---|
3850 | C OUTPUT.. |
---|
3851 | C (BR,BI) = SOLUTION VECTOR, X . |
---|
3852 | C----------------------------------------------------------------------- |
---|
3853 | IF (N .EQ. 1) GO TO 50 |
---|
3854 | NM1 = N - 1 |
---|
3855 | DO 20 K = 1,NM1 |
---|
3856 | KP1 = K + 1 |
---|
3857 | M = IP(K) |
---|
3858 | TR = BR(M) |
---|
3859 | TI = BI(M) |
---|
3860 | BR(M) = BR(K) |
---|
3861 | BI(M) = BI(K) |
---|
3862 | BR(K) = TR |
---|
3863 | BI(K) = TI |
---|
3864 | DO 10 I = KP1,N |
---|
3865 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
3866 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
3867 | BR(I) = BR(I) + PRODR |
---|
3868 | BI(I) = BI(I) + PRODI |
---|
3869 | 10 CONTINUE |
---|
3870 | 20 CONTINUE |
---|
3871 | DO 40 KB = 1,NM1 |
---|
3872 | KM1 = N - KB |
---|
3873 | K = KM1 + 1 |
---|
3874 | DEN=AR(K,K)*AR(K,K)+AI(K,K)*AI(K,K) |
---|
3875 | PRODR=BR(K)*AR(K,K)+BI(K)*AI(K,K) |
---|
3876 | PRODI=BI(K)*AR(K,K)-BR(K)*AI(K,K) |
---|
3877 | BR(K)=PRODR/DEN |
---|
3878 | BI(K)=PRODI/DEN |
---|
3879 | TR = -BR(K) |
---|
3880 | TI = -BI(K) |
---|
3881 | DO 30 I = 1,KM1 |
---|
3882 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
3883 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
3884 | BR(I) = BR(I) + PRODR |
---|
3885 | BI(I) = BI(I) + PRODI |
---|
3886 | 30 CONTINUE |
---|
3887 | 40 CONTINUE |
---|
3888 | 50 CONTINUE |
---|
3889 | DEN=AR(1,1)*AR(1,1)+AI(1,1)*AI(1,1) |
---|
3890 | PRODR=BR(1)*AR(1,1)+BI(1)*AI(1,1) |
---|
3891 | PRODI=BI(1)*AR(1,1)-BR(1)*AI(1,1) |
---|
3892 | BR(1)=PRODR/DEN |
---|
3893 | BI(1)=PRODI/DEN |
---|
3894 | RETURN |
---|
3895 | C----------------------- END OF SUBROUTINE SOLC ------------------------ |
---|
3896 | END |
---|
3897 | C |
---|
3898 | C |
---|
3899 | SUBROUTINE DECHC (N, NDIM, AR, AI, LB, IP, IER) |
---|
3900 | C VERSION COMPLEX KPP_REAL |
---|
3901 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
3902 | INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J |
---|
3903 | DIMENSION AR(NDIM,N), AI(NDIM,N), IP(N) |
---|
3904 | C----------------------------------------------------------------------- |
---|
3905 | C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION |
---|
3906 | C ------ MODIFICATION FOR COMPLEX MATRICES -------- |
---|
3907 | C INPUT.. |
---|
3908 | C N = ORDER OF MATRIX. |
---|
3909 | C NDIM = DECLARED DIMENSION OF ARRAYS AR AND AI . |
---|
3910 | C (AR, AI) = MATRIX TO BE TRIANGULARIZED. |
---|
3911 | C OUTPUT.. |
---|
3912 | C AR(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; REAL PART. |
---|
3913 | C AI(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; IMAGINARY PART. |
---|
3914 | C AR(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. |
---|
3915 | C REAL PART. |
---|
3916 | C AI(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. |
---|
3917 | C IMAGINARY PART. |
---|
3918 | C LB = LOWER BANDWIDTH OF A (DIAGONAL NOT COUNTED), LB.GE.1. |
---|
3919 | C IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. |
---|
3920 | C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O . |
---|
3921 | C IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE |
---|
3922 | C SINGULAR AT STAGE K. |
---|
3923 | C USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. |
---|
3924 | C IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. |
---|
3925 | C |
---|
3926 | C REFERENCE.. |
---|
3927 | C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR, |
---|
3928 | C C.A.C.M. 15 (1972), P. 274. |
---|
3929 | C----------------------------------------------------------------------- |
---|
3930 | IER = 0 |
---|
3931 | IP(N) = 1 |
---|
3932 | IF (LB .EQ. 0) GO TO 70 |
---|
3933 | IF (N .EQ. 1) GO TO 70 |
---|
3934 | NM1 = N - 1 |
---|
3935 | DO 60 K = 1,NM1 |
---|
3936 | KP1 = K + 1 |
---|
3937 | M = K |
---|
3938 | NA = MIN0(N,LB+K) |
---|
3939 | DO 10 I = KP1,NA |
---|
3940 | IF (DABS(AR(I,K))+DABS(AI(I,K)) .GT. |
---|
3941 | & DABS(AR(M,K))+DABS(AI(M,K))) M = I |
---|
3942 | 10 CONTINUE |
---|
3943 | IP(K) = M |
---|
3944 | TR = AR(M,K) |
---|
3945 | TI = AI(M,K) |
---|
3946 | IF (M .EQ. K) GO TO 20 |
---|
3947 | IP(N) = -IP(N) |
---|
3948 | AR(M,K) = AR(K,K) |
---|
3949 | AI(M,K) = AI(K,K) |
---|
3950 | AR(K,K) = TR |
---|
3951 | AI(K,K) = TI |
---|
3952 | 20 CONTINUE |
---|
3953 | IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 80 |
---|
3954 | DEN=TR*TR+TI*TI |
---|
3955 | TR=TR/DEN |
---|
3956 | TI=-TI/DEN |
---|
3957 | DO 30 I = KP1,NA |
---|
3958 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
3959 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
3960 | AR(I,K)=-PRODR |
---|
3961 | AI(I,K)=-PRODI |
---|
3962 | 30 CONTINUE |
---|
3963 | DO 50 J = KP1,N |
---|
3964 | TR = AR(M,J) |
---|
3965 | TI = AI(M,J) |
---|
3966 | AR(M,J) = AR(K,J) |
---|
3967 | AI(M,J) = AI(K,J) |
---|
3968 | AR(K,J) = TR |
---|
3969 | AI(K,J) = TI |
---|
3970 | IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 48 |
---|
3971 | IF (TI .EQ. 0.D0) THEN |
---|
3972 | DO 40 I = KP1,NA |
---|
3973 | PRODR=AR(I,K)*TR |
---|
3974 | PRODI=AI(I,K)*TR |
---|
3975 | AR(I,J) = AR(I,J) + PRODR |
---|
3976 | AI(I,J) = AI(I,J) + PRODI |
---|
3977 | 40 CONTINUE |
---|
3978 | GO TO 48 |
---|
3979 | END IF |
---|
3980 | IF (TR .EQ. 0.D0) THEN |
---|
3981 | DO 45 I = KP1,NA |
---|
3982 | PRODR=-AI(I,K)*TI |
---|
3983 | PRODI=AR(I,K)*TI |
---|
3984 | AR(I,J) = AR(I,J) + PRODR |
---|
3985 | AI(I,J) = AI(I,J) + PRODI |
---|
3986 | 45 CONTINUE |
---|
3987 | GO TO 48 |
---|
3988 | END IF |
---|
3989 | DO 47 I = KP1,NA |
---|
3990 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
3991 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
3992 | AR(I,J) = AR(I,J) + PRODR |
---|
3993 | AI(I,J) = AI(I,J) + PRODI |
---|
3994 | 47 CONTINUE |
---|
3995 | 48 CONTINUE |
---|
3996 | 50 CONTINUE |
---|
3997 | 60 CONTINUE |
---|
3998 | 70 K = N |
---|
3999 | IF (DABS(AR(N,N))+DABS(AI(N,N)) .EQ. 0.D0) GO TO 80 |
---|
4000 | RETURN |
---|
4001 | 80 IER = K |
---|
4002 | IP(N) = 0 |
---|
4003 | RETURN |
---|
4004 | C----------------------- END OF SUBROUTINE DECHC ----------------------- |
---|
4005 | END |
---|
4006 | C |
---|
4007 | C |
---|
4008 | SUBROUTINE SOLHC (N, NDIM, AR, AI, LB, BR, BI, IP) |
---|
4009 | C VERSION COMPLEX KPP_REAL |
---|
4010 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
4011 | INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1 |
---|
4012 | DIMENSION AR(NDIM,N), AI(NDIM,N), BR(N), BI(N), IP(N) |
---|
4013 | C----------------------------------------------------------------------- |
---|
4014 | C SOLUTION OF LINEAR SYSTEM, A*X = B . |
---|
4015 | C INPUT.. |
---|
4016 | C N = ORDER OF MATRIX. |
---|
4017 | C NDIM = DECLARED DIMENSION OF ARRAYS AR AND AI. |
---|
4018 | C (AR,AI) = TRIANGULARIZED MATRIX OBTAINED FROM DEC. |
---|
4019 | C (BR,BI) = RIGHT HAND SIDE VECTOR. |
---|
4020 | C LB = LOWER BANDWIDTH OF A. |
---|
4021 | C IP = PIVOT VECTOR OBTAINED FROM DEC. |
---|
4022 | C DO NOT USE IF DEC HAS SET IER .NE. 0. |
---|
4023 | C OUTPUT.. |
---|
4024 | C (BR,BI) = SOLUTION VECTOR, X . |
---|
4025 | C----------------------------------------------------------------------- |
---|
4026 | IF (N .EQ. 1) GO TO 50 |
---|
4027 | NM1 = N - 1 |
---|
4028 | IF (LB .EQ. 0) GO TO 25 |
---|
4029 | DO 20 K = 1,NM1 |
---|
4030 | KP1 = K + 1 |
---|
4031 | M = IP(K) |
---|
4032 | TR = BR(M) |
---|
4033 | TI = BI(M) |
---|
4034 | BR(M) = BR(K) |
---|
4035 | BI(M) = BI(K) |
---|
4036 | BR(K) = TR |
---|
4037 | BI(K) = TI |
---|
4038 | DO 10 I = KP1,MIN0(N,LB+K) |
---|
4039 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
4040 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
4041 | BR(I) = BR(I) + PRODR |
---|
4042 | BI(I) = BI(I) + PRODI |
---|
4043 | 10 CONTINUE |
---|
4044 | 20 CONTINUE |
---|
4045 | 25 CONTINUE |
---|
4046 | DO 40 KB = 1,NM1 |
---|
4047 | KM1 = N - KB |
---|
4048 | K = KM1 + 1 |
---|
4049 | DEN=AR(K,K)*AR(K,K)+AI(K,K)*AI(K,K) |
---|
4050 | PRODR=BR(K)*AR(K,K)+BI(K)*AI(K,K) |
---|
4051 | PRODI=BI(K)*AR(K,K)-BR(K)*AI(K,K) |
---|
4052 | BR(K)=PRODR/DEN |
---|
4053 | BI(K)=PRODI/DEN |
---|
4054 | TR = -BR(K) |
---|
4055 | TI = -BI(K) |
---|
4056 | DO 30 I = 1,KM1 |
---|
4057 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
4058 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
4059 | BR(I) = BR(I) + PRODR |
---|
4060 | BI(I) = BI(I) + PRODI |
---|
4061 | 30 CONTINUE |
---|
4062 | 40 CONTINUE |
---|
4063 | 50 CONTINUE |
---|
4064 | DEN=AR(1,1)*AR(1,1)+AI(1,1)*AI(1,1) |
---|
4065 | PRODR=BR(1)*AR(1,1)+BI(1)*AI(1,1) |
---|
4066 | PRODI=BI(1)*AR(1,1)-BR(1)*AI(1,1) |
---|
4067 | BR(1)=PRODR/DEN |
---|
4068 | BI(1)=PRODI/DEN |
---|
4069 | RETURN |
---|
4070 | C----------------------- END OF SUBROUTINE SOLHC ----------------------- |
---|
4071 | END |
---|
4072 | C |
---|
4073 | SUBROUTINE DECB (N, NDIM, A, ML, MU, IP, IER) |
---|
4074 | KPP_REAL A,T |
---|
4075 | DIMENSION A(NDIM,N), IP(N) |
---|
4076 | C----------------------------------------------------------------------- |
---|
4077 | C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED |
---|
4078 | C MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU |
---|
4079 | C INPUT.. |
---|
4080 | C N ORDER OF THE ORIGINAL MATRIX A. |
---|
4081 | C NDIM DECLARED DIMENSION OF ARRAY A. |
---|
4082 | C A CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS |
---|
4083 | C OF THE MATRIX ARE STORED IN THE COLUMNS OF A AND |
---|
4084 | C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS |
---|
4085 | C ML+1 THROUGH 2*ML+MU+1 OF A. |
---|
4086 | C ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). |
---|
4087 | C MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). |
---|
4088 | C OUTPUT.. |
---|
4089 | C A AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND |
---|
4090 | C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. |
---|
4091 | C IP INDEX VECTOR OF PIVOT INDICES. |
---|
4092 | C IP(N) (-1)**(NUMBER OF INTERCHANGES) OR O . |
---|
4093 | C IER = 0 IF MATRIX A IS NONSINGULAR, OR = K IF FOUND TO BE |
---|
4094 | C SINGULAR AT STAGE K. |
---|
4095 | C USE SOLB TO OBTAIN SOLUTION OF LINEAR SYSTEM. |
---|
4096 | C DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N) WITH MD=ML+MU+1. |
---|
4097 | C IF IP(N)=O, A IS SINGULAR, SOLB WILL DIVIDE BY ZERO. |
---|
4098 | C |
---|
4099 | C REFERENCE.. |
---|
4100 | C THIS IS A MODIFICATION OF |
---|
4101 | C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR, |
---|
4102 | C C.A.C.M. 15 (1972), P. 274. |
---|
4103 | C----------------------------------------------------------------------- |
---|
4104 | IER = 0 |
---|
4105 | IP(N) = 1 |
---|
4106 | MD = ML + MU + 1 |
---|
4107 | MD1 = MD + 1 |
---|
4108 | JU = 0 |
---|
4109 | IF (ML .EQ. 0) GO TO 70 |
---|
4110 | IF (N .EQ. 1) GO TO 70 |
---|
4111 | IF (N .LT. MU+2) GO TO 7 |
---|
4112 | DO 5 J = MU+2,N |
---|
4113 | DO 5 I = 1,ML |
---|
4114 | 5 A(I,J) = 0.D0 |
---|
4115 | 7 NM1 = N - 1 |
---|
4116 | DO 60 K = 1,NM1 |
---|
4117 | KP1 = K + 1 |
---|
4118 | M = MD |
---|
4119 | MDL = MIN(ML,N-K) + MD |
---|
4120 | DO 10 I = MD1,MDL |
---|
4121 | IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I |
---|
4122 | 10 CONTINUE |
---|
4123 | IP(K) = M + K - MD |
---|
4124 | T = A(M,K) |
---|
4125 | IF (M .EQ. MD) GO TO 20 |
---|
4126 | IP(N) = -IP(N) |
---|
4127 | A(M,K) = A(MD,K) |
---|
4128 | A(MD,K) = T |
---|
4129 | 20 CONTINUE |
---|
4130 | IF (T .EQ. 0.D0) GO TO 80 |
---|
4131 | T = 1.D0/T |
---|
4132 | DO 30 I = MD1,MDL |
---|
4133 | 30 A(I,K) = -A(I,K)*T |
---|
4134 | JU = MIN0(MAX0(JU,MU+IP(K)),N) |
---|
4135 | MM = MD |
---|
4136 | IF (JU .LT. KP1) GO TO 55 |
---|
4137 | DO 50 J = KP1,JU |
---|
4138 | M = M - 1 |
---|
4139 | MM = MM - 1 |
---|
4140 | T = A(M,J) |
---|
4141 | IF (M .EQ. MM) GO TO 35 |
---|
4142 | A(M,J) = A(MM,J) |
---|
4143 | A(MM,J) = T |
---|
4144 | 35 CONTINUE |
---|
4145 | IF (T .EQ. 0.D0) GO TO 45 |
---|
4146 | JK = J - K |
---|
4147 | DO 40 I = MD1,MDL |
---|
4148 | IJK = I - JK |
---|
4149 | 40 A(IJK,J) = A(IJK,J) + A(I,K)*T |
---|
4150 | 45 CONTINUE |
---|
4151 | 50 CONTINUE |
---|
4152 | 55 CONTINUE |
---|
4153 | 60 CONTINUE |
---|
4154 | 70 K = N |
---|
4155 | IF (A(MD,N) .EQ. 0.D0) GO TO 80 |
---|
4156 | RETURN |
---|
4157 | 80 IER = K |
---|
4158 | IP(N) = 0 |
---|
4159 | RETURN |
---|
4160 | C----------------------- END OF SUBROUTINE DECB ------------------------ |
---|
4161 | END |
---|
4162 | C |
---|
4163 | C |
---|
4164 | SUBROUTINE SOLB (N, NDIM, A, ML, MU, B, IP) |
---|
4165 | KPP_REAL A,B,T |
---|
4166 | DIMENSION A(NDIM,N), B(N), IP(N) |
---|
4167 | C----------------------------------------------------------------------- |
---|
4168 | C SOLUTION OF LINEAR SYSTEM, A*X = B . |
---|
4169 | C INPUT.. |
---|
4170 | C N ORDER OF MATRIX A. |
---|
4171 | C NDIM DECLARED DIMENSION OF ARRAY A . |
---|
4172 | C A TRIANGULARIZED MATRIX OBTAINED FROM DECB. |
---|
4173 | C ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). |
---|
4174 | C MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). |
---|
4175 | C B RIGHT HAND SIDE VECTOR. |
---|
4176 | C IP PIVOT VECTOR OBTAINED FROM DECB. |
---|
4177 | C DO NOT USE IF DECB HAS SET IER .NE. 0. |
---|
4178 | C OUTPUT.. |
---|
4179 | C B SOLUTION VECTOR, X . |
---|
4180 | C----------------------------------------------------------------------- |
---|
4181 | MD = ML + MU + 1 |
---|
4182 | MD1 = MD + 1 |
---|
4183 | MDM = MD - 1 |
---|
4184 | NM1 = N - 1 |
---|
4185 | IF (ML .EQ. 0) GO TO 25 |
---|
4186 | IF (N .EQ. 1) GO TO 50 |
---|
4187 | DO 20 K = 1,NM1 |
---|
4188 | M = IP(K) |
---|
4189 | T = B(M) |
---|
4190 | B(M) = B(K) |
---|
4191 | B(K) = T |
---|
4192 | MDL = MIN(ML,N-K) + MD |
---|
4193 | DO 10 I = MD1,MDL |
---|
4194 | IMD = I + K - MD |
---|
4195 | 10 B(IMD) = B(IMD) + A(I,K)*T |
---|
4196 | 20 CONTINUE |
---|
4197 | 25 CONTINUE |
---|
4198 | DO 40 KB = 1,NM1 |
---|
4199 | K = N + 1 - KB |
---|
4200 | B(K) = B(K)/A(MD,K) |
---|
4201 | T = -B(K) |
---|
4202 | KMD = MD - K |
---|
4203 | LM = MAX0(1,KMD+1) |
---|
4204 | DO 30 I = LM,MDM |
---|
4205 | IMD = I - KMD |
---|
4206 | 30 B(IMD) = B(IMD) + A(I,K)*T |
---|
4207 | 40 CONTINUE |
---|
4208 | 50 B(1) = B(1)/A(MD,1) |
---|
4209 | RETURN |
---|
4210 | C----------------------- END OF SUBROUTINE SOLB ------------------------ |
---|
4211 | END |
---|
4212 | C |
---|
4213 | SUBROUTINE DECBC (N, NDIM, AR, AI, ML, MU, IP, IER) |
---|
4214 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
4215 | DIMENSION AR(NDIM,N), AI(NDIM,N), IP(N) |
---|
4216 | C----------------------------------------------------------------------- |
---|
4217 | C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED COMPLEX |
---|
4218 | C MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU |
---|
4219 | C INPUT.. |
---|
4220 | C N ORDER OF THE ORIGINAL MATRIX A. |
---|
4221 | C NDIM DECLARED DIMENSION OF ARRAY A. |
---|
4222 | C AR, AI CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS |
---|
4223 | C OF THE MATRIX ARE STORED IN THE COLUMNS OF AR (REAL |
---|
4224 | C PART) AND AI (IMAGINARY PART) AND |
---|
4225 | C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS |
---|
4226 | C ML+1 THROUGH 2*ML+MU+1 OF AR AND AI. |
---|
4227 | C ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). |
---|
4228 | C MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). |
---|
4229 | C OUTPUT.. |
---|
4230 | C AR, AI AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND |
---|
4231 | C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. |
---|
4232 | C IP INDEX VECTOR OF PIVOT INDICES. |
---|
4233 | C IP(N) (-1)**(NUMBER OF INTERCHANGES) OR O . |
---|
4234 | C IER = 0 IF MATRIX A IS NONSINGULAR, OR = K IF FOUND TO BE |
---|
4235 | C SINGULAR AT STAGE K. |
---|
4236 | C USE SOLBC TO OBTAIN SOLUTION OF LINEAR SYSTEM. |
---|
4237 | C DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N) WITH MD=ML+MU+1. |
---|
4238 | C IF IP(N)=O, A IS SINGULAR, SOLBC WILL DIVIDE BY ZERO. |
---|
4239 | C |
---|
4240 | C REFERENCE.. |
---|
4241 | C THIS IS A MODIFICATION OF |
---|
4242 | C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR, |
---|
4243 | C C.A.C.M. 15 (1972), P. 274. |
---|
4244 | C----------------------------------------------------------------------- |
---|
4245 | IER = 0 |
---|
4246 | IP(N) = 1 |
---|
4247 | MD = ML + MU + 1 |
---|
4248 | MD1 = MD + 1 |
---|
4249 | JU = 0 |
---|
4250 | IF (ML .EQ. 0) GO TO 70 |
---|
4251 | IF (N .EQ. 1) GO TO 70 |
---|
4252 | IF (N .LT. MU+2) GO TO 7 |
---|
4253 | DO 5 J = MU+2,N |
---|
4254 | DO 5 I = 1,ML |
---|
4255 | AR(I,J) = 0.D0 |
---|
4256 | AI(I,J) = 0.D0 |
---|
4257 | 5 CONTINUE |
---|
4258 | 7 NM1 = N - 1 |
---|
4259 | DO 60 K = 1,NM1 |
---|
4260 | KP1 = K + 1 |
---|
4261 | M = MD |
---|
4262 | MDL = MIN(ML,N-K) + MD |
---|
4263 | DO 10 I = MD1,MDL |
---|
4264 | IF (DABS(AR(I,K))+DABS(AI(I,K)) .GT. |
---|
4265 | & DABS(AR(M,K))+DABS(AI(M,K))) M = I |
---|
4266 | 10 CONTINUE |
---|
4267 | IP(K) = M + K - MD |
---|
4268 | TR = AR(M,K) |
---|
4269 | TI = AI(M,K) |
---|
4270 | IF (M .EQ. MD) GO TO 20 |
---|
4271 | IP(N) = -IP(N) |
---|
4272 | AR(M,K) = AR(MD,K) |
---|
4273 | AI(M,K) = AI(MD,K) |
---|
4274 | AR(MD,K) = TR |
---|
4275 | AI(MD,K) = TI |
---|
4276 | 20 IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 80 |
---|
4277 | DEN=TR*TR+TI*TI |
---|
4278 | TR=TR/DEN |
---|
4279 | TI=-TI/DEN |
---|
4280 | DO 30 I = MD1,MDL |
---|
4281 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
4282 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
4283 | AR(I,K)=-PRODR |
---|
4284 | AI(I,K)=-PRODI |
---|
4285 | 30 CONTINUE |
---|
4286 | JU = MIN0(MAX0(JU,MU+IP(K)),N) |
---|
4287 | MM = MD |
---|
4288 | IF (JU .LT. KP1) GO TO 55 |
---|
4289 | DO 50 J = KP1,JU |
---|
4290 | M = M - 1 |
---|
4291 | MM = MM - 1 |
---|
4292 | TR = AR(M,J) |
---|
4293 | TI = AI(M,J) |
---|
4294 | IF (M .EQ. MM) GO TO 35 |
---|
4295 | AR(M,J) = AR(MM,J) |
---|
4296 | AI(M,J) = AI(MM,J) |
---|
4297 | AR(MM,J) = TR |
---|
4298 | AI(MM,J) = TI |
---|
4299 | 35 CONTINUE |
---|
4300 | IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 48 |
---|
4301 | JK = J - K |
---|
4302 | IF (TI .EQ. 0.D0) THEN |
---|
4303 | DO 40 I = MD1,MDL |
---|
4304 | IJK = I - JK |
---|
4305 | PRODR=AR(I,K)*TR |
---|
4306 | PRODI=AI(I,K)*TR |
---|
4307 | AR(IJK,J) = AR(IJK,J) + PRODR |
---|
4308 | AI(IJK,J) = AI(IJK,J) + PRODI |
---|
4309 | 40 CONTINUE |
---|
4310 | GO TO 48 |
---|
4311 | END IF |
---|
4312 | IF (TR .EQ. 0.D0) THEN |
---|
4313 | DO 45 I = MD1,MDL |
---|
4314 | IJK = I - JK |
---|
4315 | PRODR=-AI(I,K)*TI |
---|
4316 | PRODI=AR(I,K)*TI |
---|
4317 | AR(IJK,J) = AR(IJK,J) + PRODR |
---|
4318 | AI(IJK,J) = AI(IJK,J) + PRODI |
---|
4319 | 45 CONTINUE |
---|
4320 | GO TO 48 |
---|
4321 | END IF |
---|
4322 | DO 47 I = MD1,MDL |
---|
4323 | IJK = I - JK |
---|
4324 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
4325 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
4326 | AR(IJK,J) = AR(IJK,J) + PRODR |
---|
4327 | AI(IJK,J) = AI(IJK,J) + PRODI |
---|
4328 | 47 CONTINUE |
---|
4329 | 48 CONTINUE |
---|
4330 | 50 CONTINUE |
---|
4331 | 55 CONTINUE |
---|
4332 | 60 CONTINUE |
---|
4333 | 70 K = N |
---|
4334 | IF (DABS(AR(MD,N))+DABS(AI(MD,N)) .EQ. 0.D0) GO TO 80 |
---|
4335 | RETURN |
---|
4336 | 80 IER = K |
---|
4337 | IP(N) = 0 |
---|
4338 | RETURN |
---|
4339 | C----------------------- END OF SUBROUTINE DECBC ------------------------ |
---|
4340 | END |
---|
4341 | C |
---|
4342 | C |
---|
4343 | SUBROUTINE SOLBC (N, NDIM, AR, AI, ML, MU, BR, BI, IP) |
---|
4344 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
4345 | DIMENSION AR(NDIM,N), AI(NDIM,N), BR(N), BI(N), IP(N) |
---|
4346 | C----------------------------------------------------------------------- |
---|
4347 | C SOLUTION OF LINEAR SYSTEM, A*X = B , |
---|
4348 | C VERSION BANDED AND COMPLEX-KPP_REAL. |
---|
4349 | C INPUT.. |
---|
4350 | C N ORDER OF MATRIX A. |
---|
4351 | C NDIM DECLARED DIMENSION OF ARRAY A . |
---|
4352 | C AR, AI TRIANGULARIZED MATRIX OBTAINED FROM DECB (REAL AND IMAG. PART). |
---|
4353 | C ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). |
---|
4354 | C MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). |
---|
4355 | C BR, BI RIGHT HAND SIDE VECTOR (REAL AND IMAG. PART). |
---|
4356 | C IP PIVOT VECTOR OBTAINED FROM DECBC. |
---|
4357 | C DO NOT USE IF DECB HAS SET IER .NE. 0. |
---|
4358 | C OUTPUT.. |
---|
4359 | C BR, BI SOLUTION VECTOR, X (REAL AND IMAG. PART). |
---|
4360 | C----------------------------------------------------------------------- |
---|
4361 | MD = ML + MU + 1 |
---|
4362 | MD1 = MD + 1 |
---|
4363 | MDM = MD - 1 |
---|
4364 | NM1 = N - 1 |
---|
4365 | IF (ML .EQ. 0) GO TO 25 |
---|
4366 | IF (N .EQ. 1) GO TO 50 |
---|
4367 | DO 20 K = 1,NM1 |
---|
4368 | M = IP(K) |
---|
4369 | TR = BR(M) |
---|
4370 | TI = BI(M) |
---|
4371 | BR(M) = BR(K) |
---|
4372 | BI(M) = BI(K) |
---|
4373 | BR(K) = TR |
---|
4374 | BI(K) = TI |
---|
4375 | MDL = MIN(ML,N-K) + MD |
---|
4376 | DO 10 I = MD1,MDL |
---|
4377 | IMD = I + K - MD |
---|
4378 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
4379 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
4380 | BR(IMD) = BR(IMD) + PRODR |
---|
4381 | BI(IMD) = BI(IMD) + PRODI |
---|
4382 | 10 CONTINUE |
---|
4383 | 20 CONTINUE |
---|
4384 | 25 CONTINUE |
---|
4385 | DO 40 KB = 1,NM1 |
---|
4386 | K = N + 1 - KB |
---|
4387 | DEN=AR(MD,K)*AR(MD,K)+AI(MD,K)*AI(MD,K) |
---|
4388 | PRODR=BR(K)*AR(MD,K)+BI(K)*AI(MD,K) |
---|
4389 | PRODI=BI(K)*AR(MD,K)-BR(K)*AI(MD,K) |
---|
4390 | BR(K)=PRODR/DEN |
---|
4391 | BI(K)=PRODI/DEN |
---|
4392 | TR = -BR(K) |
---|
4393 | TI = -BI(K) |
---|
4394 | KMD = MD - K |
---|
4395 | LM = MAX0(1,KMD+1) |
---|
4396 | DO 30 I = LM,MDM |
---|
4397 | IMD = I - KMD |
---|
4398 | PRODR=AR(I,K)*TR-AI(I,K)*TI |
---|
4399 | PRODI=AI(I,K)*TR+AR(I,K)*TI |
---|
4400 | BR(IMD) = BR(IMD) + PRODR |
---|
4401 | BI(IMD) = BI(IMD) + PRODI |
---|
4402 | 30 CONTINUE |
---|
4403 | 40 CONTINUE |
---|
4404 | DEN=AR(MD,1)*AR(MD,1)+AI(MD,1)*AI(MD,1) |
---|
4405 | PRODR=BR(1)*AR(MD,1)+BI(1)*AI(MD,1) |
---|
4406 | PRODI=BI(1)*AR(MD,1)-BR(1)*AI(MD,1) |
---|
4407 | BR(1)=PRODR/DEN |
---|
4408 | BI(1)=PRODI/DEN |
---|
4409 | 50 CONTINUE |
---|
4410 | RETURN |
---|
4411 | C----------------------- END OF SUBROUTINE SOLBC ------------------------ |
---|
4412 | END |
---|
4413 | c |
---|
4414 | C |
---|
4415 | subroutine Elmhes(nm,n,low,igh,a,int) |
---|
4416 | C |
---|
4417 | integer i,j,m,n,la,nm,igh,kp1,low,mm1,mp1 |
---|
4418 | real*8 a(nm,n) |
---|
4419 | real*8 x,y |
---|
4420 | real*8 dabs |
---|
4421 | integer int(igh) |
---|
4422 | C |
---|
4423 | C this subroutine is a translation of the algol procedure elmhes, |
---|
4424 | C num. math. 12, 349-368(1968) by martin and wilkinson. |
---|
4425 | C handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). |
---|
4426 | C |
---|
4427 | C given a real general matrix, this subroutine |
---|
4428 | C reduces a submatrix situated in rows and columns |
---|
4429 | C low through igh to upper hessenberg form by |
---|
4430 | C stabilized elementary similarity transformations. |
---|
4431 | C |
---|
4432 | C on input: |
---|
4433 | C |
---|
4434 | C nm must be set to the row dimension of two-dimensional |
---|
4435 | C array parameters as declared in the calling program |
---|
4436 | C dimension statement; |
---|
4437 | C |
---|
4438 | C n is the order of the matrix; |
---|
4439 | C |
---|
4440 | C low and igh are integers determined by the balancing |
---|
4441 | C subroutine balanc. if balanc has not been used, |
---|
4442 | C set low=1, igh=n; |
---|
4443 | C |
---|
4444 | C a contains the input matrix. |
---|
4445 | C |
---|
4446 | C on output: |
---|
4447 | C |
---|
4448 | C a contains the hessenberg matrix. the multipliers |
---|
4449 | C which were used in the reduction are stored in the |
---|
4450 | C remaining triangle under the hessenberg matrix; |
---|
4451 | C |
---|
4452 | C int contains information on the rows and columns |
---|
4453 | C interchanged in the reduction. |
---|
4454 | C only elements low through igh are used. |
---|
4455 | C |
---|
4456 | C questions and comments should be directed to b. s. garbow, |
---|
4457 | C applied mathematics division, argonne national laboratory |
---|
4458 | C |
---|
4459 | C ------------------------------------------------------------------ |
---|
4460 | C |
---|
4461 | la = igh - 1 |
---|
4462 | kp1 = low + 1 |
---|
4463 | if (la .lt. kp1) go to 200 |
---|
4464 | C |
---|
4465 | do 180 m = kp1, la |
---|
4466 | mm1 = m - 1 |
---|
4467 | x = 0.0d0 |
---|
4468 | i = m |
---|
4469 | C |
---|
4470 | do 100 j = m, igh |
---|
4471 | if (dabs(a(j,mm1)) .le. dabs(x)) go to 100 |
---|
4472 | x = a(j,mm1) |
---|
4473 | i = j |
---|
4474 | 100 continue |
---|
4475 | C |
---|
4476 | int(m) = i |
---|
4477 | if (i .eq. m) go to 130 |
---|
4478 | C :::::::::: interchange rows and columns of a :::::::::: |
---|
4479 | do 110 j = mm1, n |
---|
4480 | y = a(i,j) |
---|
4481 | a(i,j) = a(m,j) |
---|
4482 | a(m,j) = y |
---|
4483 | 110 continue |
---|
4484 | C |
---|
4485 | do 120 j = 1, igh |
---|
4486 | y = a(j,i) |
---|
4487 | a(j,i) = a(j,m) |
---|
4488 | a(j,m) = y |
---|
4489 | 120 continue |
---|
4490 | C :::::::::: end interchange :::::::::: |
---|
4491 | 130 if (x .eq. 0.0d0) go to 180 |
---|
4492 | mp1 = m + 1 |
---|
4493 | C |
---|
4494 | do 160 i = mp1, igh |
---|
4495 | y = a(i,mm1) |
---|
4496 | if (y .eq. 0.0d0) go to 160 |
---|
4497 | y = y / x |
---|
4498 | a(i,mm1) = y |
---|
4499 | C |
---|
4500 | do 140 j = m, n |
---|
4501 | 140 a(i,j) = a(i,j) - y * a(m,j) |
---|
4502 | C |
---|
4503 | do 150 j = 1, igh |
---|
4504 | 150 a(j,m) = a(j,m) + y * a(j,i) |
---|
4505 | C |
---|
4506 | 160 continue |
---|
4507 | C |
---|
4508 | 180 continue |
---|
4509 | C |
---|
4510 | 200 return |
---|
4511 | C :::::::::: last card of elmhes :::::::::: |
---|
4512 | end |
---|
4513 | |
---|
4514 | SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN) |
---|
4515 | C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "CONTR5" |
---|
4516 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
4517 | DIMENSION Y(N),CONT(LRC) |
---|
4518 | COMMON /INTERN/XOUT |
---|
4519 | IF (NR.EQ.1) THEN |
---|
4520 | WRITE (6,99) X,Y(1),Y(2),NR-1 |
---|
4521 | XOUT=0.2D0 |
---|
4522 | ELSE |
---|
4523 | 10 CONTINUE |
---|
4524 | IF (X.GE.XOUT) THEN |
---|
4525 | C --- CONTINUOUS OUTPUT FOR RADAU5 |
---|
4526 | WRITE (6,99) XOUT,CONTR5(1,XOUT,CONT,LRC), |
---|
4527 | & CONTR5(2,XOUT,CONT,LRC),NR-1 |
---|
4528 | XOUT=XOUT+0.2D0 |
---|
4529 | GOTO 10 |
---|
4530 | END IF |
---|
4531 | END IF |
---|
4532 | 99 FORMAT(1X,'X =',F5.2,' Y =',2E18.10,' NSTEP =',I4) |
---|
4533 | RETURN |
---|
4534 | END |
---|
4535 | |
---|
4536 | SUBROUTINE FUNC_CHEM (N, T, V, FCT, RPAR, IPAR) |
---|
4537 | IMPLICIT NONE |
---|
4538 | INCLUDE 'KPP_ROOT_Parameters.h' |
---|
4539 | INCLUDE 'KPP_ROOT_Global.h' |
---|
4540 | KPP_REAL V(NVAR), FCT(NVAR) |
---|
4541 | KPP_REAL T, TOLD, RPAR(1) |
---|
4542 | INTEGER N, IPAR(1) |
---|
4543 | TOLD = TIME |
---|
4544 | TIME = T |
---|
4545 | CALL Update_SUN() |
---|
4546 | CALL Update_RCONST() |
---|
4547 | TIME = TOLD |
---|
4548 | CALL Fun(V, FIX, RCONST, FCT) |
---|
4549 | RETURN |
---|
4550 | END |
---|
4551 | |
---|
4552 | SUBROUTINE JAC_CHEM (N, T, V, JV, NROWPD, RPAR, IPAR) |
---|
4553 | IMPLICIT NONE |
---|
4554 | INCLUDE 'KPP_ROOT_Parameters.h' |
---|
4555 | INCLUDE 'KPP_ROOT_Global.h' |
---|
4556 | KPP_REAL V(NVAR), JV(NVAR,NVAR) |
---|
4557 | KPP_REAL T, TOLD, RPAR(1) |
---|
4558 | INTEGER N, IPAR(1), NROWPD, i, j |
---|
4559 | TOLD = TIME |
---|
4560 | TIME = T |
---|
4561 | CALL Update_SUN() |
---|
4562 | CALL Update_RCONST() |
---|
4563 | TIME = TOLD |
---|
4564 | DO i=1,NVAR |
---|
4565 | DO j=1,NVAR |
---|
4566 | JV(i,j) = 0.D0 |
---|
4567 | END DO |
---|
4568 | END DO |
---|
4569 | CALL Jac(V, FIX, RCONST, JV) |
---|
4570 | RETURN |
---|
4571 | END |
---|
4572 | |
---|