GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
sprepj.f
Go to the documentation of this file.
1  SUBROUTINE sprepj (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
2  1 F, JAC)
3 C***BEGIN PROLOGUE SPREPJ
4 C***SUBSIDIARY
5 C***PURPOSE Compute and process Newton iteration matrix.
6 C***TYPE SINGLE PRECISION (SPREPJ-S, DPREPJ-D)
7 C***AUTHOR Hindmarsh, Alan C., (LLNL)
8 C***DESCRIPTION
9 C
10 C SPREPJ is called by SSTODE to compute and process the matrix
11 C P = I - h*el(1)*J , where J is an approximation to the Jacobian.
12 C Here J is computed by the user-supplied routine JAC if
13 C MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5.
14 C If MITER = 3, a diagonal approximation to J is used.
15 C J is stored in WM and replaced by P. If MITER .ne. 3, P is then
16 C subjected to LU decomposition in preparation for later solution
17 C of linear systems with P as coefficient matrix. This is done
18 C by SGETRF if MITER = 1 or 2, and by SGBTRF if MITER = 4 or 5.
19 C
20 C In addition to variables described in SSTODE and SLSODE prologues,
21 C communication with SPREPJ uses the following:
22 C Y = array containing predicted values on entry.
23 C FTEM = work array of length N (ACOR in SSTODE).
24 C SAVF = array containing f evaluated at predicted y.
25 C WM = real work space for matrices. On output it contains the
26 C inverse diagonal matrix if MITER = 3 and the LU decomposition
27 C of P if MITER is 1, 2 , 4, or 5.
28 C Storage of matrix elements starts at WM(3).
29 C WM also contains the following matrix-related data:
30 C WM(1) = SQRT(UROUND), used in numerical Jacobian increments.
31 C WM(2) = H*EL0, saved for later use if MITER = 3.
32 C IWM = integer work space containing pivot information, starting at
33 C IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band
34 C parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
35 C EL0 = EL(1) (input).
36 C IERPJ = output error flag, = 0 if no trouble, .gt. 0 if
37 C P matrix found to be singular.
38 C JCUR = output flag = 1 to indicate that the Jacobian matrix
39 C (or approximation) is now current.
40 C This routine also uses the COMMON variables EL0, H, TN, UROUND,
41 C MITER, N, NFE, and NJE.
42 C
43 C***SEE ALSO SLSODE
44 C***ROUTINES CALLED SGBTRF, SGETRF, SVNORM
45 C***COMMON BLOCKS SLS001
46 C***REVISION HISTORY (YYMMDD)
47 C 791129 DATE WRITTEN
48 C 890501 Modified prologue to SLATEC/LDOC format. (FNF)
49 C 890504 Minor cosmetic changes. (FNF)
50 C 930809 Renamed to allow single/double precision versions. (ACH)
51 C 010412 Reduced size of Common block /SLS001/. (ACH)
52 C 031105 Restored 'own' variables to Common block /SLS001/, to
53 C enable interrupt/restart feature. (ACH)
54 C***END PROLOGUE SPREPJ
55 C**End
56  EXTERNAL f, jac
57  INTEGER NEQ, NYH, IWM
58  REAL Y, YH, EWT, FTEM, SAVF, WM
59  dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
60  1 wm(*), iwm(*)
61  INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, CNYH,
62  1 ialth, ipup, lmax, meo, nqnyh, nslp,
63  1 icf, ierpj, iersl, jcur, jstart, kflag, l,
64  2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
65  3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
66  REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
67  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
68  COMMON /sls001/ conit, crate, el(13), elco(13,12),
69  1 hold, rmax, tesco(3,12),
70  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
71  2 init, mxstep, mxhnil, nhnil, nslast, cnyh,
72  3 ialth, ipup, lmax, meo, nqnyh, nslp,
73  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
74  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
75  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
76  INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP,
77  1 mba, mband, meb1, meband, ml, ml3, mu, np1
78  REAL CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ,
79  1 svnorm
80 C
81 C***FIRST EXECUTABLE STATEMENT SPREPJ
82  nje = nje + 1
83  ierpj = 0
84  jcur = 1
85  hl0 = h*el0
86  GO TO (100, 200, 300, 400, 500), miter
87 C If MITER = 1, call JAC and multiply by scalar. -----------------------
88  100 lenp = n*n
89  DO 110 i = 1,lenp
90  110 wm(i+2) = 0.0e0
91  CALL jac (neq, tn, y, 0, 0, wm(3), n)
92  con = -hl0
93  DO 120 i = 1,lenp
94  120 wm(i+2) = wm(i+2)*con
95  GO TO 240
96 C If MITER = 2, make N calls to F to approximate J. --------------------
97  200 fac = svnorm(n, savf, ewt)
98  r0 = 1000.0e0*abs(h)*uround*n*fac
99  IF (r0 .EQ. 0.0e0) r0 = 1.0e0
100  srur = wm(1)
101  j1 = 2
102  DO 230 j = 1,n
103  yj = y(j)
104  r = max(srur*abs(yj),r0/ewt(j))
105  y(j) = y(j) + r
106  fac = -hl0/r
107  CALL f (neq, tn, y, ftem)
108  DO 220 i = 1,n
109  220 wm(i+j1) = (ftem(i) - savf(i))*fac
110  y(j) = yj
111  j1 = j1 + n
112  230 CONTINUE
113  nfe = nfe + n
114 C Add identity matrix. -------------------------------------------------
115  240 j = 3
116  np1 = n + 1
117  DO 250 i = 1,n
118  wm(j) = wm(j) + 1.0e0
119  250 j = j + np1
120 C Do LU decomposition on P. --------------------------------------------
121  CALL sgetrf (n, n, wm(3), n, iwm(21), ier)
122  IF (ier .NE. 0) ierpj = 1
123  RETURN
124 C If MITER = 3, construct a diagonal approximation to J and P. ---------
125  300 wm(2) = hl0
126  r = el0*0.1e0
127  DO 310 i = 1,n
128  310 y(i) = y(i) + r*(h*savf(i) - yh(i,2))
129  CALL f (neq, tn, y, wm(3))
130  nfe = nfe + 1
131  DO 320 i = 1,n
132  r0 = h*savf(i) - yh(i,2)
133  di = 0.1e0*r0 - h*(wm(i+2) - savf(i))
134  wm(i+2) = 1.0e0
135  IF (abs(r0) .LT. uround/ewt(i)) GO TO 320
136  IF (abs(di) .EQ. 0.0e0) GO TO 330
137  wm(i+2) = 0.1e0*r0/di
138  320 CONTINUE
139  RETURN
140  330 ierpj = 1
141  RETURN
142 C If MITER = 4, call JAC and multiply by scalar. -----------------------
143  400 ml = iwm(1)
144  mu = iwm(2)
145  ml3 = ml + 3
146  mband = ml + mu + 1
147  meband = mband + ml
148  lenp = meband*n
149  DO 410 i = 1,lenp
150  410 wm(i+2) = 0.0e0
151  CALL jac (neq, tn, y, ml, mu, wm(ml3), meband)
152  con = -hl0
153  DO 420 i = 1,lenp
154  420 wm(i+2) = wm(i+2)*con
155  GO TO 570
156 C If MITER = 5, make MBAND calls to F to approximate J. ----------------
157  500 ml = iwm(1)
158  mu = iwm(2)
159  mband = ml + mu + 1
160  mba = min(mband,n)
161  meband = mband + ml
162  meb1 = meband - 1
163  srur = wm(1)
164  fac = svnorm(n, savf, ewt)
165  r0 = 1000.0e0*abs(h)*uround*n*fac
166  IF (r0 .EQ. 0.0e0) r0 = 1.0e0
167  DO 560 j = 1,mba
168  DO 530 i = j,n,mband
169  yi = y(i)
170  r = max(srur*abs(yi),r0/ewt(i))
171  530 y(i) = y(i) + r
172  CALL f (neq, tn, y, ftem)
173  DO 550 jj = j,n,mband
174  y(jj) = yh(jj,1)
175  yjj = y(jj)
176  r = max(srur*abs(yjj),r0/ewt(jj))
177  fac = -hl0/r
178  i1 = max(jj-mu,1)
179  i2 = min(jj+ml,n)
180  ii = jj*meb1 - ml + 2
181  DO 540 i = i1,i2
182  540 wm(ii+i) = (ftem(i) - savf(i))*fac
183  550 CONTINUE
184  560 CONTINUE
185  nfe = nfe + mba
186 C Add identity matrix. -------------------------------------------------
187  570 ii = mband + 2
188  DO 580 i = 1,n
189  wm(ii) = wm(ii) + 1.0e0
190  580 ii = ii + meband
191 C Do LU decomposition of P. --------------------------------------------
192  CALL sgbtrf ( n, n, ml, mu, wm(3), meband, iwm(21), ier)
193  IF (ier .NE. 0) ierpj = 1
194  RETURN
195 C----------------------- END OF SUBROUTINE SPREPJ ----------------------
196  END
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
static T abs(T x)
Definition: pr-output.cc:1678
subroutine sprepj(NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC)
Definition: sprepj.f:3
real function svnorm(N, V, W)
Definition: svnorm.f:2