V3FIT
prepj.f
1  SUBROUTINE prepj(neq, y, yh, nyh, ewt, ftem, savf, wm, iwm,
2  1 f, jac)
3  USE stel_kinds
4  USE liprec, ONLY: li_gbfa, li_gefa
5  IMPLICIT NONE
6 C-----------------------------------------------
7 C D u m m y A r g u m e n t s
8 C-----------------------------------------------
9  INTEGER nyh
10  INTEGER, DIMENSION(*) :: neq, iwm
11  REAL(rprec), DIMENSION(*) :: y
12  REAL(rprec), DIMENSION(nyh,*) :: yh
13  REAL(rprec), DIMENSION(*) :: ewt, ftem, savf, wm
14 C-----------------------------------------------
15 C C o m m o n B l o c k s
16 C-----------------------------------------------
17 C... /ls0001/
18  COMMON /ls0001/ rowns(209), ccmax, el0, h, hmin, hmxi, hu, rc, tn
19  1 , uround, iownd(14), iowns(6), icf, ierpj, iersl, jcur, jstart
20  2 , kflag, l, meth, miter, maxord, maxcor, msbp, mxncf, n, nq,
21  3 nst, nfe, nje, nqu
22  INTEGER iownd, iowns, icf, ierpj, iersl, jcur, jstart, kflag, l
23  1 , meth, miter, maxord, maxcor, msbp, mxncf, n, nq, nst, nfe,
24  2 nje, nqu
25  REAL(rprec) ::rowns, ccmax, el0, h, hmin, hmxi, hu, rc, tn,
26  1 uround
27 C-----------------------------------------------
28 C L o c a l V a r i a b l e s
29 C-----------------------------------------------
30  INTEGER :: i, i1, i2, ier, ii, j, j1, jj, lenp, mba, mband, meb1,
31  1 meband, ml, ml3, mu, np1
32  REAL(rprec) :: con, di, fac, hl0, r, r0, srur, yj, yjj
33 C-----------------------------------------------
34 C E x t e r n a l F u n c t i o n s
35 C-----------------------------------------------
36  REAL(rprec) , EXTERNAL :: vnorm
37  EXTERNAL f, jac
38 C-----------------------------------------------
39 clll. optimize
40 c-----------------------------------------------------------------------
41 c prepj is called by stode to compute and process the matrix
42 c p = i - h*el(1)*j , where j is an approximation to the jacobian.
43 c here j is computed by the user-supplied routine jac if
44 c miter = 1 or 4, or by finite differencing if miter = 2, 3, or 5.
45 c if miter = 3, a diagonal approximation to j is used.
46 c j is stored in wm and replaced by p. if miter .ne. 3, p is then
47 c subjected to lu decomposition in preparation for later solution
48 c of linear systems with p as coefficient matrix. this is done
49 c by sgefa if miter = 1 or 2, and by sgbfa if miter = 4 or 5.
50 c
51 c in addition to variables described previously, communication
52 c with prepj uses the following..
53 c y = array containing predicted values on entry.
54 c ftem = work array of length n (acor in stode).
55 c savf = array containing f evaluated at predicted y.
56 c wm = real work space for matrices. on output it contains the
57 c inverse diagonal matrix if miter = 3 and the lu decomposition
58 c of p if miter is 1, 2 , 4, or 5.
59 c storage of matrix elements starts at wm(3).
60 c wm also contains the following matrix-related data..
61 c wm(1) = sqrt(uround), used in numerical jacobian increments.
62 c wm(2) = h*el0, saved for later use if miter = 3.
63 c iwm = integer work space containing pivot information, starting at
64 c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band
65 c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
66 c el0 = el(1) (input).
67 c ierpj = output error flag, = 0 if no trouble, .gt. 0 if
68 c p matrix found to be singular.
69 c jcur = output flag = 1 to indicate that the jacobian matrix
70 c (or approximation) is now current.
71 c this routine also uses the common variables el0, h, tn, uround,
72 c miter, n, nfe, and nje.
73 c-----------------------------------------------------------------------
74  nje = nje + 1
75  ierpj = 0
76  jcur = 1
77  hl0 = h*el0
78  GOTO (100,200,300,400,500) miter
79 c IF miter = 1, CALL jac and multiply by scalar. -----------------------
80  100 CONTINUE
81  lenp = n*n
82  wm(3:lenp+2) = 0.0_dp
83  CALL jac (neq, tn, y, 0, 0, wm(3), n)
84  con = -hl0
85  wm(3:lenp+2) = wm(3:lenp+2)*con
86  GOTO 240
87 c IF miter = 2, make n calls to f to approximate j. --------------------
88  200 CONTINUE
89  fac = vnorm(n,savf,ewt)
90  r0 = 1000.0_dp*abs(h)*uround*n*fac
91  IF (r0 == 0.0_dp) r0 = 1.0_dp
92  srur = wm(1)
93  j1 = 2
94  DO j = 1, n
95  yj = y(j)
96  r = max(srur*abs(yj),r0/ewt(j))
97  y(j) = y(j) + r
98  fac = -hl0/r
99  CALL f (neq, tn, y, ftem)
100  wm(1+j1:n+j1) = (ftem(:n)-savf(:n))*fac
101  y(j) = yj
102  j1 = j1 + n
103  END DO
104  nfe = nfe + n
105 c add identity matrix. -------------------------------------------------
106  240 CONTINUE
107  j = 3
108  np1 = n + 1
109  wm(j:(n-1)*np1+j:np1) = wm(j:(n-1)*np1+j:np1) + 1.0_dp
110 c DO lu decomposition on p. --------------------------------------------
111  CALL li_gefa (wm(3:3+n*n), n, n, iwm(21:21+n), ier)
112  IF (ier /= 0) ierpj = 1
113  RETURN
114 c IF miter = 3, construct a diagonal approximation to j and p. ---------
115  300 CONTINUE
116  wm(2) = hl0
117  r = el0*0.1_dp
118  y(:n) = y(:n) + r*(h*savf(:n)-yh(:n,2))
119  CALL f (neq, tn, y, wm(3))
120  nfe = nfe + 1
121  DO i = 1, n
122  r0 = h*savf(i) - yh(i,2)
123  di = 0.1_dp*r0 - h*(wm(i+2)-savf(i))
124  wm(i+2) = 1.0_dp
125  IF (abs(r0) >= uround/ewt(i)) THEN
126  IF (abs(di) == 0.0_dp) GOTO 330
127  wm(i+2) = 0.1_dp*r0/di
128  END IF
129  END DO
130  RETURN
131  330 CONTINUE
132  ierpj = 1
133  RETURN
134 c IF miter = 4, CALL jac and multiply by scalar. -----------------------
135  400 CONTINUE
136  ml = iwm(1)
137  mu = iwm(2)
138  ml3 = ml + 3
139  mband = ml + mu + 1
140  meband = mband + ml
141  lenp = meband*n
142  wm(3:lenp+2) = 0.0_dp
143  CALL jac (neq, tn, y, ml, mu, wm(ml3), meband)
144  con = -hl0
145  wm(3:lenp+2) = wm(3:lenp+2)*con
146  GOTO 570
147 c IF miter = 5, make mband calls to f to approximate j. ----------------
148  500 CONTINUE
149  ml = iwm(1)
150  mu = iwm(2)
151  mband = ml + mu + 1
152  mba = min(mband,n)
153  meband = mband + ml
154  meb1 = meband - 1
155  srur = wm(1)
156  fac = vnorm(n,savf,ewt)
157  r0 = 1000.0_dp*abs(h)*uround*n*fac
158  IF (r0 == 0.0_dp) r0 = 1.0_dp
159  DO j = 1, mba
160  y(j:n:mband) = y(j:n:mband) + max(srur*abs(y(j:n:mband)),r0/
161  1 ewt(j:n:mband))
162  CALL f (neq, tn, y, ftem)
163  DO jj = j, n, mband
164  y(jj) = yh(jj,1)
165  yjj = y(jj)
166  r = max(srur*abs(yjj),r0/ewt(jj))
167  fac = -hl0/r
168  i1 = max(jj - mu,1)
169  i2 = min(jj + ml,n)
170  ii = jj*meb1 - ml + 2
171  wm(ii+i1:i2+ii) = (ftem(i1:i2)-savf(i1:i2))*fac
172  END DO
173  END DO
174  nfe = nfe + mba
175 c add identity matrix. -------------------------------------------------
176  570 CONTINUE
177  ii = mband + 2
178  wm(ii:(n-1)*meband+ii:meband) = wm(ii:(n-1)*meband+ii:meband) +
179  1 1.0_dp
180 c DO lu decomposition of p. --------------------------------------------
181  CALL li_gbfa (wm(3:3+meband*n), meband, n, ml, mu,
182  1 iwm(21:21+n), ier)
183  IF (ier /= 0) ierpj = 1
184 
185  END SUBROUTINE prepj
liprec::li_gefa
Definition: liprec.f:73
liprec::li_gbfa
Definition: liprec.f:5
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17