GNU Octave  4.0.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
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 IOWND, IOWNS,
62  1 icf, ierpj, iersl, jcur, jstart, kflag, l,
63  2 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
64  3 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
65  REAL ROWNS,
66  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
67  COMMON /sls001/ rowns(209),
68  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
69  2 iownd(6), iowns(6),
70  3 icf, ierpj, iersl, jcur, jstart, kflag, l,
71  4 lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
72  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
73  INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP,
74  1 mba, mband, meb1, meband, ml, ml3, mu, np1
75  REAL CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ,
76  1 svnorm
77 C
78 C***FIRST EXECUTABLE STATEMENT SPREPJ
79  nje = nje + 1
80  ierpj = 0
81  jcur = 1
82  hl0 = h*el0
83  go to(100, 200, 300, 400, 500), miter
84 C If MITER = 1, call JAC and multiply by scalar. -----------------------
85  100 lenp = n*n
86  DO 110 i = 1,lenp
87  110 wm(i+2) = 0.0e0
88  CALL jac(neq, tn, y, 0, 0, wm(3), n)
89  con = -hl0
90  DO 120 i = 1,lenp
91  120 wm(i+2) = wm(i+2)*con
92  go to 240
93 C If MITER = 2, make N calls to F to approximate J. --------------------
94  200 fac = svnorm(n, savf, ewt)
95  r0 = 1000.0e0*abs(h)*uround*n*fac
96  IF (r0 .EQ. 0.0e0) r0 = 1.0e0
97  srur = wm(1)
98  j1 = 2
99  DO 230 j = 1,n
100  yj = y(j)
101  r = max(srur*abs(yj),r0/ewt(j))
102  y(j) = y(j) + r
103  fac = -hl0/r
104  CALL f(neq, tn, y, ftem)
105  DO 220 i = 1,n
106  220 wm(i+j1) = (ftem(i) - savf(i))*fac
107  y(j) = yj
108  j1 = j1 + n
109  230 CONTINUE
110  nfe = nfe + n
111 C Add identity matrix. -------------------------------------------------
112  240 j = 3
113  np1 = n + 1
114  DO 250 i = 1,n
115  wm(j) = wm(j) + 1.0e0
116  250 j = j + np1
117 C Do LU decomposition on P. --------------------------------------------
118  CALL sgetrf(n, n, wm(3), n, iwm(21), ier)
119  IF (ier .NE. 0) ierpj = 1
120  RETURN
121 C If MITER = 3, construct a diagonal approximation to J and P. ---------
122  300 wm(2) = hl0
123  r = el0*0.1e0
124  DO 310 i = 1,n
125  310 y(i) = y(i) + r*(h*savf(i) - yh(i,2))
126  CALL f(neq, tn, y, wm(3))
127  nfe = nfe + 1
128  DO 320 i = 1,n
129  r0 = h*savf(i) - yh(i,2)
130  di = 0.1e0*r0 - h*(wm(i+2) - savf(i))
131  wm(i+2) = 1.0e0
132  IF (abs(r0) .LT. uround/ewt(i)) go to 320
133  IF (abs(di) .EQ. 0.0e0) go to 330
134  wm(i+2) = 0.1e0*r0/di
135  320 CONTINUE
136  RETURN
137  330 ierpj = 1
138  RETURN
139 C If MITER = 4, call JAC and multiply by scalar. -----------------------
140  400 ml = iwm(1)
141  mu = iwm(2)
142  ml3 = ml + 3
143  mband = ml + mu + 1
144  meband = mband + ml
145  lenp = meband*n
146  DO 410 i = 1,lenp
147  410 wm(i+2) = 0.0e0
148  CALL jac(neq, tn, y, ml, mu, wm(ml3), meband)
149  con = -hl0
150  DO 420 i = 1,lenp
151  420 wm(i+2) = wm(i+2)*con
152  go to 570
153 C If MITER = 5, make MBAND calls to F to approximate J. ----------------
154  500 ml = iwm(1)
155  mu = iwm(2)
156  mband = ml + mu + 1
157  mba = min(mband,n)
158  meband = mband + ml
159  meb1 = meband - 1
160  srur = wm(1)
161  fac = svnorm(n, savf, ewt)
162  r0 = 1000.0e0*abs(h)*uround*n*fac
163  IF (r0 .EQ. 0.0e0) r0 = 1.0e0
164  DO 560 j = 1,mba
165  DO 530 i = j,n,mband
166  yi = y(i)
167  r = max(srur*abs(yi),r0/ewt(i))
168  530 y(i) = y(i) + r
169  CALL f(neq, tn, y, ftem)
170  DO 550 jj = j,n,mband
171  y(jj) = yh(jj,1)
172  yjj = y(jj)
173  r = max(srur*abs(yjj),r0/ewt(jj))
174  fac = -hl0/r
175  i1 = max(jj-mu,1)
176  i2 = min(jj+ml,n)
177  ii = jj*meb1 - ml + 2
178  DO 540 i = i1,i2
179  540 wm(ii+i) = (ftem(i) - savf(i))*fac
180  550 CONTINUE
181  560 CONTINUE
182  nfe = nfe + mba
183 C Add identity matrix. -------------------------------------------------
184  570 ii = mband + 2
185  DO 580 i = 1,n
186  wm(ii) = wm(ii) + 1.0e0
187  580 ii = ii + meband
188 C Do LU decomposition of P. --------------------------------------------
189  CALL sgbtrf( n, n, ml, mu, wm(3), meband, iwm(21), ier)
190  IF (ier .NE. 0) ierpj = 1
191  RETURN
192 C----------------------- END OF SUBROUTINE SPREPJ ----------------------
193  END
std::string dimension(void) const
F77_RET_T const double const double * f
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:233
subroutine sprepj(NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC)
Definition: sprepj.f:1
T abs(T x)
Definition: pr-output.cc:3062
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:210