subroutine stodi (neq, y, yh, nyh, yh1, ewt, savf, savr, 1 acor, wm, iwm, res, adda, jac, pjac, slvs ) clll. optimize external res, adda, jac, pjac, slvs integer neq, nyh, iwm integer iownd, ialth, ipup, lmax, meo, nqnyh, nslp, 1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 2 maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu integer i, i1, iredo, ires, iret, j, jb, kgo, m, ncf, newq double precision y, yh, yh1, ewt, savf, savr, acor, wm double precision conit, crate, el, elco, hold, rmax, tesco, 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround double precision dcon, ddn, del, delp, dsm, dup, 1 eljh, el1h, exdn, exsm, exup, 2 r, rh, rhdn, rhsm, rhup, told, vnorm dimension neq(1), y(1), yh(nyh,*), yh1(1), ewt(1), savf(1), 1 savr(1), acor(1), wm(*), iwm(*) common /ls0001/ conit, crate, el(13), elco(13,12), 1 hold, rmax, tesco(3,12), 2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround, iownd(14), 3 ialth, ipup, lmax, meo, nqnyh, nslp, 4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter, 5 maxord, maxcor, msbp, mxncf, n, nq, nst, nre, nje, nqu c----------------------------------------------------------------------- c stodi performs one step of the integration of an initial value c problem for a system of ordinary differential equations. c note.. stodi is independent of the value of the iteration method c indicator miter, and hence is independent c of the type of chord method used, or the jacobian structure. c communication with stodi is done with the following variables.. c c neq = integer array containing problem size in neq(1), and c passed as the neq argument in all calls to res, adda, c and jac. c y = an array of length .ge. n used as the y argument in c all calls to res, jac, and adda. c neq = integer array containing problem size in neq(1), and c passed as the neq argument in all calls to res, g, adda, c and jac c yh = an nyh by lmax array containing the dependent variables c and their approximate scaled derivatives, where c lmax = maxord + 1. yh(i,j+1) contains the approximate c j-th derivative of y(i), scaled by h**j/factorial(j) c (j = 0,1,...,nq). on entry for the first step, the first c two columns of yh must et from the initial values. c nyh = a constant integer .ge. n, the first dimension of yh. c yh1 = a one-dimensional array occupying the same space as yh. c ewt = an array of length n containing multiplicative weights c for local error measurements. local errors in y(i) are c compared to 1.0/ewt(i) in various error tests. c savf = an array of working storage, of length n. also used for c input of yh(*,maxord+2) when jstart = -1 and maxord is less c than the current order nq. c same as ydoti in driver. c savr = an array of working storage, of length n. c acor = a work array of length n used for the accumulated c corrections. on a succesful return, acor(i) contains c the estimated one-step local error in y(i). c wm,iwm = real and integer work arrays associated with matrix c operations in chord iteration. c pjac = name of routine to evaluate and preprocess jacobian matrix. c slvs = name of routine to solve linear system in chord iteration. c ccmax = maximum relative change in h*el0 before pjac is called. c h = the step size to be attempted on the next step. c h is altered by the error control algorithm during the c problem. h can be either positive or negative, but its c sign must remain constant throughout the problem. c hmin = the minimum absolute value of the step size h to be used. c hmxi = inverse of the maximum absolute value of h to be used. c hmxi = 0.0 is allowed and corresponds to an infinite hmax. c hmin and hmxi may be changed at any time, but will not c take effect until the next change of h is considered. c tn = the independent variable. tn is updated on each step taken. c jstart = an integer used for input only, with the following c values and meanings.. c 0 perform the first step. c .gt.0 take a new step continuing from the last. c -1 take the next step with a new value of h, maxord, c n, meth, miter, and/or matrix parameters. c -2 take the next step with a new value of h, c but with other inputs unchanged. c on return, jstart is set to 1 to facilitate continuation. c kflag = a completion code with the following meanings.. c 0 the step was succesful. c -1 the requested error could not be achieved. c -2 corrector convergence could not be achieved. c -3 res ordered immediate return. c -4 error condition from res could not be avoided. c -5 fatal error in pjac or slvs. c a return with kflag = -1, -2, or -4 means either c abs(h) = hmin or 10 consecutive failures occurred. c on a return with kflag negative, the values of tn and c the yh array are as of the beginning of the last c step, and h is the last step size attempted. c maxord = the maximum order of integration method to be allowed. c maxcor = the maximum number of corrector iterations allowed. c msbp = maximum number of steps between pjac calls. c mxncf = maximum number of convergence failures allowed. c meth/miter = the method flags. see description in driver. c n = the number of first-order differential equations. c----------------------------------------------------------------------- kflag = 0 told = tn ncf = 0 ierpj = 0 iersl = 0 jcur = 0 icf = 0 delp = 0.0d0 if (jstart .gt. 0) go to 200 if (jstart .eq. -1) go to 100 if (jstart .eq. -2) go to 160 c----------------------------------------------------------------------- c on the first call, the order is set to 1, and other variables are c initialized. rmax is the maximum ratio by which h can be increased c in a single step. it is initially 1.e4 to compensate for the small c initial h, but then is normally equal to 10. if a failure c occurs (in corrector convergence or error test), rmax is set at 2 c for the next increase. c----------------------------------------------------------------------- lmax = maxord + 1 nq = 1 l = 2 ialth = 2 rmax = 10000.0d0 rc = 0.0d0 el0 = 1.0d0 crate = 0.7d0 hold = h meo = meth nslp = 0 ipup = miter iret = 3 go to 140 c----------------------------------------------------------------------- c the following block handles preliminaries needed when jstart = -1. c ipup is set to miter to force a matrix update. c if an order increase is about to be considered (ialth = 1), c ialth is reset to 2 to postpone consideration one more step. c if the caller has changed meth, cfode is called to reset c the coefficients of the method. c if the caller has changed maxord to a value less than the current c order nq, nq is reduced to maxord, and a new h chosen accordingly. c if h is to be changed, yh must be rescaled. c if h or meth is being changed, ialth is reset to l = nq + 1 c to prevent further changes in h for that many steps. c----------------------------------------------------------------------- 100 ipup = miter lmax = maxord + 1 if (ialth .eq. 1) ialth = 2 if (meth .eq. meo) go to 110 call cfode (meth, elco, tesco) meo = meth if (nq .gt. maxord) go to 120 ialth = l iret = 1 go to 150 110 if (nq .le. maxord) go to 160 120 nq = maxord l = lmax do 125 i = 1,l 125 el(i) = elco(i,nq) nqnyh = nq*nyh rc = rc*el(1)/el0 el0 = el(1) conit = 0.5d0/dfloat(nq+2) ddn = vnorm (n, savf, ewt)/tesco(1,l) exdn = 1.0d0/dfloat(l) rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0) rh = dmin1(rhdn,1.0d0) iredo = 3 if (h .eq. hold) go to 170 rh = dmin1(rh,dabs(h/hold)) h = hold go to 175 c----------------------------------------------------------------------- c cfode is called to get all the integration coefficients for the c current meth. then the el vector and related constants are reset c whenever the order nq is changed, or at the start of the problem. c----------------------------------------------------------------------- 140 call cfode (meth, elco, tesco) 150 do 155 i = 1,l 155 el(i) = elco(i,nq) nqnyh = nq*nyh rc = rc*el(1)/el0 el0 = el(1) conit = 0.5d0/dfloat(nq+2) go to (160, 170, 200), iret c----------------------------------------------------------------------- c if h is being changed, the h ratio rh is checked against c rmax, hmin, and hmxi, and the yh array rescaled. ialth is set to c l = nq + 1 to prevent a change of h for that many steps, unless c forced by a convergence or error test failure. c----------------------------------------------------------------------- 160 if (h .eq. hold) go to 200 rh = h/hold h = hold iredo = 3 go to 175 170 rh = dmax1(rh,hmin/dabs(h)) 175 rh = dmin1(rh,rmax) rh = rh/dmax1(1.0d0,dabs(h)*hmxi*rh) r = 1.0d0 do 180 j = 2,l r = r*rh do 180 i = 1,n 180 yh(i,j) = yh(i,j)*r h = h*rh rc = rc*rh ialth = l if (iredo .eq. 0) go to 690 c----------------------------------------------------------------------- c this section computes the predicted values by effectively c multiplying the yh array by the pascal triangle matrix. c rc is the ratio of new to old values of the coefficient h*el(1). c when rc differs from 1 by more than ccmax, ipup is set to miter c to force pjac to be called. c in any case, pjac is called at least every msbp steps. c----------------------------------------------------------------------- 200 if (dabs(rc-1.0d0) .gt. ccmax) ipup = miter if (nst .ge. nslp+msbp) ipup = miter tn = tn + h i1 = nqnyh + 1 do 215 jb = 1,nq i1 = i1 - nyh cdir$ ivdep do 210 i = i1,nqnyh 210 yh1(i) = yh1(i) + yh1(i+nyh) 215 continue c----------------------------------------------------------------------- c up to maxcor corrector iterations are taken. a convergence test is c made on the r.m.s. norm of each correction, weighted by h and the c error weight vector ewt. the sum of the corrections is accumulated c in acor(i). the yh array is not altered in the corrector loop. c----------------------------------------------------------------------- 220 m = 0 do 230 i = 1,n savf(i) = yh(i,2) / h 230 y(i) = yh(i,1) if (ipup .le. 0) go to 240 c----------------------------------------------------------------------- c if indicated, the matrix p = a - h*el(1)*dr/dy is reevaluated and c preprocessed before starting the corrector iteration. ipup is set c to 0 as an indicator that this has been done. c----------------------------------------------------------------------- call pjac (neq, y, yh, nyh, ewt, acor, savr, savf, wm, iwm, 1 res, jac, adda ) ipup = 0 rc = 1.0d0 nslp = nst crate = 0.7d0 if (ierpj .eq. 0) go to 250 ires = ierpj go to (430, 435, 430), ires c get residual at predicted values, if not already done in pjac. ------- 240 ires = 1 call res ( neq, tn, y, savf, savr, ires ) nre = nre + 1 kgo = iabs(ires) go to ( 250, 435, 430 ) , kgo 250 do 260 i = 1,n 260 acor(i) = 0.0d0 c----------------------------------------------------------------------- c solve the linear system with the current residual as c right-hand side and p as coefficient matrix. c----------------------------------------------------------------------- 270 continue call slvs (wm, iwm, savr, savf) if (iersl .lt. 0) go to 430 if (iersl .gt. 0) go to 410 el1h = el(1) * h del = vnorm (n, savr, ewt) * dabs(h) do 380 i = 1,n acor(i) = acor(i) + savr(i) savf(i) = acor(i) + yh(i,2)/h 380 y(i) = yh(i,1) + el1h*acor(i) c----------------------------------------------------------------------- c test for convergence. if m.gt.0, an estimate of the convergence c rate constant is stored in crate, and this is used in the test. c----------------------------------------------------------------------- if (m .ne. 0) crate = dmax1(0.2d0*crate,del/delp) dcon = del*dmin1(1.0d0,1.5d0*crate)/(tesco(2,nq)*conit) if (dcon .le. 1.0d0) go to 460 m = m + 1 if (m .eq. maxcor) go to 410 if (m .ge. 2 .and. del .gt. 2.0d0*delp) go to 410 delp = del ires = 1 call res ( neq, tn, y, savf, savr, ires ) nre = nre + 1 kgo = iabs(ires) go to ( 270, 435, 410 ) , kgo c----------------------------------------------------------------------- c the correctors failed to converge, or res has returned abnormally. c on a convergence failure, if the jacobian is out of date, pjac is c called for the next try. otherwise the yh array is retracted to its c values before prediction, and h is reduced, if possible. c take an error exit if ires = 2, or h cannot be reduced, or mxncf c failures have occurred, or a fatal error occurred in pjac or slvs. c----------------------------------------------------------------------- 410 icf = 1 if (jcur .eq. 1) go to 430 ipup = miter go to 220 430 icf = 2 ncf = ncf + 1 rmax = 2.0d0 435 tn = told i1 = nqnyh + 1 do 445 jb = 1,nq i1 = i1 - nyh cdir$ ivdep do 440 i = i1,nqnyh 440 yh1(i) = yh1(i) - yh1(i+nyh) 445 continue if (ires .eq. 2) go to 680 if (ierpj .lt. 0 .or. iersl .lt. 0) go to 685 if (dabs(h) .le. hmin*1.00001d0) go to 450 if (ncf .eq. mxncf) go to 450 rh = 0.25d0 ipup = miter iredo = 1 go to 170 450 if (ires .eq. 3) go to 680 go to 670 c----------------------------------------------------------------------- c the corrector has converged. jcur is set to 0 c to signal that the jacobian involved may need updating later. c the local error test is made and control passes to statement 500 c if it fails. c----------------------------------------------------------------------- 460 jcur = 0 if (m .eq. 0) dsm = del/tesco(2,nq) if (m .gt. 0) dsm = dabs(h) * vnorm (n, acor, ewt)/tesco(2,nq) if (dsm .gt. 1.0d0) go to 500 c----------------------------------------------------------------------- c after a successful step, update the yh array. c consider changing h if ialth = 1. otherwise decrease ialth by 1. c if ialth is then 1 and nq .lt. maxord, then acor is saved for c use in a possible order increase on the next step. c if a change in h is considered, an increase or decrease in order c by one is considered also. a change in h is made only if it is by a c factor of at least 1.1. if not, ialth is set to 3 to prevent c testing for that many steps. c----------------------------------------------------------------------- kflag = 0 iredo = 0 nst = nst + 1 hu = h nqu = nq do 470 j = 1,l eljh = el(j)*h do 470 i = 1,n 470 yh(i,j) = yh(i,j) + eljh*acor(i) ialth = ialth - 1 if (ialth .eq. 0) go to 520 if (ialth .gt. 1) go to 700 if (l .eq. lmax) go to 700 do 490 i = 1,n 490 yh(i,lmax) = acor(i) go to 700 c----------------------------------------------------------------------- c the error test failed. kflag keeps track of multiple failures. c restore tn and the yh array to their previous values, and prepare c to try the step again. compute the optimum step size for this or c one lower order. after 2 or more failures, h is forced to decrease c by a factor of 0.1 or less. c----------------------------------------------------------------------- 500 kflag = kflag - 1 tn = told i1 = nqnyh + 1 do 515 jb = 1,nq i1 = i1 - nyh cdir$ ivdep do 510 i = i1,nqnyh 510 yh1(i) = yh1(i) - yh1(i+nyh) 515 continue rmax = 2.0d0 if (dabs(h) .le. hmin*1.00001d0) go to 660 if (kflag .le. -7) go to 660 iredo = 2 rhup = 0.0d0 go to 540 c----------------------------------------------------------------------- c regardless of the success or failure of the step, factors c rhdn, rhsm, and rhup are computed, by which h could be multiplied c at order nq - 1, order nq, or order nq + 1, respectively. c in the case of failure, rhup = 0.0 to avoid an order increase. c the largest of these is determined and the new order chosen c accordingly. if the order is to be increased, we compute one c additional scaled derivative. c----------------------------------------------------------------------- 520 rhup = 0.0d0 if (l .eq. lmax) go to 540 do 530 i = 1,n 530 savf(i) = acor(i) - yh(i,lmax) dup = dabs(h) * vnorm (n, savf, ewt)/tesco(3,nq) exup = 1.0d0/dfloat(l+1) rhup = 1.0d0/(1.4d0*dup**exup + 0.0000014d0) 540 exsm = 1.0d0/dfloat(l) rhsm = 1.0d0/(1.2d0*dsm**exsm + 0.0000012d0) rhdn = 0.0d0 if (nq .eq. 1) go to 560 ddn = vnorm (n, yh(1,l), ewt)/tesco(1,nq) exdn = 1.0d0/dfloat(nq) rhdn = 1.0d0/(1.3d0*ddn**exdn + 0.0000013d0) 560 if (rhsm .ge. rhup) go to 570 if (rhup .gt. rhdn) go to 590 go to 580 570 if (rhsm .lt. rhdn) go to 580 newq = nq rh = rhsm go to 620 580 newq = nq - 1 rh = rhdn if (kflag .lt. 0 .and. rh .gt. 1.0d0) rh = 1.0d0 go to 620 590 newq = l rh = rhup if (rh .lt. 1.1d0) go to 610 r = h*el(l)/dfloat(l) do 600 i = 1,n 600 yh(i,newq+1) = acor(i)*r go to 630 610 ialth = 3 go to 700 620 if ((kflag .eq. 0) .and. (rh .lt. 1.1d0)) go to 610 if (kflag .le. -2) rh = dmin1(rh,0.1d0) c----------------------------------------------------------------------- c if there is a change of order, reset nq, l, and the coefficients. c in any case h is reset according to rh and the yh array is rescaled. c then exit from 690 if the step was ok, or redo the step otherwise. c----------------------------------------------------------------------- if (newq .eq. nq) go to 170 630 nq = newq l = nq + 1 iret = 2 go to 150 c----------------------------------------------------------------------- c all returns are made through this section. h is saved in hold c to allow the caller to change h on the next step. c----------------------------------------------------------------------- 660 kflag = -1 go to 720 670 kflag = -2 go to 720 680 kflag = -1 - ires go to 720 685 kflag = -5 go to 720 690 rmax = 10.0d0 700 r = h/tesco(2,nqu) do 710 i = 1,n 710 acor(i) = acor(i)*r 720 hold = h jstart = 1 return c----------------------- end of subroutine stodi ----------------------- end