GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
prepj.f
Go to the documentation of this file.
1  SUBROUTINE prepj (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
2  1 F, JAC, IERR)
3 CLLL. OPTIMIZE
4  EXTERNAL f, jac
5  INTEGER NEQ, NYH, IWM
6  INTEGER ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
7  1 mxstep, mxhnil, nhnil, ntrep, nslast, cnyh,
8  2 ialth, ipup, lmax, meo, nqnyh, nslp
9  INTEGER ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
10  2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
11  INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP,
12  1 mba, mband, meb1, meband, ml, ml3, mu, np1
13  DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WM
14  DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
15  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
16  DOUBLE PRECISION CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ,
17  1 vnorm
18  dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
19  1 wm(*), iwm(*)
20  COMMON /ls0001/ conit, crate, el(13), elco(13,12),
21  1 hold, rmax, tesco(3,12),
22  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
23  2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
24  3 mxstep, mxhnil, nhnil, ntrep, nslast, cnyh,
25  3 ialth, ipup, lmax, meo, nqnyh, nslp,
26  4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
27  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
28 C-----------------------------------------------------------------------
29 C PREPJ IS CALLED BY STODE TO COMPUTE AND PROCESS THE MATRIX
30 C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
31 C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF
32 C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
33 C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
34 C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN
35 C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
36 C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
37 C BY DGETRF IF MITER = 1 OR 2, AND BY DGBTRF IF MITER = 4 OR 5.
38 C
39 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
40 C WITH PREPJ USES THE FOLLOWING..
41 C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
42 C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STODE).
43 C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
44 C WM = REAL WORK SPACE FOR MATRICES. ON OUTPUT IT CONTAINS THE
45 C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
46 C OF P IF MITER IS 1, 2 , 4, OR 5.
47 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
48 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
49 C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
50 C WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
51 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
52 C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND
53 C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
54 C EL0 = EL(1) (INPUT).
55 C IERPJ = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .GT. 0 IF
56 C P MATRIX FOUND TO BE SINGULAR.
57 C JCUR = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX
58 C (OR APPROXIMATION) IS NOW CURRENT.
59 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
60 C MITER, N, NFE, AND NJE.
61 C-----------------------------------------------------------------------
62  nje = nje + 1
63  ierpj = 0
64  jcur = 1
65  hl0 = h*el0
66  GO TO (100, 200, 300, 400, 500), miter
67 C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
68  100 lenp = n*n
69  DO 110 i = 1,lenp
70  110 wm(i+2) = 0.0d0
71  CALL jac (neq, tn, y, 0, 0, wm(3), n)
72  con = -hl0
73  DO 120 i = 1,lenp
74  120 wm(i+2) = wm(i+2)*con
75  GO TO 240
76 C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
77  200 fac = vnorm(n, savf, ewt)
78  r0 = 1000.0d0*dabs(h)*uround*dble(n)*fac
79  IF (r0 .EQ. 0.0d0) r0 = 1.0d0
80  srur = wm(1)
81  j1 = 2
82  DO 230 j = 1,n
83  yj = y(j)
84  r = dmax1(srur*dabs(yj),r0/ewt(j))
85  y(j) = y(j) + r
86  fac = -hl0/r
87  ierr = 0
88  CALL f (neq, tn, y, ftem, ierr)
89  IF (ierr .LT. 0) RETURN
90  DO 220 i = 1,n
91  220 wm(i+j1) = (ftem(i) - savf(i))*fac
92  y(j) = yj
93  j1 = j1 + n
94  230 CONTINUE
95  nfe = nfe + n
96 C ADD IDENTITY MATRIX. -------------------------------------------------
97  240 j = 3
98  np1 = n + 1
99  DO 250 i = 1,n
100  wm(j) = wm(j) + 1.0d0
101  250 j = j + np1
102 C DO LU DECOMPOSITION ON P. --------------------------------------------
103  CALL dgetrf ( n, n, wm(3), n, iwm(21), ier)
104  IF (ier .NE. 0) ierpj = 1
105  RETURN
106 C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
107  300 wm(2) = hl0
108  r = el0*0.1d0
109  DO 310 i = 1,n
110  310 y(i) = y(i) + r*(h*savf(i) - yh(i,2))
111  ierr = 0
112  CALL f (neq, tn, y, wm(3), ierr)
113  IF (ierr .LT. 0) RETURN
114  nfe = nfe + 1
115  DO 320 i = 1,n
116  r0 = h*savf(i) - yh(i,2)
117  di = 0.1d0*r0 - h*(wm(i+2) - savf(i))
118  wm(i+2) = 1.0d0
119  IF (dabs(r0) .LT. uround/ewt(i)) GO TO 320
120  IF (dabs(di) .EQ. 0.0d0) GO TO 330
121  wm(i+2) = 0.1d0*r0/di
122  320 CONTINUE
123  RETURN
124  330 ierpj = 1
125  RETURN
126 C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
127  400 ml = iwm(1)
128  mu = iwm(2)
129  ml3 = ml + 3
130  mband = ml + mu + 1
131  meband = mband + ml
132  lenp = meband*n
133  DO 410 i = 1,lenp
134  410 wm(i+2) = 0.0d0
135  CALL jac (neq, tn, y, ml, mu, wm(ml3), meband)
136  con = -hl0
137  DO 420 i = 1,lenp
138  420 wm(i+2) = wm(i+2)*con
139  GO TO 570
140 C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
141  500 ml = iwm(1)
142  mu = iwm(2)
143  mband = ml + mu + 1
144  mba = min0(mband,n)
145  meband = mband + ml
146  meb1 = meband - 1
147  srur = wm(1)
148  fac = vnorm(n, savf, ewt)
149  r0 = 1000.0d0*dabs(h)*uround*dble(n)*fac
150  IF (r0 .EQ. 0.0d0) r0 = 1.0d0
151  DO 560 j = 1,mba
152  DO 530 i = j,n,mband
153  yi = y(i)
154  r = dmax1(srur*dabs(yi),r0/ewt(i))
155  530 y(i) = y(i) + r
156  ierr = 0
157  CALL f (neq, tn, y, ftem, ierr)
158  IF (ierr .LT. 0) RETURN
159  DO 550 jj = j,n,mband
160  y(jj) = yh(jj,1)
161  yjj = y(jj)
162  r = dmax1(srur*dabs(yjj),r0/ewt(jj))
163  fac = -hl0/r
164  i1 = max0(jj-mu,1)
165  i2 = min0(jj+ml,n)
166  ii = jj*meb1 - ml + 2
167  DO 540 i = i1,i2
168  540 wm(ii+i) = (ftem(i) - savf(i))*fac
169  550 CONTINUE
170  560 CONTINUE
171  nfe = nfe + mba
172 C ADD IDENTITY MATRIX. -------------------------------------------------
173  570 ii = mband + 2
174  DO 580 i = 1,n
175  wm(ii) = wm(ii) + 1.0d0
176  580 ii = ii + meband
177 C DO LU DECOMPOSITION OF P. --------------------------------------------
178  CALL dgbtrf ( n, n, ml, mu, wm(3), meband, iwm(21), ier)
179  IF (ier .NE. 0) ierpj = 1
180  RETURN
181 C----------------------- END OF SUBROUTINE PREPJ -----------------------
182  END
subroutine prepj(NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC, IERR)
Definition: prepj.f:3
double precision function vnorm(N, V, W)
Definition: vnorm.f:2