|
2329
|
1 SUBROUTINE DDAJAC (NEQ, X, Y, YPRIME, DELTA, CJ, H, |
|
|
2 + IER, WT, E, WM, IWM, RES, IRES, UROUND, JAC, RPAR, |
|
|
3 + IPAR, NTEMP) |
|
|
4 C***BEGIN PROLOGUE DDAJAC |
|
|
5 C***SUBSIDIARY |
|
|
6 C***PURPOSE Compute the iteration matrix for DDASSL and form the |
|
|
7 C LU-decomposition. |
|
|
8 C***LIBRARY SLATEC (DASSL) |
|
|
9 C***TYPE DOUBLE PRECISION (SDAJAC-S, DDAJAC-D) |
|
|
10 C***AUTHOR PETZOLD, LINDA R., (LLNL) |
|
|
11 C***DESCRIPTION |
|
|
12 C----------------------------------------------------------------------- |
|
|
13 C THIS ROUTINE COMPUTES THE ITERATION MATRIX |
|
|
14 C PD=DG/DY+CJ*DG/DYPRIME (WHERE G(X,Y,YPRIME)=0). |
|
|
15 C HERE PD IS COMPUTED BY THE USER-SUPPLIED |
|
|
16 C ROUTINE JAC IF IWM(MTYPE) IS 1 OR 4, AND |
|
|
17 C IT IS COMPUTED BY NUMERICAL FINITE DIFFERENCING |
|
|
18 C IF IWM(MTYPE)IS 2 OR 5 |
|
|
19 C THE PARAMETERS HAVE THE FOLLOWING MEANINGS. |
|
|
20 C Y = ARRAY CONTAINING PREDICTED VALUES |
|
|
21 C YPRIME = ARRAY CONTAINING PREDICTED DERIVATIVES |
|
|
22 C DELTA = RESIDUAL EVALUATED AT (X,Y,YPRIME) |
|
|
23 C (USED ONLY IF IWM(MTYPE)=2 OR 5) |
|
|
24 C CJ = SCALAR PARAMETER DEFINING ITERATION MATRIX |
|
|
25 C H = CURRENT STEPSIZE IN INTEGRATION |
|
|
26 C IER = VARIABLE WHICH IS .NE. 0 |
|
|
27 C IF ITERATION MATRIX IS SINGULAR, |
|
|
28 C AND 0 OTHERWISE. |
|
|
29 C WT = VECTOR OF WEIGHTS FOR COMPUTING NORMS |
|
|
30 C E = WORK SPACE (TEMPORARY) OF LENGTH NEQ |
|
|
31 C WM = REAL WORK SPACE FOR MATRICES. ON |
|
|
32 C OUTPUT IT CONTAINS THE LU DECOMPOSITION |
|
|
33 C OF THE ITERATION MATRIX. |
|
|
34 C IWM = INTEGER WORK SPACE CONTAINING |
|
|
35 C MATRIX INFORMATION |
|
|
36 C RES = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE |
|
|
37 C TO EVALUATE THE RESIDUAL FUNCTION G(X,Y,YPRIME) |
|
|
38 C IRES = FLAG WHICH IS EQUAL TO ZERO IF NO ILLEGAL VALUES |
|
|
39 C IN RES, AND LESS THAN ZERO OTHERWISE. (IF IRES |
|
|
40 C IS LESS THAN ZERO, THE MATRIX WAS NOT COMPLETED) |
|
|
41 C IN THIS CASE (IF IRES .LT. 0), THEN IER = 0. |
|
|
42 C UROUND = THE UNIT ROUNDOFF ERROR OF THE MACHINE BEING USED. |
|
|
43 C JAC = NAME OF THE EXTERNAL USER-SUPPLIED ROUTINE |
|
|
44 C TO EVALUATE THE ITERATION MATRIX (THIS ROUTINE |
|
|
45 C IS ONLY USED IF IWM(MTYPE) IS 1 OR 4) |
|
|
46 C----------------------------------------------------------------------- |
|
4329
|
47 C***ROUTINES CALLED DGBTRF, DGETRF |
|
2329
|
48 C***REVISION HISTORY (YYMMDD) |
|
|
49 C 830315 DATE WRITTEN |
|
|
50 C 901009 Finished conversion to SLATEC 4.0 format (F.N.Fritsch) |
|
|
51 C 901010 Modified three MAX calls to be all on one line. (FNF) |
|
|
52 C 901019 Merged changes made by C. Ulrich with SLATEC 4.0 format. |
|
|
53 C 901026 Added explicit declarations for all variables and minor |
|
|
54 C cosmetic changes to prologue. (FNF) |
|
|
55 C 901101 Corrected PURPOSE. (FNF) |
|
4329
|
56 C 020204 Convert to use LAPACK |
|
2329
|
57 C***END PROLOGUE DDAJAC |
|
|
58 C |
|
|
59 INTEGER NEQ, IER, IWM(*), IRES, IPAR(*), NTEMP |
|
|
60 DOUBLE PRECISION |
|
|
61 * X, Y(*), YPRIME(*), DELTA(*), CJ, H, WT(*), E(*), WM(*), |
|
|
62 * UROUND, RPAR(*) |
|
|
63 EXTERNAL RES, JAC |
|
|
64 C |
|
4329
|
65 EXTERNAL DGBTRF, DGETRF |
|
2329
|
66 C |
|
|
67 INTEGER I, I1, I2, II, IPSAVE, ISAVE, J, K, L, LENPD, LIPVT, |
|
|
68 * LML, LMTYPE, LMU, MBA, MBAND, MEB1, MEBAND, MSAVE, MTYPE, N, |
|
|
69 * NPD, NPDM1, NROW |
|
|
70 DOUBLE PRECISION DEL, DELINV, SQUR, YPSAVE, YSAVE |
|
|
71 C |
|
|
72 PARAMETER (NPD=1) |
|
|
73 PARAMETER (LML=1) |
|
|
74 PARAMETER (LMU=2) |
|
|
75 PARAMETER (LMTYPE=4) |
|
4429
|
76 PARAMETER (LIPVT=22) |
|
2329
|
77 C |
|
|
78 C***FIRST EXECUTABLE STATEMENT DDAJAC |
|
|
79 IER = 0 |
|
|
80 NPDM1=NPD-1 |
|
|
81 MTYPE=IWM(LMTYPE) |
|
|
82 GO TO (100,200,300,400,500),MTYPE |
|
|
83 C |
|
|
84 C |
|
|
85 C DENSE USER-SUPPLIED MATRIX |
|
|
86 100 LENPD=NEQ*NEQ |
|
|
87 DO 110 I=1,LENPD |
|
|
88 110 WM(NPDM1+I)=0.0D0 |
|
|
89 CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR) |
|
|
90 GO TO 230 |
|
|
91 C |
|
|
92 C |
|
|
93 C DENSE FINITE-DIFFERENCE-GENERATED MATRIX |
|
|
94 200 IRES=0 |
|
|
95 NROW=NPDM1 |
|
|
96 SQUR = SQRT(UROUND) |
|
|
97 DO 210 I=1,NEQ |
|
|
98 DEL=SQUR*MAX(ABS(Y(I)),ABS(H*YPRIME(I)),ABS(WT(I))) |
|
|
99 DEL=SIGN(DEL,H*YPRIME(I)) |
|
|
100 DEL=(Y(I)+DEL)-Y(I) |
|
|
101 YSAVE=Y(I) |
|
|
102 YPSAVE=YPRIME(I) |
|
|
103 Y(I)=Y(I)+DEL |
|
|
104 YPRIME(I)=YPRIME(I)+CJ*DEL |
|
|
105 CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR) |
|
|
106 IF (IRES .LT. 0) RETURN |
|
|
107 DELINV=1.0D0/DEL |
|
|
108 DO 220 L=1,NEQ |
|
|
109 220 WM(NROW+L)=(E(L)-DELTA(L))*DELINV |
|
|
110 NROW=NROW+NEQ |
|
|
111 Y(I)=YSAVE |
|
|
112 YPRIME(I)=YPSAVE |
|
|
113 210 CONTINUE |
|
|
114 C |
|
|
115 C |
|
|
116 C DO DENSE-MATRIX LU DECOMPOSITION ON PD |
|
4329
|
117 230 CALL DGETRF( NEQ, NEQ, WM(NPD), NEQ, IWM(LIPVT), IER) |
|
2329
|
118 RETURN |
|
|
119 C |
|
|
120 C |
|
|
121 C DUMMY SECTION FOR IWM(MTYPE)=3 |
|
|
122 300 RETURN |
|
|
123 C |
|
|
124 C |
|
|
125 C BANDED USER-SUPPLIED MATRIX |
|
|
126 400 LENPD=(2*IWM(LML)+IWM(LMU)+1)*NEQ |
|
|
127 DO 410 I=1,LENPD |
|
|
128 410 WM(NPDM1+I)=0.0D0 |
|
|
129 CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR) |
|
|
130 MEBAND=2*IWM(LML)+IWM(LMU)+1 |
|
|
131 GO TO 550 |
|
|
132 C |
|
|
133 C |
|
|
134 C BANDED FINITE-DIFFERENCE-GENERATED MATRIX |
|
|
135 500 MBAND=IWM(LML)+IWM(LMU)+1 |
|
|
136 MBA=MIN(MBAND,NEQ) |
|
|
137 MEBAND=MBAND+IWM(LML) |
|
|
138 MEB1=MEBAND-1 |
|
|
139 MSAVE=(NEQ/MBAND)+1 |
|
|
140 ISAVE=NTEMP-1 |
|
|
141 IPSAVE=ISAVE+MSAVE |
|
|
142 IRES=0 |
|
|
143 SQUR=SQRT(UROUND) |
|
|
144 DO 540 J=1,MBA |
|
|
145 DO 510 N=J,NEQ,MBAND |
|
|
146 K= (N-J)/MBAND + 1 |
|
|
147 WM(ISAVE+K)=Y(N) |
|
|
148 WM(IPSAVE+K)=YPRIME(N) |
|
|
149 DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),ABS(WT(N))) |
|
|
150 DEL=SIGN(DEL,H*YPRIME(N)) |
|
|
151 DEL=(Y(N)+DEL)-Y(N) |
|
|
152 Y(N)=Y(N)+DEL |
|
|
153 510 YPRIME(N)=YPRIME(N)+CJ*DEL |
|
|
154 CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR) |
|
|
155 IF (IRES .LT. 0) RETURN |
|
|
156 DO 530 N=J,NEQ,MBAND |
|
|
157 K= (N-J)/MBAND + 1 |
|
|
158 Y(N)=WM(ISAVE+K) |
|
|
159 YPRIME(N)=WM(IPSAVE+K) |
|
|
160 DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),ABS(WT(N))) |
|
|
161 DEL=SIGN(DEL,H*YPRIME(N)) |
|
|
162 DEL=(Y(N)+DEL)-Y(N) |
|
|
163 DELINV=1.0D0/DEL |
|
|
164 I1=MAX(1,(N-IWM(LMU))) |
|
|
165 I2=MIN(NEQ,(N+IWM(LML))) |
|
|
166 II=N*MEB1-IWM(LML)+NPDM1 |
|
|
167 DO 520 I=I1,I2 |
|
|
168 520 WM(II+I)=(E(I)-DELTA(I))*DELINV |
|
|
169 530 CONTINUE |
|
|
170 540 CONTINUE |
|
|
171 C |
|
|
172 C |
|
|
173 C DO LU DECOMPOSITION OF BANDED PD |
|
4329
|
174 550 CALL DGBTRF(NEQ, NEQ, IWM(LML), IWM(LMU), WM(NPD), MEBAND, |
|
|
175 * IWM(LIPVT), IER) |
|
2329
|
176 RETURN |
|
|
177 C------END OF SUBROUTINE DDAJAC------ |
|
|
178 END |
