subroutine epsode (diffun, pederv, n, t0, h0, y0, tout, eps,
1 ierror, mf, index)
c this is the june 24, 1975 version of
c episode.. experimental package for integration of
c systems of ordinary differential equations,
c dy/dt = f(y,t), y = (y(1),y(2),...,y(n)) transpose,
c given the initial value of y.
c this code is for the ibm 370/195 at argonne national laboratory
c and is a modification of earlier versions by g.d.byrne
c and a.c.hindmarsh.
c references
c 1. g. d. byrne and a. c. hindmarsh, a polyalgorithm for the
c numerical solution of ordinary differential equations,
c ucrl-75652, lawrence livermore laboratory, p. o. box 808,
c livermore, ca 94550, april 1974. also in acm transactions
c on mathematical software, 1 (1975), pp. 71-96.
c
c 2. a. c. hindmarsh and g. d. byrne, episode.. an experimental
c package for the integration of systems of ordinary
c differential equations, ucid-30112, l.l.l., may, 1975.
c
c 3. a. c. hindmarsh, gear.. ordinary differential equation
c system solver, ucid-30001, rev. 3, l.l.l., december, 1974.
c
c-----------------------------------------------------------------------
c epsode is a driver subroutine for the episode package.
c epsode is to be called once for each output value of t.
c it then makes repeated calls to the core integrator
c subroutine, tstep.
c
c the input parameters are as follows.
c diffun= name of double precision subroutine diffun(n,t,y,ydot)
c which the user must supply and declare external. see
c note below.
c pederv= name of double precision subroutine pederv(n,t,y,pd,n0)
c which the user must supply and declare external. see
c note below.
c n = the number of differential equations (used only on
c first call, unless index = -1). n must never be
c increased during a given problem.
c t0 = the initial value of t, the independent variable
c (used for input only on first call).
c h0 = the step size in t (used for input only on the
c first call, unless index = 3 on input). when
c index = 3, h0 is the maximum absolute value of
c the step size to be used.
c y0 = a vector of length n containing the initial values of
c y (used for input only on first call).
c tout = the value of t at which output is desired next.
c integration will normally go beyond tout and
c interpolate to t = tout. (used only for input.)
c eps = the relative error bound (used only on first call,
c unless index = -1). this bound is used as follows.
c let r(i) denote the estimated relative local error
c in y(i), i.e. the error relative to ymax(i), as
c measured per step (of size h) or per ss units of t.
c then eps is a bound on the root-mean-square norm
c of the vector r, i.e.
c n
c sqrt ( sum ( r(i)**2 )/n ) .lt. eps.
c i=1
c the vector ymax is computed in epsode as described
c under ierror below.
c if error control per ss units of t is desired, set ss
c to a positive number after statement 10 (where it is
c now set to zero) and update it after statement 60.
c see also the comments on ss and ymax below.
c ierror = the error flag with values and meanings as follow.
c 1 absolute error is controlled. ymax(i) = 1.0.
c 2 error relative to abs(y) is controlled. if y(i) = 0.0
c a divide error will occur. ymax(i) = abs(y(i)).
c 3 error relative to the largest value of abs(y(i)) seen
c so far is controlled. if the initial value of y(i) is
c 0.0, then ymax(i) is set to 1.0 initially and remains
c at least 1.0.
c mf = the method flag (used only on first call, unless
c index = -1). allowed values are 10, 11, 12, 13,
c 20, 21, 22, 23. mf is an integer with two decimal
c digits, meth and miter (mf = 10*meth + miter). (mf
c can be thought of as the ordered pair (meth,miter).)
c meth is the basic method indicator.
c meth = 1 indicates variable-step size, variable-
c order adams method, suitable for non-
c stiff problems.
c meth = 2 indicates variable-step size, variable-
c order backward differentiation method,
c suitable for stiff problems.
c miter indicates the method of iterative correction
c (nonlinear system solution).
c miter = 0 indicates functional iteration (no
c partial derivatives needed).
c miter = 1 indicates a chord or semi-stationary
c newton method with closed form (exact)
c jacobian, which is computed in the
c user supplied subroutine
c pederv(n,t,y,pd,n0) described below.
c miter = 2 indicates a chord or semi-stationary
c newton method with an internally
c computed finite difference approximation
c to the jacobian.
c miter = 3 indicates a chord or semi-stationary
c newton method with an internally
c computed diagonal matrix approximation
c to the jacobian, based on a directional
c derivative.
c index = integer used on input to indicate type of call,
c with the following values and meanings..
c 1 this is the first call for this problem.
c 0 this is not the first call for this problem,
c and integration is to continue.
c -1 this is not the first call for the problem,
c and the user has reset n, eps, and/or mf.
c 2 same as 0 except that tout is to be hit
c exactly (no interpolation is done).
c assumes tout .ge. the current t.
c 3 same as 0 except control returns to calling
c program after one step. tout is ignored.
c since the normal output value of index is 0,
c it need not be reset for normal continuation.
c
c after the initial call, if a normal return occurred and a normal
c continuation is desired, simply reset tout and call again.
c all other parameters will be ready for the next call.
c a change of parameters with index = -1 can be made after
c either a successful or an unsuccessful return.
c
c the output parameters are as follows.
c t0 = the output value of t. if integration was successful,
c t0 = tout. otherwise, t0 is the last value of t
c reached successfully.
c h0 = the step size h used last, whether successfully or not.
c y0 = the computed values of y at t = t0.
c index = integer used on output to indicate results,
c with the following values and meanings..
c 0 integration was completed to tout or beyond.
c -1 the integration was halted after failing to pass the
c error test even after reducing h by a factor of
c 1.e10 from its initial value.
c -2 after some initial success, the integration was
c halted either by repeated error test failures or
c by a test on eps. possibly too much accuracy has
c been requested, or a bad choice of mf was made.
c -3 the integration was halted after failing to achieve
c corrector convergence even after reducing h by a
c factor of 1.e10 from its initial value.
c -4 immediate halt because of illegal values of input
c parameters. see printed message.
c -5 index was -1 on input, but the desired changes of
c parameters were not implemented because tout
c was not beyond t. interpolation to t = tout was
c performed as on a normal return. to continue,
c simply call again with index = -1 and a new tout.
c -6 index was 2 on input, but tout was not beyond t.
c no action was taken.
c
c in addition to epsode, the following subroutines are used by and
c provided in this package:
c interp(tout,y,n0,y0) interpolates to give output values at
c t = tout by using data in the y array.
c tstep(diffun,pederv,y,n0) is the core integration subroutine, which
c integrates over a single step and does associated error
c control.
c coset sets coefficients for use in tstep.
c adjust(y,n0) adjusts the history array y on reduction of order.
c pset(diffun,pederv,y,n0,con,miter,ier) computes and processes the
c jacobian matrix, j = df/dy.
c dec(n,n0,a,ip,ier) performs the lu decomposition of a matrix.
c sol(n,n0,a,b,ip) solves a linear system a*x = b, after dec
c has been called for the matrix a.
c note: pset, dec, and sol are called if and only if miter = 1
c or miter = 2.
c
c the user must furnish the following double precision subroutines
c and declare them external in his calling program:
c diffun(n,t,y,ydot) computes the function ydot = f(y,t),
c the right hand side of the ordinary
c differential equation system, where y
c and ydot are vectors of length n.
c pederv(n,t,y,pd,n0) computes the n by n jacobian matrix of
c partial derivatives and stores it in pd as
c an n0 by n0 array. pd(i,j) is to be set
c to the partial derivative of ydot(i) with
c respect to y(j). pederv is called if and
c only if miter = 1. for other values of
c miter, pederv can be a dummy subroutine.
c
c caution: at the present time the maximum number of differential
c equations, which can be solved by episode, is 20. to
c change this number to a new value, say nmax, change
c y(20,13) to y(nmax,13), ymax(20) to ymax(nmax),
c error(20) to error(nmax), save1(20) to save1(nmax),
c save2(20) to save2(nmax), pw(400) to pw(nmax*nmax),
c and ipiv(20) to ipiv(nmax) in the common and dimension
c statements below. also change the argument in the
c if...go to 440 statement (after the common statements)
c from 20 to nmax. no other changes need to be made to
c any other subroutine in this package when the maximum
c number of equations is changed. elsewhere, the column
c length of the y array is n0 instead of 20. the row
c length of y can be reduced from 13 to 6 if meth = 2.
c the array ipiv is used if and only if miter = 1 or
c miter = 2. the size of the pw array can be reduced
c to 1 if miter = 0 or to n if miter = 3.
c
c the common block epcom2 can be accessed externally by the user,
c if he desires to modify the error control behaviour of the package.
c the block contains ymax(20) used to control the error. see above
c and reference (2).
c
c the common block epcom9 can be accessed externally by the user,
c if he desires. it contains the step size last used successfully
c (hused), the order last used successfully (nqused), the
c number of steps taken so far (nstep), the number of function
c evaluations (diffun calls) so far (nfe), and the number of
c jacobian evaluations so far (nje).
c
c in a data statement below, lout is set to the logical unit number
c for the output of messages during integration. currently, lout
c = 6.
c-----------------------------------------------------------------------
c
external diffun, pederv
integer ierror, index, mf, n
integer ipiv, jstart, kflag, mfc, nc, nfe, nje,
1 nqused, nsq, nstep
integer i, kgo, nhcut, n0
integer lout
double precision eps, h0, tout, t0, y0
double precision epsc, epsj, error, hmax, h, hmin, hused,
1 pw, save1, save2, ss, t, uround, ymax
double precision ayi, d, epsj, h0, t0p, y
double precision hcut
double precision four, hundrd, one, ten, zero
dimension y(20,13)
dimension y0(n)
c
common /epcom1/t,h,hmin,hmax,epsc,ss,uround,nc,mfc,kflag,jstart
common /epcom2/ ymax(20)
common /epcom3/ error(20)
common /epcom4/ save1(20)
common /epcom5/ save2(20)
common /epcom6/ pw(400)
common /epcom7/ ipiv(20)
common /epcom8/ epsj,nsq
common /epcom9/ hused,nqused,nstep,nfe,nje
c
data lout /6/
data hcut /0.1d0/
data four /4.0d0/, hundrd /1.0d2/, one /1.0d0/,
1 ten /1.0d1/, zero /0.0d0/
if (index .eq. 0) go to 20
if (index .eq. 2) go to 25
if (index .eq. -1) go to 30
if (index .eq. 3) go to 40
if (index .ne. 1) go to 430
if (eps .le. zero) go to 400
if (n .le. 0) go to 410
if (n .gt. 20) go to 440
if ((t0-tout)*h0 .ge. zero) go to 420
c-----------------------------------------------------------------------
c if initial values for ymax other than those below are desired,
c they should be set here. all ymax(i) must be positive. if
c values for hmin or hmax, the bounds on the absolute value of h,
c other than those below, are desired, they also should be set here.
c if error per ss units of t is to be controlled, ss should be set
c to a positive value below. error per unit step is controlled
c when ss = 1. the default value for ss is 0 and yields control
c of error per step.
c-----------------------------------------------------------------------
c set uround, the machine roundoff constant, here.
c use statement below for short precision on ibm 360 or 370.
c uround = 9.53674e-7
c use statement below for single precision on cdc 7600 or 6600.
c uround = 7.105427406e-15
c use statement below for long precision on ibm 360 or 370.
c uround = 2.220446047d-16
uround = d1mach(4)
do 10 i = 1,n
go to (5, 6, 7), ierror
c ierror = 1, 2, 3 ------------------------------------------------
5 ymax(i) = one
go to 10
6 ymax(i) = dabs(y0(i))
go to 10
7 ymax(i) = dabs(y0(i))
if (ymax(i) .eq. zero) ymax(i) = one
10 y(i,1) = y0(i)
nc = n
t = t0
h = h0
if ((t+h) .eq. t) write(lout,15) t
15 format(/47h--- message from subroutine epsode in episode,,
1 24h the o.d.e. solver. ---/22h warning.. t + h = t =,
2 e18.8,18h in the next step./)
hmin = dabs(h0)
hmax = dabs(t0 - tout)*ten
epsc = eps
mfc = mf
jstart = 0
ss = zero
n0 = n
nsq = n0*n0
epsj = dsqrt(uround)
nhcut = 0
go to 50
c t0p is the previous output value of t0 for use in hmax. --------------
20 hmax = dabs(tout - t0p)*ten
go to 80
25 hmax = dabs(tout - t0p)*ten
c
if ((t-tout)*h .ge. zero) go to 460
go to 85
c
30 if ((t-tout)*h .ge. zero) go to 450
if (mf .ne. mfc) jstart = -1
nc = n
epsc = eps
mfc = mf
go to 45
c
40 hmax = h0
c
45 if ((t+h) .eq. t) write(lout,15) t
c
50 call tstep (diffun, pederv, y, n0)
c
kgo = 1 - kflag
go to (60, 100, 200, 300), kgo
c kflag = 0, -1, -2, -3 -----------------------------------------
c
60 continue
c-----------------------------------------------------------------------
c normal return from tstep.
c
c the weights ymax(i) are updated. if different values are desired,
c they should be set here. if ss is to be updated for control of
c error per ss units of t, it should also be done here. a test is
c made to determine if eps is too small for machine precision.
c
c any other tests or calculations that are required after each step
c should be inserted here.
c
c if index = 3, y0 is set to the current y values on return.
c if index = 2, h is controlled to hit tout (within roundoff
c error), and then the current y values are put in y0 on
c return. for any other value of index, control returns to
c the integrator unless tout has been reached. then
c interpolated values of y are computed and stored in y0 on
c return.
c if interpolation is not desired, the call to interp should
c be deleted and control transferred to statement 500 instead
c of 520.
c-----------------------------------------------------------------------
d = zero
do 70 i = 1,n
ayi = dabs(y(i,1))
go to (70, 66, 67), ierror
c ierror = 1, 2, 3 ---------------------------------------------
66 ymax(i) = ayi
go to 70
67 ymax(i) = dmax1(ymax(i), ayi)
70 d = d + (ayi/ymax(i))**2
d = d*(uround/eps)**2
if (d .gt. dfloat(n)) go to 250
if (index .eq. 3) go to 500
if (index .eq. 2) go to 85
80 if ((t-tout)*h .lt. zero) go to 45
call interp (tout, y, n0, y0)
t0 = tout
go to 520
85 if (((t+h)-tout)*h .le. zero) go to 45
if (dabs(t-tout) .le. hundrd*uround*hmax) go to 500
if ((t-tout)*h .ge. zero) go to 500
h = (tout - t)*(one - four*uround)
jstart = -1
go to 45
c-----------------------------------------------------------------------
c on an error return from tstep, an immediate return occurs if
c kflag = -2, and recovery attempts are made otherwise.
c to recover, h and hmin are reduced by a factor of .1 up to 10
c times before giving up.
c-----------------------------------------------------------------------
100 write (lout,101)
101 format(/47h--- message from subroutine epsode in episode,,
1 24h the o.d.e. solver. ---/)
write(lout,105) t,hmin
105 format(//35h kflag = -1 from integrator at t = ,e18.8/
1 40h error test failed with abs(h) = hmin =,e18.8/)
110 if (nhcut .eq. 10) go to 150
nhcut = nhcut + 1
hmin = hcut*hmin
h = hcut*h
write (lout,115) h
115 format(24h h has been reduced to ,e18.8,
1 26h and step will be retried//)
jstart = -1
go to 45
c
150 write (lout,155)
155 format(//44h problem appears unsolvable with given input//)
go to 500
c
200 write (lout,101)
write (lout,205) t,h,eps
205 format(//16h kflag = -2 t =,e18.8,4h h =,e18.8,6h eps =,e18.8/
1 50h the requested error is too small for integrator.//)
go to 500
c
250 write (lout,101)
write (lout,255) t,eps
255 format(//47h integration halted by subroutine epsode at t =,
1 e18.8/43h eps is too small for machine precision and/
2 29h problem being solved. eps =,e18.8//)
kflag = -2
go to 500
c
300 write (lout,101)
write (lout,305) t
305 format(//34h kflag = -3 from integrator at t =,e18.8/
1 45h corrector convergence could not be achieved/)
go to 110
c
400 write (lout,101)
write (lout,405) eps
405 format(//35h illegal input.. eps .le. 0. eps = ,e18.8//)
index = -4
return
c
410 write (lout,101)
write (lout,415) n
415 format(//31h illegal input.. n .le. 0. n = ,i8//)
index = -4
return
c
420 write (lout,101)
write (lout,425) t0,tout,h0
425 format(//39h illegal input.. (t0 - tout)*h0 .ge. 0./
1 5h t0 =,e18.8,7h tout =,e18.8,5h h0 =,e18.8//)
index = -4
return
c
430 write (lout,101)
write (lout,435) index
435 format(//24h illegal input.. index =,i8//)
index = -4
return
c
440 write (lout,101)
write (lout,445) n
445 format (//39h illegal input. the number of ordinary/
1 43h differential equations being solved is n =, i6/
2 43h storage allocation in subroutine epsode is/
3 47h too small. see writeups#epsode public......../)
index = -4
return
c
450 write (lout,101)
write (lout,455) t,tout,h
455 format(//46h index = -1 on input with (t - tout)*h .ge. 0./
1 44h interpolation was done as on normal return./
2 41h desired parameter changes were not made./
3 4h t =,e18.8,7h tout =,e18.8,4h h =,e18.8//)
call interp (tout, y, n0, y0)
t0 = tout
index = -5
return
c
460 write (lout,101)
write (lout,465) t,tout,h
465 format(//45h index = 2 on input with (t - tout)*h .ge. 0./
1 4h t =,e18.8,7h tout =,e18.8,4h h =,e18.8//)
index = -6
return
c
500 t0 = t
do 510 i = 1,n
510 y0(i) = y(i,1)
520 index = kflag
t0p = t0
h0 = hused
if (kflag .ne. 0) h0 = h
return
c----------------------- end of subroutine epsode ---------------------
end
subroutine interp (tout, y, n0, y0)
c-----------------------------------------------------------------------
c subroutine interp computes interpolated values of the dependent
c variable y and stores them in y0. the interpolation is to the
c point t = tout and uses the nordsieck history array y as follows..
c nq
c y0(i) = sum y(i,j+1)*s**j ,
c j=0
c where s = -(t-tout)/h.
c-----------------------------------------------------------------------
c caution: not all members of epcom1 are used in this subroutine.
c-----------------------------------------------------------------------
integer n0
integer jstart, kflag, mf, n
integer i, j, l
double precision tout, y, y0
double precision eps, h, hmax, hmin, ss, t, uround
double precision s, s1
double precision one
dimension y0(n0),y(n0,13)
c
common /epcom1/ t,h,hmin,hmax,eps,ss,uround,n,mf,kflag,jstart
data one /1.0d0/
do 10 i = 1,n
10 y0(i) = y(i,1)
l = jstart + 1
s = (tout - t)/h
s1 = one
do 30 j = 2,l
s1 = s1*s
do 20 i = 1,n
20 y0(i) = y0(i) + s1*y(i,j)
30 continue
return
c----------------------- end of subroutine interp --------------------
end
subroutine tstep (diffun, pederv, y, n0)
c-----------------------------------------------------------------------
c tstep performs one step of the integration of an initial value
c problem for a system of ordinary differential equations.
c communication with tstep is via the following variables..
c
c y an n0 by lmax array containing the dependent variables
c and their scaled derivatives. lmax is currently 6 for
c the variable step backward differentiation formulas,
c and 13 for the variable step adams formulas.
c (lmax -1) = maxder, the maximum order used.
c see subroutine coset. y(i,j+1) contains the
c j-th derivative of y(i), scaled by h**j/factorial(j)
c for j = 0,1,...,nq, where nq is the current order.
c n0 a constant integer .ge. n, used for dimensioning
c purposes.
c t the independent variable, updated on each step taken.
c h the step size to be attempted on the next step.
c h is altered by the error control algorithm during
c the solution of the problem. h can be either positive
c or negative, but its sign must remain constant
c throughout the problem run.
c hmin, the minimum and maximum absolute values of the step
c hmax size to be used for the step. these may be changed at
c any time, but the change will not take effect until the
c next change in h is made.
c eps the relative error bound. see description in
c subroutine epsode.
c ss the size of the time interval to be used for error
c control. a default value of 0 is used to produce
c control of error per step. see subroutine epsode.
c uround the unit of roundoff for the computer being used.
c n the number of first order ordinary differential
c equations being solved.
c mf the method flag. see description in subroutine epsode.
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 with abs(h) = hmin.
c -2 the requested error is smaller than can
c be handled for this problem.
c -3 corrector convergence could not be
c achieved for abs(h) = hmin.
c on a return with kflag negative, the values of t and
c the y array are as of the beginning of the last
c step and h is the last step size attempted.
c jstart an integer used on input and output.
c on input, it has the following values and meanings..
c 0 perform the first step.
c .gt.0 take a new step continuing from the last.
c .lt.0 take the next step with a new value of
c h and/or mf.
c on exit, jstart is set to nq, the current order of the
c method.
c ymax an array of n elements with which the estimated local
c errors in y are compared.
c error an array of n elements. error(i)/tq(2) is the
c estimated local error in y(i) per ss units of
c t or per step (of size h).
c save1, two arrays for working storage,
c save2 each of length n.
c pw a block of locations used for the partial derivatives
c of f with respect to y, if miter is not o. see
c description in subroutine epsode.
c ipiv an integer array of length n, which is used for pivot
c information for the linear algebraic system in the
c correction process, when miter = 1 or 2.
c
c the common block epcm10, declared below, is primarily intended
c for internal use, but it can be accessed externally.
c-----------------------------------------------------------------------
external pederv
integer n0
integer ipiv, jstart, kflag, l, lmax, meth, mf, n, nfe, nje,
1 nq, nqindx, nqused, nstep
integer i, iback, ier, iredo, j, j1, j2, m, mfold, mio,
1 miter, miter1, newj, nstepj
integer istepj, kfc, kfh, maxcor
double precision y
double precision el, eps, error, h, hmax, hmin, hused, pw,
1 save1, save2, ss, t, tau, tq, uround, ymax
double precision bnd, cnquot, con, conp, crate, d, drc,
1 d1, e, edn, eta, etamax, etamin, etaq, etaqm1,
2 etaqp1, eup, flotl, flotn, hold, hrl1, phrl1,
3 prl1, r, rc, rl1, r0, r1, told
double precision addon, bias1, bias2, bias3, crdown, delrc,
1 etacf, etamin, etamxf, etamx1, etamx2,
2 etamx3, onepsm, short, thresh
double precision one, pt5, zero
dimension y(n0,13)
c
common /epcom1/ t,h,hmin,hmax,eps,ss,uround,n,mf,kflag,jstart
common /epcom2/ ymax(1)
common /epcom3/ error(1)
common /epcom4/ save1(1)
common /epcom5/ save2(1)
common /epcom6/ pw(1)
common /epcom7/ ipiv(1)
common /epcom9/ hused,nqused,nstep,nfe,nje
common /epcm10/ tau(13),el(13),tq(5),lmax,meth,nq,l,nqindx
c
data istepj /20/, kfc /-3/, kfh /-7/, maxcor /3/
data addon /1.0d-6/, bias1 /2.5d1/, bias2 /2.5d1/,
1 bias3 /1.0d2/, crdown /0.1d0/, delrc /0.3d0/,
2 etacf /0.25d0/, etamin /0.1d0/, etamxf /0.2d0/,
3 etamx1 /1.0d4/, etamx2 /1.0d1/, etamx3 /1.5d0/,
4 onepsm /1.00001d0/, short /0.1d0/, thresh /1.3d0/
data one /1.0d0/, pt5 /0.5d0/, zero /0.0d0/
kflag = 0
told = t
flotn = dfloat(n)
if (jstart .gt. 0) go to 200
if (jstart .ne. 0) go to 150
c-----------------------------------------------------------------------
c on the first call, the order is set to 1 and the initial
c derivatives are calculated. etamax is the maximum ratio by
c which h can be increased in a single step. it is 1.e04 for the
c first step to compensate for the small initial h, then 10 for
c the next 10 steps, and then 1.5 thereafter. if a failure
c occurs (in corrector convergence or error test), etamax is set at 1
c for the next increase. etamin = .1 is the minimum ratio by which
c h can be reduced on any retry of a step.
c-----------------------------------------------------------------------
call diffun (n, t, y, save1)
do 110 i = 1,n
110 y(i,2) = h*save1(i)
meth = mf/10
miter = mf - 10*meth
miter1 = miter + 1
mfold = mf
nq = 1
l = 2
tau(1) = h
prl1 = one
rc = zero
etamax = etamx1
nqindx = 2
nstep = 0
nstepj = 0
nfe = 1
nje = 0
go to 200
c-----------------------------------------------------------------------
c if the user has changed h, then y must be rescaled. if the
c user has changed miter, then newj is set to miter to force
c the partial derivativees to be updated, if they are being used.
c-----------------------------------------------------------------------
150 if (mf .eq. mfold) go to 170
mio = miter
meth = mf/10
miter = mf - 10*meth
mfold = mf
if (miter .eq. mio) go to 170
newj = miter
miter1 = miter + 1
170 if (h .eq. hold) go to 200
eta = h/hold
h = hold
iredo = 3
go to 185
180 eta = dmax1(eta,hmin/dabs(h),etamin)
185 eta = dmin1(eta,hmax/dabs(h),etamax)
r1 = one
do 190 j = 2,l
r1 = r1*eta
do 190 i = 1,n
190 y(i,j) = y(i,j)*r1
h = h*eta
rc = rc*eta
if (iredo .eq. 0) go to 690
c-----------------------------------------------------------------------
c this section computes the predicted values by effectively
c multiplying the y array by the pascal triangle matrix. then
c coset is called to obtain el, the vector of coefficients of
c length nq + 1. rc is the ratio of new to old values of the
c coefficient h/el(2). when rc differs from 1 by more than
c delrc, newj is set to miter to force the partial derivatives
c to be updated, if used. delrc is 0.3. in any case, the partial
c derivatives are updated at least every 20-th step.
c-----------------------------------------------------------------------
200 t = t + h
do 210 j1 = 1,nq
do 210 j2 = j1,nq
j = (nq + j1) - j2
do 210 i = 1,n
210 y(i,j) = y(i,j) + y(i,j+1)
call coset
bnd = flotn*(tq(4)*eps)**2
rl1 = one/el(2)
rc = rc*(rl1/prl1)
prl1 = rl1
if (nstep .ge. nstepj+istepj) newj = miter
drc = dabs(rc-one)
if (drc .le. delrc) go to 215
newj = miter
crate = one
rc = one
go to 220
215 if ((miter .ne. 0) .and. (drc .ne. zero)) crate = one
c-----------------------------------------------------------------------
c up to 3 corrector iterations are taken. a convergence test is made
c on the root mean square norm of each correction, using bnd, which
c is dependent on eps. the sum of the corrections is accumulated in
c the vector error. the y array is not altered in the corrector
c loop. the updated y vector is stored temporarily in save1.
c-----------------------------------------------------------------------
220 do 230 i = 1,n
230 error(i) = zero
m = 0
call diffun (n, t, y, save2)
nfe = nfe + 1
if (newj .le. 0) go to 290
c-----------------------------------------------------------------------
c if indicated, the matrix p = i - h*rl1*j is reevaluated before
c starting the corrector iteration. newj is set to 0 as an
c indicator that this has been done. if miter = 1 or 2, p is
c computed and processed in pset. if miter = 3, the matrix is
c p = i - h*rl1*d, where d is a diagonal matrix. rl1 is 1/el(2).
c-----------------------------------------------------------------------
newj = 0
rc = one
nje = nje + 1
nstepj = nstep
go to (250, 240, 260), miter
240 nfe = nfe + n
250 con = -h*rl1
call pset(diffun, pederv, y, n0, con, miter, ier)
if (ier .ne. 0) go to 420
go to 350
260 r = rl1*short
do 270 i = 1,n
270 pw(i) = y(i,1) + r*(h*save2(i) - y(i,2))
call diffun(n, t, pw, save1)
nfe = nfe + 1
hrl1 = h*rl1
do 280 i = 1,n
r0 = h*save2(i) - y(i,2)
pw(i) = one
d = short*r0 - h*(save1(i) - save2(i))
save1(i) = zero
if (dabs(r0) .lt. uround*ymax(i)) go to 280
if (dabs(d) .eq. zero) go to 420
pw(i) = short*r0/d
save1(i) = pw(i)*rl1*r0
280 continue
go to 370
290 go to (295, 350, 350, 310), miter1
c-----------------------------------------------------------------------
c in the case of functional iteration, y is updated directly from
c the result of the last diffun call.
c-----------------------------------------------------------------------
295 d = zero
do 300 i = 1,n
r = rl1*(h*save2(i) - y(i,2))
d = d + ((r - error(i))/ymax(i))**2
save1(i) = y(i,1) + r
300 error(i) = r
go to 400
c-----------------------------------------------------------------------
c in the case of a chord method, the residual -g(y sub n(m))
c is computed and the linear system with that as right-hand side
c and p as coefficient matrix is solved, using the lu decomposition
c of p if miter = 1 or 2. if miter = 3 the scalar h*rl1 is updated.
c-----------------------------------------------------------------------
310 phrl1 = hrl1
hrl1 = h*rl1
if (hrl1 .eq. phrl1) go to 330
r = hrl1/phrl1
do 320 i = 1,n
d = one - r*(one - one/pw(i))
if (dabs(d) .eq. zero) go to 440
320 pw(i) = one/d
330 do 340 i = 1,n
340 save1(i) = pw(i)*(rl1*h*save2(i) - (rl1*y(i,2) + error(i)))
go to 370
350 do 360 i = 1,n
360 save1(i) = rl1*h*save2(i) - (rl1*y(i,2) + error(i))
call sol (n, n0, pw, save1, ipiv)
370 d = zero
do 380 i = 1,n
error(i) = error(i) + save1(i)
d = d + (save1(i)/ymax(i))**2
380 save1(i) = y(i,1) + error(i)
c-----------------------------------------------------------------------
c test for convergence. if m .gt. 0, an estimate of the square of
c the convergence rate constant is stored in crate, and this is used
c in the test.
c-----------------------------------------------------------------------
400 if (m .ne. 0) crate = dmax1(crdown*crate,d/d1)
if (d*dmin1(one,crate) .le. bnd) go to 450
d1 = d
m = m + 1
if (m .eq. maxcor) go to 410
call diffun (n, t, save1, save2)
go to (295, 350, 350, 310), miter1
c-----------------------------------------------------------------------
c the corrector iteration failed to converge in 3 tries. if partial
c derivatives are involved but are not up to date, they are
c reevaluated for the next try. otherwise the y array is restored
c to its values before prediction, and h is reduced,
c if possible. if not, a no-convergence exit is taken.
c-----------------------------------------------------------------------
410 nfe = nfe + maxcor - 1
if (newj .eq. -1) go to 440
420 t = told
etamax = one
do 430 j1 = 1,nq
do 430 j2 = j1,nq
j = (nq + j1) - j2
do 430 i = 1,n
430 y(i,j) = y(i,j) - y(i,j+1)
if (dabs(h) .le. hmin*onepsm) go to 680
eta = etacf
iredo = 1
go to 180
440 newj = miter
go to 220
c-----------------------------------------------------------------------
c the corrector has converged. newj is set to -1 if partial
c derivatives were used, to signal that they may need updating on
c subsequent steps. the error test is made and control passes to
c statement 500 if it fails.
c-----------------------------------------------------------------------
450 if (miter .ne. 0) newj = -1
nfe = nfe + m
d = zero
do 460 i = 1,n
460 d = d + (error(i)/ymax(i))**2
e = flotn*(tq(2)*eps)**2
if (d .gt. e) go to 500
c-----------------------------------------------------------------------
c after a successful step, the y array, tau, nstep, and nqindx are
c updated, and a new value of h at order nq is computed.
c the vector tau contains the nq + 1 most recent values of h.
c a change in nq up or down by 1 is considered if nqindx = 0.
c if nqindx = 1 and nq .lt. maxder, then error is saved
c for use in a possible order increase on the next step.
c a change in h or nq is made only of the increase in h
c is by a factor of at least 1.3.
c if not, nqindx is set to 2 to prevent testing for that many
c steps. if nq is changed, nqindx is set to nq + 1 (new value).
c-----------------------------------------------------------------------
kflag = 0
iredo = 0
nstep = nstep + 1
hused = h
nqused = nq
do 470 iback = 1,nq
i = l - iback
470 tau(i+1) = tau(i)
tau(1) = h
do 480 j = 1,l
do 480 i = 1,n
480 y(i,j) = y(i,j) + error(i)*el(j)
nqindx = nqindx - 1
if ((l .eq. lmax) .or. (nqindx .ne. 1)) go to 495
do 490 i = 1,n
490 y(i,lmax) = error(i)
conp = tq(5)
495 if (etamax .ne. one) go to 520
if (nqindx .lt. 2) nqindx = 2
go to 690
c-----------------------------------------------------------------------
c the error test failed. kflag keeps track of multiple failures.
c t and the y array are restored to their previous values. a new
c h for a retry of the step is computed. the order is kept fixed.
c-----------------------------------------------------------------------
500 kflag = kflag - 1
t = told
do 510 j1 = 1,nq
do 510 j2 = j1,nq
j = (nq + j1) - j2
do 510 i = 1,n
510 y(i,j) = y(i,j) - y(i,j+1)
newj = miter
etamax = one
if (dabs(h) .le. hmin*onepsm) go to 660
if (kflag .le. kfc) go to 630
iredo = 2
c compute ratio of new h to current h at the current order. ------------
520 flotl = dfloat(l)
etaq = one/((bias2*d/e)**(pt5/flotl) + addon)
if ((nqindx .ne. 0) .or. (kflag .ne. 0)) go to 580
etaqm1 = zero
if (nq .eq. 1) go to 540
c compute ratio of new h to current h at the current order less one. ---
d = zero
do 530 i = 1,n
530 d = d + (y(i,l)/ymax(i))**2
edn = flotn*(tq(1)*eps)**2
etaqm1 = one/((bias1*d/edn)**(pt5/(flotl - one)) + addon)
540 etaqp1 = zero
if (l .eq. lmax) go to 560
c compute ratio of new h to current h at current order plus one. -------
cnquot = (tq(5)/conp)*(h/tau(2))**l
d = zero
do 550 i = 1,n
550 d = d + ((error(i) - cnquot*y(i,lmax))/ymax(i))**2
eup = flotn*(tq(3)*eps)**2
etaqp1 = one/((bias3*d/eup)**(pt5/(flotl + one)) + addon)
560 nqindx = 2
if (etaq .ge. etaqp1) go to 570
if (etaqp1 .gt. etaqm1) go to 600
go to 590
570 if (etaq .lt. etaqm1) go to 590
580 if ((etaq .lt. thresh) .and. (kflag .eq. 0)) go to 690
eta = etaq
if ((kflag .le. -2) .and. (eta .gt. etamxf)) eta = etamxf
go to 180
590 if (etaqm1 .lt. thresh) go to 690
call adjust (y, n0)
l = nq
nq = nq - 1
eta = etaqm1
nqindx = l
go to 180
600 if (etaqp1 .lt. thresh) go to 690
nq = l
eta = etaqp1
l = l + 1
do 610 i = 1,n
610 y(i,l) = zero
nqindx = l
go to 180
c-----------------------------------------------------------------------
c control reaches this section if 3 or more consecutive failures
c have occurred. it is assumed that the elements of the y array
c have accumulated errors of the wrong order. the order is reduced
c by one, if possible. then h is reduced by a factor of 0.1 and
c the step is retried. after a total of 7 consecutive failures,
c an exit is taken with kflag = -2.
c-----------------------------------------------------------------------
630 if (kflag .eq. kfh) go to 670
if (nq .eq. 1) go to 640
eta = etamin
call adjust (y, n0)
l = nq
nq = nq - 1
nqindx = l
go to 180
640 eta = dmax1(etamin,hmin/dabs(h))
h = h*eta
call diffun (n, t, y, save1)
nfe = nfe + 1
do 650 i = 1,n
650 y(i,2) = h*save1(i)
nqindx = 10
go to 200
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 700
670 kflag = -2
go to 700
680 kflag = -3
go to 700
690 etamax = etamx3
if (nstep .le. 10) etamax = etamx2
700 hold = h
jstart = nq
return
c----------------------- end of subroutine tstep ---------------------
end
subroutine coset
c-----------------------------------------------------------------------
c coset is called by tstep and sets coefficients for use there.
c
c for each order nq, the coefficients in el are calculated by use of
c the generating polynomial lambda(x), with coefficients el(i):
c lambda(x) = el(1) + el(2)*x + ... + el(nq+1)*(x**nq).
c for the backward differentiation formulas,
c nq
c lambda(x) = product (1 + x/xi(i) ) .
c i = 1
c for the adams formulas,
c nq-1
c (d/dx) lambda(x) = c * product (1 + x/xi(i) ) ,
c i = 1
c lambda(-1) = 0, lambda(0) = 1,
c where c is a normalization constant.
c in both cases, xi(i) is defined by
c h*xi(i) = t sub n - t sub (n-i)
c = h + tau(1) + tau(2) + ... tau(i-1).
c
c coset also sets maxder, the maximum order of the formulas
c available. currently this is 5 for the backward differentiation
c formulas, and 12 for the adams formulas. to use different
c values (.le. 13), change the numbers in statements 1 and 2 below.
c
c in addition to variables described previously, communication
c with coset uses the following..
c tau = a vector of length 13 containing the past nq values
c of h.
c el = a vector of length 13 in which coset stores the
c coefficients for the corrector formula.
c tq = a vector of length 5 in which coset stores constants
c used for the convergence test, the error test, and
c selection of h at a new order.
c lmax = maxder + 1, where maxder is the maximum order
c available. lmax is the maximum number of columns
c of the y array to be used.
c meth = the basic method indicator.
c nq = the current order.
c l = nq + 1, the length of the vector stored in el, and
c the number of columns of the y array being used.
c nqindx = a counter controlling the frequency of order changes.
c an order change is about to be considered if
c nqindx = 1.
c-----------------------------------------------------------------------
c caution: not all members of epcom1 are used in this subroutine.
c-----------------------------------------------------------------------
integer jstart, kflag, l, lmax, meth, mf, n, nq, nqindx
integer i, iback, j, jp1, jstart, kflag, l, lmax, maxder,
1 meth, mf,n, nq, nqindx, nqm1
double precision el, eps, h, hmax, hmin, ss, t, tau, tq,
1 uround
double precision ahdss, cnqm1, csum, elp, em, em0, floti,
1 flotl, flotnq, hsum, hsum1, prod, rxi, s, xi
double precision cortes
double precision one, six, two, zero
dimension em(13)
c
common /epcom1/ t,h,hmin,hmax,eps,ss,uround,n,mf,kflag,jstart
common /epcm10/ tau(13),el(13),tq(5),lmax,meth,nq,l,nqindx
data cortes /0.1d0/
data one /1.0d0/, six /6.0d0/, two /2.0d0/, zero /0.0d0/
ahdss = one
if (ss .ne. zero) ahdss = dabs(h)/ss
flotl = dfloat(l)
nqm1 = nq - 1
go to (1, 2), meth
1 maxder = 12
go to 100
c
2 maxder = 5
go to 200
c
100 if (nq .ne. 1) go to 110
el(1) = one
el(2) = one
tq(1) = one
tq(2) = two*ahdss
tq(3) = six*tq(2)
tq(5) = one
go to 300
110 hsum = h
em(1) = one
flotnq = flotl - one
do 115 i = 2,l
115 em(i) = zero
do 150 j = 1,nqm1
if ((j .ne. nqm1) .or. (nqindx .ne. 1)) go to 130
s = one
csum = zero
do 120 i = 1,nqm1
csum = csum + s*em(i)/dfloat(i+1)
120 s = -s
tq(1) = ahdss*em(nqm1)/(flotnq*csum)
130 rxi = h/hsum
do 140 iback = 1,j
i = (j + 2) - iback
140 em(i) = em(i) + em(i-1)*rxi
150 hsum = hsum + tau(j)
c compute integral from -1 to 0 of polynomial and of x times it. -------
s = one
em0 = zero
csum = zero
do 160 i = 1,nq
floti = dfloat(i)
em0 = em0 + s*em(i)/floti
csum = csum + s*em(i)/(floti+1)
160 s = -s
c in el, form coefficients of normalized integrated polynomial. --------
s = one/em0
el(1) = one
do 170 i = 1,nq
170 el(i+1) = s*em(i)/dfloat(i)
xi = hsum/h
tq(2) = ahdss*xi*em0/csum
tq(5) = xi/el(l)
if (nqindx .ne. 1) go to 300
c for higher order control constant, multiply polynomial by 1+x/xi(q). -
rxi = one/xi
do 180 iback = 1,nq
i = (l + 1) - iback
180 em(i) = em(i) + em(i-1)*rxi
c compute integral of polynomial. --------------------------------------
s = one
csum = zero
do 190 i = 1,l
csum = csum + s*em(i)/dfloat(i+1)
190 s = -s
tq(3) = ahdss*flotl*em0/csum
go to 300
c
200 do 210 i = 3,l
210 el(i) = zero
el(1) = one
el(2) = one
hsum = h
hsum1 = zero
prod = one
rxi = one
if (nq .eq. 1) go to 240
do 230 j = 1,nqm1
c in el, construct coefficients of (1+x/xi(1))*...*(1+x/xi(j+1)). ------
hsum = hsum + tau(j)
hsum1 = hsum1 + tau(j)
prod = prod*(hsum/hsum1)
rxi = h/hsum
jp1 = j + 1
do 220 iback = 1,jp1
i = (j + 3) - iback
220 el(i) = el(i) + el(i-1)*rxi
230 continue
240 tq(2) = ahdss*el(2)*(one + prod)
tq(5) = (one + prod)/el(l)
if (nqindx .ne. 1) go to 300
cnqm1 = rxi/el(l)
elp = el(2) - rxi
tq(1) = ahdss*elp/cnqm1
hsum = hsum + tau(nq)
rxi = h/hsum
elp = el(2) + rxi
tq(3) = ahdss*elp*rxi*(one + prod)*(flotl + one)
300 tq(4) = cortes*tq(2)
lmax = maxder + 1
return
c----------------------- end of subroutine coset ---------------------
end
subroutine adjust (y, n0)
c-----------------------------------------------------------------------
c this subroutine adjusts the y array on reduction of order.
c see reference 1 for details.
c-----------------------------------------------------------------------
c caution: not all members of epcom1 are used in this subroutine.
c-----------------------------------------------------------------------
integer n0
integer jstart, kflag, l, lmax, meth, mf, n, nq, nqindx
integer i, iback, j, jp1, nqm1, nqm2
double precision y
double precision el, eps, h, hmax, hmin, ss, t, tau, tq, uround
double precision hsum, xi
double precision one, zero
dimension y(n0,13)
c
common /epcom1/ t,h,hmin,hmax,eps,ss,uround,n,mf,kflag,jstart
common /epcm10/ tau(13),el(13),tq(5),lmax,meth,nq,l,nqindx
data one /1.0d0/, zero /0.0d0/
if (nq .eq. 2) return
nqm1 = nq - 1
nqm2 = nq - 2
go to (100, 200), meth
c
100 do 110 j = 1,lmax
110 el(j) = zero
el(2) = one
hsum = zero
do 130 j = 1,nqm2
c construct coefficients of x*(x+xi(1))*...*(x+xi(j)). -----------------
hsum = hsum + tau(j)
xi = hsum/h
jp1 = j + 1
do 120 iback = 1,jp1
i = (j + 3) - iback
120 el(i) = el(i)*xi + el(i-1)
130 continue
c construct coefficients of integrated polynomial. ---------------------
do 140 j = 2,nqm1
140 el(j+1) = dfloat(nq)*el(j)/dfloat(j)
go to 300
c
200 do 210 j = 1,lmax
210 el(j) = zero
el(3) = one
hsum = zero
do 230 j = 1,nqm2
c construct coefficients of x*x*(x+xi(1))*...*(x+xi(j)). ---------------
hsum = hsum + tau(j)
xi = hsum/h
jp1 = j + 1
do 220 iback = 1,jp1
i = (j + 4) - iback
220 el(i) = el(i)*xi + el(i-1)
230 continue
c
c subtract correction terms from y array. ------------------------------
300 do 320 j = 3,nq
do 310 i = 1,n
310 y(i,j) = y(i,j) - y(i,l)*el(j)
320 continue
return
c----------------------- end of subroutine adjust --------------------
end
subroutine pset (diffun, pederv, y, n0, con, miter, ier)
c-----------------------------------------------------------------------
c pset is called by tstep to compute and to process the matrix
c p = i - (h/el(2))*j, where j is an approximation to the
c jacobian. j is computed by either the user supplied
c subroutine pederv, when miter = 1, or by finite differences,
c when miter = 2. j is stored in pw and replaced by p, using
c con = -h/el(2). then p is subjected to an lu decomposition
c for later solution of linear algebraic systems with p as the
c coefficient matrix.
c
c in addition to variables described previously, communication
c with pset uses the following..
c epsj = sqrt(uround), used in the numerical jacobian increments.
c nsq = n0**2.
c-----------------------------------------------------------------------
c caution: not all epcom1 variables are used inthis subroutine.
c-----------------------------------------------------------------------
integer ier, miter, n0
integer ipiv, jstart, kflag, mf, n, nsq
integer i, ier, j, j1, n
double precision con, y
double precision eps, epsj, h, hmax, hmin, pw, save1, save2,
1 ss, t, uround, ymax
double precision d, r, r0, t, yj
double precision one, rep, zero
dimension y(n0,1)
c
common /epcom1/ t,h,hmin,hmax,eps,ss,uround,n,mf,kflag,jstart
common /epcom2/ ymax(1)
common /epcom4/ save1(1)
common /epcom5/ save2(1)
common /epcom6/ pw(1)
common /epcom7/ ipiv(1)
common /epcom8/ epsj,nsq
data one /1.0d0/, rep /1.0d-3/, zero /0.0d0/
if (miter .eq. 2) go to 20
c if miter = 1, call pederv and multiply by a scalar. -----------------
call pederv (n, t, y, pw, n0)
do 10 i = 1,nsq
10 pw(i) = pw(i)*con
go to 60
c if miter = 2, make n calls to diffun to approximate j. ---------------
20 d = zero
do 30 i = 1,n
30 d = d + save2(i)**2
r0 = dabs(h)*dsqrt(d)*uround/rep
j1 = 0
do 50 j = 1,n
yj = y(j,1)
r = epsj*ymax(j)
r = dmax1(r,r0)
y(j,1) = y(j,1) + r
d = con/r
call diffun (n, t, y, save1)
do 40 i = 1,n
40 pw(i+j1) = (save1(i) - save2(i))*d
y(j,1) = yj
j1 = j1 + n0
50 continue
c add on the identity matrix. -----------------------------------------
60 j = 1
do 70 i = 1,n
pw(j) = pw(j) + one
70 j = j + (n0 + 1)
c get lu decomposition of p. -------------------------------------------
call dec (n, n0, pw, ipiv, ier)
return
c----------------------- end of subroutine pset -----------------------
end
subroutine dec (n, ndim, a, ip, ier)
c-----------------------------------------------------------------------
c matrix triangularization by gaussian elimination.
c input..
c n = order of matrix.
c ndim = declared dimension of array a .
c a = matrix to be triangularized.
c output..
c a(i,j), i.le.j = upper triangular factor, u .
c a(i,j), i.gt.j = multipliers = lower triangular factor, i - l.
c ip(k), k.lt.n = index of k-th pivot row.
c ip(n) = (-1)**(number of interchanges) or o .
c ier = 0 if a nonsingular, or k if a found to be
c singular at stage k.
c use sol to obtain solution of linear system.
c determ(a) = ip(n)*a(1,1)*a(2,2)*...*a(n,n).
c if ip(n)=0, a is singular, sol will divide by zero.
c interchanges finished in u , only partly in l .
c
c reference..
c c. b. moler, algorithm 423, linear equation solver,
c comm. assoc. comput. mach., 15 (1972), p. 274.
c-----------------------------------------------------------------------
integer ier, ip, n, ndim
integer i, j, k, kp1, m, nm1
double precision a
double precision t
double precision one, zero
dimension a(ndim,n),ip(n)
data one /1.0d0/, zero /0.0d0/
ier = 0
ip(n) = 1
if (n .eq. 1) go to 70
nm1 = n - 1
do 60 k = 1,nm1
kp1 = k + 1
m = k
do 10 i = kp1,n
10 if (dabs(a(i,k)) .gt. dabs(a(m,k))) m = i
ip(k) = m
t = a(m,k)
if (m .eq. k) go to 20
ip(n) = -ip(n)
a(m,k) = a(k,k)
a(k,k) = t
20 if (t .eq. zero) go to 80
t = one/t
do 30 i = kp1,n
30 a(i,k) = -a(i,k)*t
do 50 j = kp1,n
t = a(m,j)
a(m,j) = a(k,j)
a(k,j) = t
if (t .eq. zero) go to 50
do 40 i = kp1,n
40 a(i,j) = a(i,j) + a(i,k)*t
50 continue
60 continue
70 k = n
if (a(n,n) .eq. zero) go to 80
return
80 ier = k
ip(n) = 0
return
c----------------------- end of subroutine dec -----------------------
end
subroutine sol (n, ndim, a, b, ip)
c-----------------------------------------------------------------------
c solution of linear system, a*x = b .
c input..
c n = order of matrix.
c ndim = declared dimension of array a .
c a = triangularized matrix obtained from dec.
c b = right hand side vector.
c ip = pivot vector obtained from dec.
c do not use if dec has set ier .ne. 0.
c output..
c b = solution vector, x .
c-----------------------------------------------------------------------
integer ip, n, ndim
integer i, k, kb, km1, kp1, m, nm1
double precision a, b
double precision t
dimension a(ndim, n), b(n), ip(n)
c
if (n .eq. 1) go to 50
nm1 = n - 1
do 20 k = 1,nm1
kp1 = k + 1
m = ip(k)
t = b(m)
b(m) = b(k)
b(k) = t
do 10 i = kp1,n
10 b(i) = b(i) + a(i,k)*t
20 continue
do 40 kb = 1,nm1
km1 = n - kb
k = km1 + 1
b(k) = b(k)/a(k,k)
t = -b(k)
do 30 i = 1,km1
30 b(i) = b(i) + a(i,k)*t
40 continue
50 b(1) = b(1)/a(1,1)
return
c----------------------- end of subroutine sol -----------------------
end