|
3911
|
1 C Work performed under the auspices of the U.S. Department of Energy |
|
|
2 C by Lawrence Livermore National Laboratory under contract number |
|
|
3 C W-7405-Eng-48. |
|
|
4 C |
|
|
5 SUBROUTINE DMATD(NEQ,X,Y,YPRIME,DELTA,CJ,H,IER,EWT,E, |
|
|
6 * WM,IWM,RES,IRES,UROUND,JACD,RPAR,IPAR) |
|
|
7 C |
|
|
8 C***BEGIN PROLOGUE DMATD |
|
|
9 C***REFER TO DDASPK |
|
|
10 C***DATE WRITTEN 890101 (YYMMDD) |
|
|
11 C***REVISION DATE 900926 (YYMMDD) |
|
|
12 C***REVISION DATE 940701 (YYMMDD) (new LIPVT) |
|
|
13 C |
|
|
14 C----------------------------------------------------------------------- |
|
|
15 C***DESCRIPTION |
|
|
16 C |
|
|
17 C This routine computes the iteration matrix |
|
|
18 C J = dG/dY+CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0). |
|
|
19 C Here J is computed by: |
|
|
20 C the user-supplied routine JACD if IWM(MTYPE) is 1 or 4, or |
|
|
21 C by numerical difference quotients if IWM(MTYPE) is 2 or 5. |
|
|
22 C |
|
|
23 C The parameters have the following meanings. |
|
|
24 C X = Independent variable. |
|
|
25 C Y = Array containing predicted values. |
|
|
26 C YPRIME = Array containing predicted derivatives. |
|
|
27 C DELTA = Residual evaluated at (X,Y,YPRIME). |
|
|
28 C (Used only if IWM(MTYPE)=2 or 5). |
|
|
29 C CJ = Scalar parameter defining iteration matrix. |
|
|
30 C H = Current stepsize in integration. |
|
|
31 C IER = Variable which is .NE. 0 if iteration matrix |
|
|
32 C is singular, and 0 otherwise. |
|
|
33 C EWT = Vector of error weights for computing norms. |
|
|
34 C E = Work space (temporary) of length NEQ. |
|
|
35 C WM = Real work space for matrices. On output |
|
|
36 C it contains the LU decomposition |
|
|
37 C of the iteration matrix. |
|
|
38 C IWM = Integer work space containing |
|
|
39 C matrix information. |
|
|
40 C RES = External user-supplied subroutine |
|
|
41 C to evaluate the residual. See RES description |
|
|
42 C in DDASPK prologue. |
|
|
43 C IRES = Flag which is equal to zero if no illegal values |
|
|
44 C in RES, and less than zero otherwise. (If IRES |
|
|
45 C is less than zero, the matrix was not completed). |
|
|
46 C In this case (if IRES .LT. 0), then IER = 0. |
|
|
47 C UROUND = The unit roundoff error of the machine being used. |
|
|
48 C JACD = Name of the external user-supplied routine |
|
|
49 C to evaluate the iteration matrix. (This routine |
|
|
50 C is only used if IWM(MTYPE) is 1 or 4) |
|
|
51 C See JAC description for the case INFO(12) = 0 |
|
|
52 C in DDASPK prologue. |
|
|
53 C RPAR,IPAR= Real and integer parameter arrays that |
|
|
54 C are used for communication between the |
|
|
55 C calling program and external user routines. |
|
|
56 C They are not altered by DMATD. |
|
|
57 C----------------------------------------------------------------------- |
|
|
58 C***ROUTINES CALLED |
|
4329
|
59 C JACD, RES, DGETRF, DGBTRF |
|
3911
|
60 C |
|
|
61 C***END PROLOGUE DMATD |
|
|
62 C |
|
|
63 C |
|
|
64 IMPLICIT DOUBLE PRECISION(A-H,O-Z) |
|
|
65 DIMENSION Y(*),YPRIME(*),DELTA(*),EWT(*),E(*) |
|
|
66 DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*) |
|
|
67 EXTERNAL RES, JACD |
|
|
68 C |
|
|
69 PARAMETER (LML=1, LMU=2, LMTYPE=4, LNRE=12, LNPD=22, LLCIWP=30) |
|
|
70 C |
|
|
71 LIPVT = IWM(LLCIWP) |
|
|
72 IER = 0 |
|
|
73 MTYPE=IWM(LMTYPE) |
|
|
74 GO TO (100,200,300,400,500),MTYPE |
|
|
75 C |
|
|
76 C |
|
|
77 C Dense user-supplied matrix. |
|
|
78 C |
|
|
79 100 LENPD=IWM(LNPD) |
|
|
80 DO 110 I=1,LENPD |
|
|
81 110 WM(I)=0.0D0 |
|
|
82 CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR) |
|
|
83 GO TO 230 |
|
|
84 C |
|
|
85 C |
|
|
86 C Dense finite-difference-generated matrix. |
|
|
87 C |
|
|
88 200 IRES=0 |
|
|
89 NROW=0 |
|
|
90 SQUR = SQRT(UROUND) |
|
|
91 DO 210 I=1,NEQ |
|
|
92 DEL=SQUR*MAX(ABS(Y(I)),ABS(H*YPRIME(I)), |
|
|
93 * ABS(1.D0/EWT(I))) |
|
|
94 DEL=SIGN(DEL,H*YPRIME(I)) |
|
|
95 DEL=(Y(I)+DEL)-Y(I) |
|
|
96 YSAVE=Y(I) |
|
|
97 YPSAVE=YPRIME(I) |
|
|
98 Y(I)=Y(I)+DEL |
|
|
99 YPRIME(I)=YPRIME(I)+CJ*DEL |
|
|
100 IWM(LNRE)=IWM(LNRE)+1 |
|
|
101 CALL RES(X,Y,YPRIME,CJ,E,IRES,RPAR,IPAR) |
|
|
102 IF (IRES .LT. 0) RETURN |
|
|
103 DELINV=1.0D0/DEL |
|
|
104 DO 220 L=1,NEQ |
|
|
105 220 WM(NROW+L)=(E(L)-DELTA(L))*DELINV |
|
|
106 NROW=NROW+NEQ |
|
|
107 Y(I)=YSAVE |
|
|
108 YPRIME(I)=YPSAVE |
|
|
109 210 CONTINUE |
|
|
110 C |
|
|
111 C |
|
|
112 C Do dense-matrix LU decomposition on J. |
|
|
113 C |
|
4329
|
114 230 CALL DGETRF( NEQ, NEQ, WM, NEQ, IWM(LIPVT), IER) |
|
3911
|
115 RETURN |
|
|
116 C |
|
|
117 C |
|
|
118 C Dummy section for IWM(MTYPE)=3. |
|
|
119 C |
|
|
120 300 RETURN |
|
|
121 C |
|
|
122 C |
|
|
123 C Banded user-supplied matrix. |
|
|
124 C |
|
|
125 400 LENPD=IWM(LNPD) |
|
|
126 DO 410 I=1,LENPD |
|
|
127 410 WM(I)=0.0D0 |
|
|
128 CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR) |
|
|
129 MEBAND=2*IWM(LML)+IWM(LMU)+1 |
|
|
130 GO TO 550 |
|
|
131 C |
|
|
132 C |
|
|
133 C Banded finite-difference-generated matrix. |
|
|
134 C |
|
|
135 500 MBAND=IWM(LML)+IWM(LMU)+1 |
|
|
136 MBA=MIN0(MBAND,NEQ) |
|
|
137 MEBAND=MBAND+IWM(LML) |
|
|
138 MEB1=MEBAND-1 |
|
|
139 MSAVE=(NEQ/MBAND)+1 |
|
|
140 ISAVE=IWM(LNPD) |
|
|
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)), |
|
|
150 * ABS(1.D0/EWT(N))) |
|
|
151 DEL=SIGN(DEL,H*YPRIME(N)) |
|
|
152 DEL=(Y(N)+DEL)-Y(N) |
|
|
153 Y(N)=Y(N)+DEL |
|
|
154 510 YPRIME(N)=YPRIME(N)+CJ*DEL |
|
|
155 IWM(LNRE)=IWM(LNRE)+1 |
|
|
156 CALL RES(X,Y,YPRIME,CJ,E,IRES,RPAR,IPAR) |
|
|
157 IF (IRES .LT. 0) RETURN |
|
|
158 DO 530 N=J,NEQ,MBAND |
|
|
159 K= (N-J)/MBAND + 1 |
|
|
160 Y(N)=WM(ISAVE+K) |
|
|
161 YPRIME(N)=WM(IPSAVE+K) |
|
|
162 DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)), |
|
|
163 * ABS(1.D0/EWT(N))) |
|
|
164 DEL=SIGN(DEL,H*YPRIME(N)) |
|
|
165 DEL=(Y(N)+DEL)-Y(N) |
|
|
166 DELINV=1.0D0/DEL |
|
|
167 I1=MAX0(1,(N-IWM(LMU))) |
|
|
168 I2=MIN0(NEQ,(N+IWM(LML))) |
|
|
169 II=N*MEB1-IWM(LML) |
|
|
170 DO 520 I=I1,I2 |
|
|
171 520 WM(II+I)=(E(I)-DELTA(I))*DELINV |
|
|
172 530 CONTINUE |
|
|
173 540 CONTINUE |
|
|
174 C |
|
|
175 C |
|
|
176 C Do LU decomposition of banded J. |
|
|
177 C |
|
4329
|
178 550 CALL DGBTRF(NEQ, NEQ, IWM(LML), IWM(LMU), WM, MEBAND, |
|
|
179 * IWM(LIPVT), IER) |
|
3911
|
180 RETURN |
|
|
181 C |
|
|
182 C------END OF SUBROUTINE DMATD------------------------------------------ |
|
|
183 END |
