GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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