GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dmatd.f
Go to the documentation of this file.
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
59 C JACD, RES, DGETRF, DGBTRF
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
114 230 CALL dgetrf( neq, neq, wm, neq, iwm(lipvt), ier)
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
178 550 CALL dgbtrf(neq, neq, iwm(lml), iwm(lmu), wm, meband,
179  * iwm(lipvt), ier)
180  RETURN
181 C
182 C------END OF SUBROUTINE DMATD------------------------------------------
183  END
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230
subroutine dmatd(NEQ, X, Y, YPRIME, DELTA, CJ, H, IER, EWT, E, WM, IWM, RES, IRES, UROUND, JACD, RPAR, IPAR)
Definition: dmatd.f:7
static T abs(T x)
Definition: pr-output.cc:1678