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 rownd, 1 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, vnorx dimension neq(1), y(1), yh(nyh,1), yh1(1), ewt(1), savf(1), 1 savr(1), acor(1), wm(1), iwm(1) common /ls0001/ rownd, 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 be set 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 lsodi. 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 = vnorx (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----------------------------------------------------------------------- ipup = 0 rc = 1.0d0 nslp = nst crate = 0.7d0 call pjac (neq, y, yh, nyh, ewt, acor, savr, savf, wm, iwm, 1 res, jac, adda ) 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 = vnorx (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) * vnorx (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) * vnorx (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 = vnorx (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