GNU Octave  3.8.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
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 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