1 SUBROUTINE prepj (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
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,
16 dimension neq(*), y(*), yh(nyh,*), ewt(*), ftem(*), savf(*),
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
61 go to(100, 200, 300, 400, 500), miter
66 CALL jac(neq, tn, y, 0, 0, wm(3), n)
69 120 wm(i+2) = wm(i+2)*con
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
79 r = dmax1(srur*dabs(yj),r0/ewt(j))
83 CALL
f(neq, tn, y, ftem, ierr)
84 IF (ierr .LT. 0)
RETURN
86 220 wm(i+j1) = (ftem(i) - savf(i))*fac
98 CALL dgetrf( n, n, wm(3), n, iwm(21), ier)
99 IF (ier .NE. 0) ierpj = 1
105 310 y(i) = y(i) + r*(h*savf(i) - yh(i,2))
107 CALL
f(neq, tn, y, wm(3), ierr)
108 IF (ierr .LT. 0)
RETURN
111 r0 = h*savf(i) - yh(i,2)
112 di = 0.1d0*r0 - h*(wm(i+2) - savf(i))
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
130 CALL jac(neq, tn, y, ml, mu, wm(ml3), meband)
133 420 wm(i+2) = wm(i+2)*con
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
149 r = dmax1(srur*dabs(yi),r0/ewt(i))
152 CALL
f(neq, tn, y, ftem, ierr)
153 IF (ierr .LT. 0)
RETURN
154 DO 550 jj = j,n,mband
157 r = dmax1(srur*dabs(yjj),r0/ewt(jj))
161 ii = jj*meb1 - ml + 2
163 540 wm(ii+i) = (ftem(i) - savf(i))*fac
170 wm(ii) = wm(ii) + 1.0d0
173 CALL dgbtrf( n, n, ml, mu, wm(3), meband, iwm(21), ier)
174 IF (ier .NE. 0) ierpj = 1
subroutine prepj(NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC, IERR)
std::string dimension(void) const
F77_RET_T const double const double * f
double precision function vnorm(N, V, W)