GNU Octave  4.0.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
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
subroutine ddajac(NEQ, X, Y, YPRIME, DELTA, CJ, H, IER, WT, E, WM, IWM, RES, IRES, UROUND, JAC, RPAR, IPAR, NTEMP)
Definition: ddajac.f:1
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:233
float_format & precision(int p)
Definition: pr-output.cc:176
T abs(T x)
Definition: pr-output.cc:3062
octave_value sqrt(void) const
Definition: ov.h:1200
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:210