GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
ddajac.f
Go to the documentation of this file.
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-----------------------------------------------------------------------
47 C***ROUTINES CALLED DGBTRF, DGETRF
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)
56 C 020204 Convert to use LAPACK
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
65  EXTERNAL dgbtrf, dgetrf
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)
76  parameter(lipvt=22)
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
117 230 CALL dgetrf( neq, neq, wm(npd), neq, iwm(lipvt), ier)
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
174 550 CALL dgbtrf(neq, neq, iwm(lml), iwm(lmu), wm(npd), meband,
175  * iwm(lipvt), ier)
176  RETURN
177 C------END OF SUBROUTINE DDAJAC------
178  END
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
subroutine ddajac(NEQ, X, Y, YPRIME, DELTA, CJ, H, IER, WT, E, WM, IWM, RES, IRES, UROUND, JAC, RPAR, IPAR, NTEMP)
Definition: ddajac.f:4
static T abs(T x)
Definition: pr-output.cc:1678