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
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 iownd, iowns,
7  1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
8  2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
9  INTEGER i, i1, i2, ier, ii, j, j1, jj, lenp,
10  1 mba, mband, meb1, meband, ml, ml3, mu, np1
11  DOUBLE PRECISION y, yh, ewt, ftem, savf, wm
12  DOUBLE PRECISION rowns,
13  1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
14  DOUBLE PRECISION con, di, fac, hl0, r, r0, srur, yi, yj, yjj,
15  1 vnorm
16  dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
17  1 wm(*), iwm(*)
18  COMMON /ls0001/ rowns(209),
19  2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
20  3 iownd(14), iowns(6),
21  4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
22  5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
23 C-----------------------------------------------------------------------
24 C PREPJ IS CALLED BY STODE TO COMPUTE AND PROCESS THE MATRIX
25 C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
26 C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF
27 C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
28 C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
29 C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN
30 C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
31 C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
32 C BY DGETRF IF MITER = 1 OR 2, AND BY DGBTRF IF MITER = 4 OR 5.
33 C
34 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
35 C WITH PREPJ USES THE FOLLOWING..
36 C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
37 C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STODE).
38 C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
39 C WM = REAL WORK SPACE FOR MATRICES. ON OUTPUT IT CONTAINS THE
40 C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
41 C OF P IF MITER IS 1, 2 , 4, OR 5.
42 C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
43 C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
44 C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
45 C WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
46 C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
47 C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND
48 C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
49 C EL0 = EL(1) (INPUT).
50 C IERPJ = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .GT. 0 IF
51 C P MATRIX FOUND TO BE SINGULAR.
52 C JCUR = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX
53 C (OR APPROXIMATION) IS NOW CURRENT.
54 C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
55 C MITER, N, NFE, AND NJE.
56 C-----------------------------------------------------------------------
57  nje = nje + 1
58  ierpj = 0
59  jcur = 1
60  hl0 = h*el0
61  go to(100, 200, 300, 400, 500), miter
62 C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
63  100 lenp = n*n
64  DO 110 i = 1,lenp
65  110 wm(i+2) = 0.0d0
66  CALL jac(neq, tn, y, 0, 0, wm(3), n)
67  con = -hl0
68  DO 120 i = 1,lenp
69  120 wm(i+2) = wm(i+2)*con
70  go to 240
71 C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
72  200 fac = vnorm(n, savf, ewt)
73  r0 = 1000.0d0*dabs(h)*uround*dble(n)*fac
74  IF (r0 .EQ. 0.0d0) r0 = 1.0d0
75  srur = wm(1)
76  j1 = 2
77  DO 230 j = 1,n
78  yj = y(j)
79  r = dmax1(srur*dabs(yj),r0/ewt(j))
80  y(j) = y(j) + r
81  fac = -hl0/r
82  ierr = 0
83  CALL f(neq, tn, y, ftem, ierr)
84  IF (ierr .LT. 0) RETURN
85  DO 220 i = 1,n
86  220 wm(i+j1) = (ftem(i) - savf(i))*fac
87  y(j) = yj
88  j1 = j1 + n
89  230 CONTINUE
90  nfe = nfe + n
91 C ADD IDENTITY MATRIX. -------------------------------------------------
92  240 j = 3
93  np1 = n + 1
94  DO 250 i = 1,n
95  wm(j) = wm(j) + 1.0d0
96  250 j = j + np1
97 C DO LU DECOMPOSITION ON P. --------------------------------------------
98  CALL dgetrf( n, n, wm(3), n, iwm(21), ier)
99  IF (ier .NE. 0) ierpj = 1
100  RETURN
101 C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
102  300 wm(2) = hl0
103  r = el0*0.1d0
104  DO 310 i = 1,n
105  310 y(i) = y(i) + r*(h*savf(i) - yh(i,2))
106  ierr = 0
107  CALL f(neq, tn, y, wm(3), ierr)
108  IF (ierr .LT. 0) RETURN
109  nfe = nfe + 1
110  DO 320 i = 1,n
111  r0 = h*savf(i) - yh(i,2)
112  di = 0.1d0*r0 - h*(wm(i+2) - savf(i))
113  wm(i+2) = 1.0d0
114  IF (dabs(r0) .LT. uround/ewt(i)) go to 320
115  IF (dabs(di) .EQ. 0.0d0) go to 330
116  wm(i+2) = 0.1d0*r0/di
117  320 CONTINUE
118  RETURN
119  330 ierpj = 1
120  RETURN
121 C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
122  400 ml = iwm(1)
123  mu = iwm(2)
124  ml3 = ml + 3
125  mband = ml + mu + 1
126  meband = mband + ml
127  lenp = meband*n
128  DO 410 i = 1,lenp
129  410 wm(i+2) = 0.0d0
130  CALL jac(neq, tn, y, ml, mu, wm(ml3), meband)
131  con = -hl0
132  DO 420 i = 1,lenp
133  420 wm(i+2) = wm(i+2)*con
134  go to 570
135 C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
136  500 ml = iwm(1)
137  mu = iwm(2)
138  mband = ml + mu + 1
139  mba = min0(mband,n)
140  meband = mband + ml
141  meb1 = meband - 1
142  srur = wm(1)
143  fac = vnorm(n, savf, ewt)
144  r0 = 1000.0d0*dabs(h)*uround*dble(n)*fac
145  IF (r0 .EQ. 0.0d0) r0 = 1.0d0
146  DO 560 j = 1,mba
147  DO 530 i = j,n,mband
148  yi = y(i)
149  r = dmax1(srur*dabs(yi),r0/ewt(i))
150  530 y(i) = y(i) + r
151  ierr = 0
152  CALL f(neq, tn, y, ftem, ierr)
153  IF (ierr .LT. 0) RETURN
154  DO 550 jj = j,n,mband
155  y(jj) = yh(jj,1)
156  yjj = y(jj)
157  r = dmax1(srur*dabs(yjj),r0/ewt(jj))
158  fac = -hl0/r
159  i1 = max0(jj-mu,1)
160  i2 = min0(jj+ml,n)
161  ii = jj*meb1 - ml + 2
162  DO 540 i = i1,i2
163  540 wm(ii+i) = (ftem(i) - savf(i))*fac
164  550 CONTINUE
165  560 CONTINUE
166  nfe = nfe + mba
167 C ADD IDENTITY MATRIX. -------------------------------------------------
168  570 ii = mband + 2
169  DO 580 i = 1,n
170  wm(ii) = wm(ii) + 1.0d0
171  580 ii = ii + meband
172 C DO LU DECOMPOSITION OF P. --------------------------------------------
173  CALL dgbtrf( n, n, ml, mu, wm(3), meband, iwm(21), ier)
174  IF (ier .NE. 0) ierpj = 1
175  RETURN
176 C----------------------- END OF SUBROUTINE PREPJ -----------------------
177  END