N.I.M.R.O.D.  

dlsode.f90

Go to the documentation of this file.
00001 !------------------------------------------------------------------------------
00002 ! N.I.M.R.O.D. - Normal approximation Inference in Models with Random
00003 !                effects based on Ordinary Differential equations
00004 !------------------------------------------------------------------------------
00005 !
00006 ! VERSION : 1.0
00007 !
00008 ! MODULE: DLSODE
00009 !
00030       SUBROUTINE DLSODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, &
00031                        ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
00032       EXTERNAL F, JAC
00033       INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW, MF
00034       DOUBLE PRECISION Y, T, TOUT, RTOL, ATOL, RWORK
00035       DIMENSION NEQ(*), Y(*), RTOL(*), ATOL(*), RWORK(LRW), IWORK(LIW)
00036 !C**BEGIN PROLOGUE  DLSODE
00037 !C**PURPOSE  Livermore Solver for Ordinary Differential Equations.
00038 !C           DLSODE solves the initial-value problem for stiff or
00039 !C           nonstiff systems of first-order ODE's,
00040 !C              dy/dt = f(t,y),   or, in component form,
00041 !C              dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(N)),  i=1,...,N.
00042 !C**CATEGORY  I1A
00043 !C**TYPE      DOUBLE PRECISION (SLSODE-S, DLSODE-D)
00044 !C**KEYWORDS  ORDINARY DIFFERENTIAL EQUATIONS, INITIAL VALUE PROBLEM,
00045 !C            STIFF, NONSTIFF
00046 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
00047 !C            Center for Applied Scientific Computing, L-561
00048 !C            Lawrence Livermore National Laboratory
00049 !C            Livermore, CA 94551.
00050 !C**DESCRIPTION
00051 !C
00052 !C    NOTE: The "Usage" and "Arguments" sections treat only a subset of
00053 !C          available options, in condensed fashion.  The options
00054 !C          covered and the information supplied will support most
00055 !C          standard uses of DLSODE.
00056 !C
00057 !C          For more sophisticated uses, full details on all options are
00058 !C          given in the concluding section, headed "Long Description."
00059 !C          A synopsis of the DLSODE Long Description is provided at the
00060 !C          beginning of that section; general topics covered are:
00061 !C          - Elements of the call sequence; optional input and output
00062 !C          - Optional supplemental routines in the DLSODE package
00063 !C          - internal COMMON block
00064 !C
00065 !C*Usage:
00066 !C    Communication between the user and the DLSODE package, for normal
00067 !C    situations, is summarized here.  This summary describes a subset
00068 !C    of the available options.  See "Long Description" for complete
00069 !C    details, including optional communication, nonstandard options,
00070 !C    and instructions for special situations.
00071 !C
00072 !C    A sample program is given in the "Examples" section.
00073 !C
00074 !C    Refer to the argument descriptions for the definitions of the
00075 !C    quantities that appear in the following sample declarations.
00076 !C
00077 !C    For MF = 10,
00078 !C       PARAMETER  (LRW = 20 + 16*NEQ,           LIW = 20)
00079 !C    For MF = 21 or 22,
00080 !C       PARAMETER  (LRW = 22 +  9*NEQ + NEQ**2,  LIW = 20 + NEQ)
00081 !C    For MF = 24 or 25,
00082 !C       PARAMETER  (LRW = 22 + 10*NEQ + (2*ML+MU)*NEQ,
00083 !C      *                                         LIW = 20 + NEQ)
00084 !C
00085 !C       EXTERNAL F, JAC
00086 !C       INTEGER  NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK(LIW),
00087 !C      *         LIW, MF
00088 !C       DOUBLE PRECISION Y(NEQ), T, TOUT, RTOL, ATOL(ntol), RWORK(LRW)
00089 !C
00090 !C       CALL DLSODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
00091 !C      *            ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
00092 !C
00093 !C*Arguments:
00094 !C    F     :EXT    Name of subroutine for right-hand-side vector f.
00095 !C                  This name must be declared EXTERNAL in calling
00096 !C                  program.  The form of F must be:
00097 !C
00098 !C                  SUBROUTINE  F (NEQ, T, Y, YDOT)
00099 !C                  INTEGER  NEQ
00100 !C                  DOUBLE PRECISION  T, Y(*), YDOT(*)
00101 !C
00102 !C                  The inputs are NEQ, T, Y.  F is to set
00103 !C
00104 !C                  YDOT(i) = f(i,T,Y(1),Y(2),...,Y(NEQ)),
00105 !C                                                    i = 1, ..., NEQ .
00106 !C
00107 !C    NEQ   :IN     Number of first-order ODE's.
00108 !C
00109 !C    Y     :INOUT  Array of values of the y(t) vector, of length NEQ.
00110 !C                  Input:  For the first call, Y should contain the
00111 !C                          values of y(t) at t = T. (Y is an input
00112 !C                          variable only if ISTATE = 1.)
00113 !C                  Output: On return, Y will contain the values at the
00114 !C                          new t-value.
00115 !C
00116 !C    T     :INOUT  Value of the independent variable.  On return it
00117 !C                  will be the current value of t (normally TOUT).
00118 !C
00119 !C    TOUT  :IN     Next point where output is desired (.NE. T).
00120 !C
00121 !C    ITOL  :IN     1 or 2 according as ATOL (below) is a scalar or
00122 !C                  an array.
00123 !C
00124 !C    RTOL  :IN     Relative tolerance parameter (scalar).
00125 !C
00126 !C    ATOL  :IN     Absolute tolerance parameter (scalar or array).
00127 !C                  If ITOL = 1, ATOL need not be dimensioned.
00128 !C                  If ITOL = 2, ATOL must be dimensioned at least NEQ.
00129 !C
00130 !C                  The estimated local error in Y(i) will be controlled
00131 !C                  so as to be roughly less (in magnitude) than
00132 !C
00133 !C                  EWT(i) = RTOL*ABS(Y(i)) + ATOL     if ITOL = 1, or
00134 !C                  EWT(i) = RTOL*ABS(Y(i)) + ATOL(i)  if ITOL = 2.
00135 !C
00136 !C                  Thus the local error test passes if, in each
00137 !C                  component, either the absolute error is less than
00138 !C                  ATOL (or ATOL(i)), or the relative error is less
00139 !C                  than RTOL.
00140 !C
00141 !C                  Use RTOL = 0.0 for pure absolute error control, and
00142 !C                  use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative
00143 !C                  error control.  Caution:  Actual (global) errors may
00144 !C                  exceed these local tolerances, so choose them
00145 !C                  conservatively.
00146 !C
00147 !C    ITASK :IN     Flag indicating the task DLSODE is to perform.
00148 !C                  Use ITASK = 1 for normal computation of output
00149 !C                  values of y at t = TOUT.
00150 !C
00151 !C    ISTATE:INOUT  Index used for input and output to specify the state
00152 !C                  of the calculation.
00153 !C                  Input:
00154 !C                   1   This is the first call for a problem.
00155 !C                   2   This is a subsequent call.
00156 !C                  Output:
00157 !C                   1   Nothing was done, because TOUT was equal to T.
00158 !C                   2   DLSODE was successful (otherwise, negative).
00159 !C                       Note that ISTATE need not be modified after a
00160 !C                       successful return.
00161 !C                  -1   Excess work done on this call (perhaps wrong
00162 !C                       MF).
00163 !C                  -2   Excess accuracy requested (tolerances too
00164 !C                       small).
00165 !C                  -3   Illegal input detected (see printed message).
00166 !C                  -4   Repeated error test failures (check all
00167 !C                       inputs).
00168 !C                  -5   Repeated convergence failures (perhaps bad
00169 !C                       Jacobian supplied or wrong choice of MF or
00170 !C                       tolerances).
00171 !C                  -6   Error weight became zero during problem
00172 !C                       (solution component i vanished, and ATOL or
00173 !C                       ATOL(i) = 0.).
00174 !C
00175 !C    IOPT  :IN     Flag indicating whether optional inputs are used:
00176 !C                  0   No.
00177 !C                  1   Yes.  (See "Optional inputs" under "Long
00178 !C                      Description," Part 1.)
00179 !C
00180 !C    RWORK :WORK   Real work array of length at least:
00181 !C                  20 + 16*NEQ                    for MF = 10,
00182 !C                  22 +  9*NEQ + NEQ**2           for MF = 21 or 22,
00183 !C                  22 + 10*NEQ + (2*ML + MU)*NEQ  for MF = 24 or 25.
00184 !C
00185 !C    LRW   :IN     Declared length of RWORK (in user's DIMENSION
00186 !C                  statement).
00187 !C
00188 !C    IWORK :WORK   Integer work array of length at least:
00189 !C                  20        for MF = 10,
00190 !C                  20 + NEQ  for MF = 21, 22, 24, or 25.
00191 !C
00192 !C                  If MF = 24 or 25, input in IWORK(1),IWORK(2) the
00193 !C                  lower and upper Jacobian half-bandwidths ML,MU.
00194 !C
00195 !C                  On return, IWORK contains information that may be
00196 !C                  of interest to the user:
00197 !C
00198 !C           Name   Location   Meaning
00199 !C           -----  ---------  -----------------------------------------
00200 !C           NST    IWORK(11)  Number of steps taken for the problem so
00201 !C                             far.
00202 !C           NFE    IWORK(12)  Number of f evaluations for the problem
00203 !C                             so far.
00204 !C           NJE    IWORK(13)  Number of Jacobian evaluations (and of
00205 !C                             matrix LU decompositions) for the problem
00206 !C                             so far.
00207 !C           NQU    IWORK(14)  Method order last used (successfully).
00208 !C           LENRW  IWORK(17)  Length of RWORK actually required.  This
00209 !C                             is defined on normal returns and on an
00210 !C                             illegal input return for insufficient
00211 !C                             storage.
00212 !C           LENIW  IWORK(18)  Length of IWORK actually required.  This
00213 !C                             is defined on normal returns and on an
00214 !C                             illegal input return for insufficient
00215 !C                             storage.
00216 !C
00217 !C    LIW   :IN     Declared length of IWORK (in user's DIMENSION
00218 !C                  statement).
00219 !C
00220 !C    JAC   :EXT    Name of subroutine for Jacobian matrix (MF =
00221 !C                  21 or 24).  If used, this name must be declared
00222 !C                  EXTERNAL in calling program.  If not used, pass a
00223 !C                  dummy name.  The form of JAC must be:
00224 !C
00225 !C                  SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
00226 !C                  INTEGER  NEQ, ML, MU, NROWPD
00227 !C                  DOUBLE PRECISION  T, Y(*), PD(NROWPD,*)
00228 !C
00229 !C                  See item c, under "Description" below for more
00230 !C                  information about JAC.
00231 !C
00232 !C    MF    :IN     Method flag.  Standard values are:
00233 !C                  10  Nonstiff (Adams) method, no Jacobian used.
00234 !C                  21  Stiff (BDF) method, user-supplied full Jacobian.
00235 !C                  22  Stiff method, internally generated full
00236 !C                      Jacobian.
00237 !C                  24  Stiff method, user-supplied banded Jacobian.
00238 !C                  25  Stiff method, internally generated banded
00239 !C                      Jacobian.
00240 !C
00241 !C*Description:
00242 !C    DLSODE solves the initial value problem for stiff or nonstiff
00243 !C    systems of first-order ODE's,
00244 !C
00245 !C       dy/dt = f(t,y) ,
00246 !C
00247 !C    or, in component form,
00248 !C
00249 !C       dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ))
00250 !C                                                 (i = 1, ..., NEQ) .
00251 !C
00252 !C    DLSODE is a package based on the GEAR and GEARB packages, and on
00253 !C    the October 23, 1978, version of the tentative ODEPACK user
00254 !C    interface standard, with minor modifications.
00255 !C
00256 !C    The steps in solving such a problem are as follows.
00257 !C
00258 !C    a. First write a subroutine of the form
00259 !C
00260 !C          SUBROUTINE  F (NEQ, T, Y, YDOT)
00261 !C          INTEGER  NEQ
00262 !C          DOUBLE PRECISION  T, Y(*), YDOT(*)
00263 !C
00264 !C       which supplies the vector function f by loading YDOT(i) with
00265 !C       f(i).
00266 !C
00267 !C    b. Next determine (or guess) whether or not the problem is stiff.
00268 !C       Stiffness occurs when the Jacobian matrix df/dy has an
00269 !C       eigenvalue whose real part is negative and large in magnitude
00270 !C       compared to the reciprocal of the t span of interest.  If the
00271 !C       problem is nonstiff, use method flag MF = 10.  If it is stiff,
00272 !C       there are four standard choices for MF, and DLSODE requires the
00273 !C       Jacobian matrix in some form.  This matrix is regarded either
00274 !C       as full (MF = 21 or 22), or banded (MF = 24 or 25).  In the
00275 !C       banded case, DLSODE requires two half-bandwidth parameters ML
00276 !C       and MU. These are, respectively, the widths of the lower and
00277 !C       upper parts of the band, excluding the main diagonal.  Thus the
00278 !C       band consists of the locations (i,j) with
00279 !C
00280 !C          i - ML <= j <= i + MU ,
00281 !C
00282 !C       and the full bandwidth is ML + MU + 1 .
00283 !C
00284 !C    c. If the problem is stiff, you are encouraged to supply the
00285 !C       Jacobian directly (MF = 21 or 24), but if this is not feasible,
00286 !C       DLSODE will compute it internally by difference quotients (MF =
00287 !C       22 or 25).  If you are supplying the Jacobian, write a
00288 !C       subroutine of the form
00289 !C
00290 !C          SUBROUTINE  JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
00291 !C          INTEGER  NEQ, ML, MU, NRWOPD
00292 !C          DOUBLE PRECISION  T, Y(*), PD(NROWPD,*)
00293 !C
00294 !C       which provides df/dy by loading PD as follows:
00295 !C       - For a full Jacobian (MF = 21), load PD(i,j) with df(i)/dy(j),
00296 !C         the partial derivative of f(i) with respect to y(j).  (Ignore
00297 !C         the ML and MU arguments in this case.)
00298 !C       - For a banded Jacobian (MF = 24), load PD(i-j+MU+1,j) with
00299 !C         df(i)/dy(j); i.e., load the diagonal lines of df/dy into the
00300 !C         rows of PD from the top down.
00301 !C       - In either case, only nonzero elements need be loaded.
00302 !C
00303 !C    d. Write a main program that calls subroutine DLSODE once for each
00304 !C       point at which answers are desired.  This should also provide
00305 !C       for possible use of logical unit 6 for output of error messages
00306 !C       by DLSODE.
00307 !C
00308 !C       Before the first call to DLSODE, set ISTATE = 1, set Y and T to
00309 !C       the initial values, and set TOUT to the first output point.  To
00310 !C       continue the integration after a successful return, simply
00311 !C       reset TOUT and call DLSODE again.  No other parameters need be
00312 !C       reset.
00313 !C
00314 !C*Examples:
00315 !C    The following is a simple example problem, with the coding needed
00316 !C    for its solution by DLSODE. The problem is from chemical kinetics,
00317 !C    and consists of the following three rate equations:
00318 !C
00319 !C       dy1/dt = -.04*y1 + 1.E4*y2*y3
00320 !C       dy2/dt = .04*y1 - 1.E4*y2*y3 - 3.E7*y2**2
00321 !C       dy3/dt = 3.E7*y2**2
00322 !C
00323 !C    on the interval from t = 0.0 to t = 4.E10, with initial conditions
00324 !C    y1 = 1.0, y2 = y3 = 0. The problem is stiff.
00325 !C
00326 !C    The following coding solves this problem with DLSODE, using
00327 !C    MF = 21 and printing results at t = .4, 4., ..., 4.E10.  It uses
00328 !C    ITOL = 2 and ATOL much smaller for y2 than for y1 or y3 because y2
00329 !C    has much smaller values.  At the end of the run, statistical
00330 !C    quantities of interest are printed.
00331 !C
00332 !C       EXTERNAL  FEX, JEX
00333 !C       INTEGER  IOPT, IOUT, ISTATE, ITASK, ITOL, IWORK(23), LIW, LRW,
00334 !C      *         MF, NEQ
00335 !C       DOUBLE PRECISION  ATOL(3), RTOL, RWORK(58), T, TOUT, Y(3)
00336 !C       NEQ = 3
00337 !C       Y(1) = 1.D0
00338 !C       Y(2) = 0.D0
00339 !C       Y(3) = 0.D0
00340 !C       T = 0.D0
00341 !C       TOUT = .4D0
00342 !C       ITOL = 2
00343 !C       RTOL = 1.D-4
00344 !C       ATOL(1) = 1.D-6
00345 !C       ATOL(2) = 1.D-10
00346 !C       ATOL(3) = 1.D-6
00347 !C       ITASK = 1
00348 !C       ISTATE = 1
00349 !C       IOPT = 0
00350 !C       LRW = 58
00351 !C       LIW = 23
00352 !C       MF = 21
00353 !C       DO 40 IOUT = 1,12
00354 !C         CALL DLSODE (FEX, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
00355 !C      *               ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEX, MF)
00356 !         WRITE(6,20)  T, Y(1), Y(2), Y(3)
00357 !   20    FORMAT(' At t =',D12.4,'   y =',3D14.6)
00358 !         IF (ISTATE .LT. 0)  GO TO 80
00359 !   40    TOUT = TOUT*10.D0
00360 !       WRITE(6,60)  IWORK(11), IWORK(12), IWORK(13)
00361 !   60  FORMAT(/' No. steps =',i4,',  No. f-s =',i4,',  No. J-s =',i4)
00362 !C       STOP
00363 !C   80  WRITE(6,90)  ISTATE
00364 !C   90  FORMAT(///' Error halt.. ISTATE =',I3)
00365 !C       STOP
00366 !C       END
00367 !C
00368 !C       SUBROUTINE  FEX (NEQ, T, Y, YDOT)
00369 !C       INTEGER  NEQ
00370 !C       DOUBLE PRECISION  T, Y(3), YDOT(3)
00371 !C       YDOT(1) = -.04D0*Y(1) + 1.D4*Y(2)*Y(3)
00372 !C       YDOT(3) = 3.D7*Y(2)*Y(2)
00373 !C       YDOT(2) = -YDOT(1) - YDOT(3)
00374 !C       RETURN
00375 !C       END
00376 !C
00377 !C       SUBROUTINE  JEX (NEQ, T, Y, ML, MU, PD, NRPD)
00378 !C       INTEGER  NEQ, ML, MU, NRPD
00379 !C       DOUBLE PRECISION  T, Y(3), PD(NRPD,3)
00380 !C       PD(1,1) = -.04D0
00381 !C       PD(1,2) = 1.D4*Y(3)
00382 !C       PD(1,3) = 1.D4*Y(2)
00383 !C       PD(2,1) = .04D0
00384 !C       PD(2,3) = -PD(1,3)
00385 !C       PD(3,2) = 6.D7*Y(2)
00386 !C       PD(2,2) = -PD(1,2) - PD(3,2)
00387 !C       RETURN
00388 !C       END
00389 !C
00390 !C    The output from this program (on a Cray-1 in single precision)
00391 !C    is as follows.
00392 !C
00393 !C    At t =  4.0000e-01   y =  9.851726e-01  3.386406e-05  1.479357e-02
00394 !C    At t =  4.0000e+00   y =  9.055142e-01  2.240418e-05  9.446344e-02
00395 !C    At t =  4.0000e+01   y =  7.158050e-01  9.184616e-06  2.841858e-01
00396 !C    At t =  4.0000e+02   y =  4.504846e-01  3.222434e-06  5.495122e-01
00397 !C    At t =  4.0000e+03   y =  1.831701e-01  8.940379e-07  8.168290e-01
00398 !C    At t =  4.0000e+04   y =  3.897016e-02  1.621193e-07  9.610297e-01
00399 !C    At t =  4.0000e+05   y =  4.935213e-03  1.983756e-08  9.950648e-01
00400 !C    At t =  4.0000e+06   y =  5.159269e-04  2.064759e-09  9.994841e-01
00401 !C    At t =  4.0000e+07   y =  5.306413e-05  2.122677e-10  9.999469e-01
00402 !C    At t =  4.0000e+08   y =  5.494530e-06  2.197825e-11  9.999945e-01
00403 !C    At t =  4.0000e+09   y =  5.129458e-07  2.051784e-12  9.999995e-01
00404 !C    At t =  4.0000e+10   y = -7.170603e-08 -2.868241e-13  1.000000e+00
00405 !C
00406 !C    No. steps = 330,  No. f-s = 405,  No. J-s = 69
00407 !C
00408 !C*Accuracy:
00409 !C    The accuracy of the solution depends on the choice of tolerances
00410 !C    RTOL and ATOL.  Actual (global) errors may exceed these local
00411 !C    tolerances, so choose them conservatively.
00412 !C
00413 !C*Cautions:
00414 !C    The work arrays should not be altered between calls to DLSODE for
00415 !C    the same problem, except possibly for the conditional and optional
00416 !C    inputs.
00417 !C
00418 !C*Portability:
00419 !C    Since NEQ is dimensioned inside DLSODE, some compilers may object
00420 !C    to a call to DLSODE with NEQ a scalar variable.  In this event,
00421 !C    use DIMENSION NEQ(1).  Similar remarks apply to RTOL and ATOL.
00422 !C
00423 !C    Note to Cray users:
00424 !C    For maximum efficiency, use the CFT77 compiler.  Appropriate
00425 !C    compiler optimization directives have been inserted for CFT77.
00426 !C
00427 !C*Reference:
00428 !C    Alan C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
00429 !C    Solvers," in Scientific Computing, R. S. Stepleman, et al., Eds.
00430 !C    (North-Holland, Amsterdam, 1983), pp. 55-64.
00431 !C
00432 !C*Long Description:
00433 !C    The following complete description of the user interface to
00434 !C    DLSODE consists of four parts:
00435 !C
00436 !C    1.  The call sequence to subroutine DLSODE, which is a driver
00437 !C        routine for the solver.  This includes descriptions of both
00438 !C        the call sequence arguments and user-supplied routines.
00439 !C        Following these descriptions is a description of optional
00440 !C        inputs available through the call sequence, and then a
00441 !C        description of optional outputs in the work arrays.
00442 !C
00443 !C    2.  Descriptions of other routines in the DLSODE package that may
00444 !C        be (optionally) called by the user.  These provide the ability
00445 !C        to alter error message handling, save and restore the internal
00446 !C        COMMON, and obtain specified derivatives of the solution y(t).
00447 !C
00448 !C    3.  Descriptions of COMMON block to be declared in overlay or
00449 !C        similar environments, or to be saved when doing an interrupt
00450 !C        of the problem and continued solution later.
00451 !C
00452 !C    4.  Description of two routines in the DLSODE package, either of
00453 !C        which the user may replace with his own version, if desired.
00454 !C        These relate to the measurement of errors.
00455 !C
00456 !C
00457 !C                        Part 1.  Call Sequence
00458 !C                        ----------------------
00459 !C
00460 !C    Arguments
00461 !C    ---------
00462 !C    The call sequence parameters used for input only are
00463 !C
00464 !C       F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, MF,
00465 !C
00466 !C    and those used for both input and output are
00467 !C
00468 !C       Y, T, ISTATE.
00469 !C
00470 !C    The work arrays RWORK and IWORK are also used for conditional and
00471 !C    optional inputs and optional outputs.  (The term output here
00472 !C    refers to the return from subroutine DLSODE to the user's calling
00473 !C    program.)
00474 !C
00475 !C    The legality of input parameters will be thoroughly checked on the
00476 !C    initial call for the problem, but not checked thereafter unless a
00477 !C    change in input parameters is flagged by ISTATE = 3 on input.
00478 !C
00479 !C    The descriptions of the call arguments are as follows.
00480 !C
00481 !C    F        The name of the user-supplied subroutine defining the ODE
00482 !C             system.  The system must be put in the first-order form
00483 !C             dy/dt = f(t,y), where f is a vector-valued function of
00484 !C             the scalar t and the vector y. Subroutine F is to compute
00485 !C             the function f. It is to have the form
00486 !C
00487 !C                SUBROUTINE F (NEQ, T, Y, YDOT)
00488 !C                DOUBLE PRECISION  T, Y(*), YDOT(*)
00489 !C
00490 !C             where NEQ, T, and Y are input, and the array YDOT =
00491 !C             f(T,Y) is output.  Y and YDOT are arrays of length NEQ.
00492 !C             Subroutine F should not alter Y(1),...,Y(NEQ).  F must be
00493 !C             declared EXTERNAL in the calling program.
00494 !C
00495 !C             Subroutine F may access user-defined quantities in
00496 !C             NEQ(2),... and/or in Y(NEQ(1)+1),..., if NEQ is an array
00497 !C             (dimensioned in F) and/or Y has length exceeding NEQ(1).
00498 !C             See the descriptions of NEQ and Y below.
00499 !C
00500 !C             If quantities computed in the F routine are needed
00501 !C             externally to DLSODE, an extra call to F should be made
00502 !C             for this purpose, for consistent and accurate results.
00503 !C             If only the derivative dy/dt is needed, use DINTDY
00504 !C             instead.
00505 !C
00506 !C    NEQ      The size of the ODE system (number of first-order
00507 !C             ordinary differential equations).  Used only for input.
00508 !C             NEQ may be decreased, but not increased, during the
00509 !C             problem.  If NEQ is decreased (with ISTATE = 3 on input),
00510 !C             the remaining components of Y should be left undisturbed,
00511 !C             if these are to be accessed in F and/or JAC.
00512 !C
00513 !C             Normally, NEQ is a scalar, and it is generally referred
00514 !C             to as a scalar in this user interface description.
00515 !C             However, NEQ may be an array, with NEQ(1) set to the
00516 !C             system size.  (The DLSODE package accesses only NEQ(1).)
00517 !C             In either case, this parameter is passed as the NEQ
00518 !C             argument in all calls to F and JAC.  Hence, if it is an
00519 !C             array, locations NEQ(2),... may be used to store other
00520 !C             integer data and pass it to F and/or JAC.  Subroutines
00521 !C             F and/or JAC must include NEQ in a DIMENSION statement
00522 !C             in that case.
00523 !C
00524 !C    Y        A real array for the vector of dependent variables, of
00525 !C             length NEQ or more.  Used for both input and output on
00526 !C             the first call (ISTATE = 1), and only for output on
00527 !C             other calls.  On the first call, Y must contain the
00528 !C             vector of initial values.  On output, Y contains the
00529 !C             computed solution vector, evaluated at T. If desired,
00530 !C             the Y array may be used for other purposes between
00531 !C             calls to the solver.
00532 !C
00533 !C             This array is passed as the Y argument in all calls to F
00534 !C             and JAC.  Hence its length may exceed NEQ, and locations
00535 !C             Y(NEQ+1),... may be used to store other real data and
00536 !C             pass it to F and/or JAC.  (The DLSODE package accesses
00537 !C             only Y(1),...,Y(NEQ).)
00538 !C
00539 !C    T        The independent variable.  On input, T is used only on
00540 !C             the first call, as the initial point of the integration.
00541 !C             On output, after each call, T is the value at which a
00542 !C             computed solution Y is evaluated (usually the same as
00543 !C             TOUT).  On an error return, T is the farthest point
00544 !C             reached.
00545 !C
00546 !C    TOUT     The next value of T at which a computed solution is
00547 !C             desired.  Used only for input.
00548 !C
00549 !C             When starting the problem (ISTATE = 1), TOUT may be equal
00550 !C             to T for one call, then should not equal T for the next
00551 !C             call.  For the initial T, an input value of TOUT .NE. T
00552 !C             is used in order to determine the direction of the
00553 !C             integration (i.e., the algebraic sign of the step sizes)
00554 !C             and the rough scale of the problem.  Integration in
00555 !C             either direction (forward or backward in T) is permitted.
00556 !C
00557 !C             If ITASK = 2 or 5 (one-step modes), TOUT is ignored
00558 !C             after the first call (i.e., the first call with
00559 !C             TOUT .NE. T).  Otherwise, TOUT is required on every call.
00560 !C
00561 !C             If ITASK = 1, 3, or 4, the values of TOUT need not be
00562 !C             monotone, but a value of TOUT which backs up is limited
00563 !C             to the current internal T interval, whose endpoints are
00564 !C             TCUR - HU and TCUR.  (See "Optional Outputs" below for
00565 !C             TCUR and HU.)
00566 !C
00567 !C
00568 !C    ITOL     An indicator for the type of error control.  See
00569 !C             description below under ATOL.  Used only for input.
00570 !C
00571 !C    RTOL     A relative error tolerance parameter, either a scalar or
00572 !C             an array of length NEQ.  See description below under
00573 !C             ATOL.  Input only.
00574 !C
00575 !C    ATOL     An absolute error tolerance parameter, either a scalar or
00576 !C             an array of length NEQ.  Input only.
00577 !C
00578 !C             The input parameters ITOL, RTOL, and ATOL determine the
00579 !C             error control performed by the solver.  The solver will
00580 !C             control the vector e = (e(i)) of estimated local errors
00581 !C             in Y, according to an inequality of the form
00582 !C
00583 !C                rms-norm of ( e(i)/EWT(i) ) <= 1,
00584 !C
00585 !C             where
00586 !C
00587 !C                EWT(i) = RTOL(i)*ABS(Y(i)) + ATOL(i),
00588 !C
00589 !C             and the rms-norm (root-mean-square norm) here is
00590 !C
00591 !C                rms-norm(v) = SQRT(sum v(i)**2 / NEQ).
00592 !C
00593 !C             Here EWT = (EWT(i)) is a vector of weights which must
00594 !C             always be positive, and the values of RTOL and ATOL
00595 !C             should all be nonnegative.  The following table gives the
00596 !C             types (scalar/array) of RTOL and ATOL, and the
00597 !C             corresponding form of EWT(i).
00598 !C
00599 !C             ITOL    RTOL      ATOL      EWT(i)
00600 !C             ----    ------    ------    -----------------------------
00601 !C             1       scalar    scalar    RTOL*ABS(Y(i)) + ATOL
00602 !C             2       scalar    array     RTOL*ABS(Y(i)) + ATOL(i)
00603 !C             3       array     scalar    RTOL(i)*ABS(Y(i)) + ATOL
00604 !C             4       array     array     RTOL(i)*ABS(Y(i)) + ATOL(i)
00605 !C
00606 !C             When either of these parameters is a scalar, it need not
00607 !C             be dimensioned in the user's calling program.
00608 !C
00609 !C             If none of the above choices (with ITOL, RTOL, and ATOL
00610 !C             fixed throughout the problem) is suitable, more general
00611 !C             error controls can be obtained by substituting
00612 !C             user-supplied routines for the setting of EWT and/or for
00613 !C             the norm calculation.  See Part 4 below.
00614 !C
00615 !C             If global errors are to be estimated by making a repeated
00616 !C             run on the same problem with smaller tolerances, then all
00617 !C             components of RTOL and ATOL (i.e., of EWT) should be
00618 !C             scaled down uniformly.
00619 !C
00620 !C    ITASK    An index specifying the task to be performed.  Input
00621 !C             only.  ITASK has the following values and meanings:
00622 !C             1   Normal computation of output values of y(t) at
00623 !C                 t = TOUT (by overshooting and interpolating).
00624 !C             2   Take one step only and return.
00625 !C             3   Stop at the first internal mesh point at or beyond
00626 !C                 t = TOUT and return.
00627 !C             4   Normal computation of output values of y(t) at
00628 !C                 t = TOUT but without overshooting t = TCRIT.  TCRIT
00629 !C                 must be input as RWORK(1).  TCRIT may be equal to or
00630 !C                 beyond TOUT, but not behind it in the direction of
00631 !C                 integration.  This option is useful if the problem
00632 !C                 has a singularity at or beyond t = TCRIT.
00633 !C             5   Take one step, without passing TCRIT, and return.
00634 !C                 TCRIT must be input as RWORK(1).
00635 !C
00636 !C             Note:  If ITASK = 4 or 5 and the solver reaches TCRIT
00637 !C             (within roundoff), it will return T = TCRIT (exactly) to
00638 !C             indicate this (unless ITASK = 4 and TOUT comes before
00639 !C             TCRIT, in which case answers at T = TOUT are returned
00640 !C             first).
00641 !C
00642 !C    ISTATE   An index used for input and output to specify the state
00643 !C             of the calculation.
00644 !C
00645 !C             On input, the values of ISTATE are as follows:
00646 !C             1   This is the first call for the problem
00647 !C                 (initializations will be done).  See "Note" below.
00648 !C             2   This is not the first call, and the calculation is to
00649 !C                 continue normally, with no change in any input
00650 !C                 parameters except possibly TOUT and ITASK.  (If ITOL,
00651 !C                 RTOL, and/or ATOL are changed between calls with
00652 !C                 ISTATE = 2, the new values will be used but not
00653 !C                 tested for legality.)
00654 !C             3   This is not the first call, and the calculation is to
00655 !C                 continue normally, but with a change in input
00656 !C                 parameters other than TOUT and ITASK.  Changes are
00657 !C                 allowed in NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF,
00658 !C                 ML, MU, and any of the optional inputs except H0.
00659 !C                 (See IWORK description for ML and MU.)
00660 !C
00661 !C             Note:  A preliminary call with TOUT = T is not counted as
00662 !C             a first call here, as no initialization or checking of
00663 !C             input is done.  (Such a call is sometimes useful for the
00664 !C             purpose of outputting the initial conditions.)  Thus the
00665 !C             first call for which TOUT .NE. T requires ISTATE = 1 on
00666 !C             input.
00667 !C
00668 !C             On output, ISTATE has the following values and meanings:
00669 !C              1  Nothing was done, as TOUT was equal to T with
00670 !C                 ISTATE = 1 on input.
00671 !C              2  The integration was performed successfully.
00672 !C             -1  An excessive amount of work (more than MXSTEP steps)
00673 !C                 was done on this call, before completing the
00674 !C                 requested task, but the integration was otherwise
00675 !C                 successful as far as T. (MXSTEP is an optional input
00676 !C                 and is normally 500.)  To continue, the user may
00677 !C                 simply reset ISTATE to a value >1 and call again (the
00678 !C                 excess work step counter will be reset to 0).  In
00679 !C                 addition, the user may increase MXSTEP to avoid this
00680 !C                 error return; see "Optional Inputs" below.
00681 !C             -2  Too much accuracy was requested for the precision of
00682 !C                 the machine being used.  This was detected before
00683 !C                 completing the requested task, but the integration
00684 !C                 was successful as far as T. To continue, the
00685 !C                 tolerance parameters must be reset, and ISTATE must
00686 !C                 be set to 3. The optional output TOLSF may be used
00687 !C                 for this purpose.  (Note:  If this condition is
00688 !C                 detected before taking any steps, then an illegal
00689 !C                 input return (ISTATE = -3) occurs instead.)
00690 !C             -3  Illegal input was detected, before taking any
00691 !C                 integration steps.  See written message for details.
00692 !C                 (Note:  If the solver detects an infinite loop of
00693 !C                 calls to the solver with illegal input, it will cause
00694 !C                 the run to stop.)
00695 !C             -4  There were repeated error-test failures on one
00696 !C                 attempted step, before completing the requested task,
00697 !C                 but the integration was successful as far as T.  The
00698 !C                 problem may have a singularity, or the input may be
00699 !C                 inappropriate.
00700 !C             -5  There were repeated convergence-test failures on one
00701 !C                 attempted step, before completing the requested task,
00702 !C                 but the integration was successful as far as T. This
00703 !C                 may be caused by an inaccurate Jacobian matrix, if
00704 !C                 one is being used.
00705 !C             -6  EWT(i) became zero for some i during the integration.
00706 !C                 Pure relative error control (ATOL(i)=0.0) was
00707 !C                 requested on a variable which has now vanished.  The
00708 !C                 integration was successful as far as T.
00709 !C
00710 !C             Note:  Since the normal output value of ISTATE is 2, it
00711 !C             does not need to be reset for normal continuation.  Also,
00712 !C             since a negative input value of ISTATE will be regarded
00713 !C             as illegal, a negative output value requires the user to
00714 !C             change it, and possibly other inputs, before calling the
00715 !C             solver again.
00716 !C
00717 !C    IOPT     An integer flag to specify whether any optional inputs
00718 !C             are being used on this call.  Input only.  The optional
00719 !C             inputs are listed under a separate heading below.
00720 !C             0   No optional inputs are being used.  Default values
00721 !C                 will be used in all cases.
00722 !C             1   One or more optional inputs are being used.
00723 !C
00724 !C    RWORK    A real working array (double precision).  The length of
00725 !C             RWORK must be at least
00726 !C
00727 !C                20 + NYH*(MAXORD + 1) + 3*NEQ + LWM
00728 !C
00729 !C             where
00730 !C                NYH = the initial value of NEQ,
00731 !C             MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
00732 !C                      smaller value is given as an optional input),
00733 !C                LWM = 0           if MITER = 0,
00734 !C                LWM = NEQ**2 + 2  if MITER = 1 or 2,
00735 !C                LWM = NEQ + 2     if MITER = 3, and
00736 !C                LWM = (2*ML + MU + 1)*NEQ + 2
00737 !C                                  if MITER = 4 or 5.
00738 !C             (See the MF description below for METH and MITER.)
00739 !C
00740 !C             Thus if MAXORD has its default value and NEQ is constant,
00741 !C             this length is:
00742 !C             20 + 16*NEQ                    for MF = 10,
00743 !C             22 + 16*NEQ + NEQ**2           for MF = 11 or 12,
00744 !C             22 + 17*NEQ                    for MF = 13,
00745 !C             22 + 17*NEQ + (2*ML + MU)*NEQ  for MF = 14 or 15,
00746 !C             20 +  9*NEQ                    for MF = 20,
00747 !C             22 +  9*NEQ + NEQ**2           for MF = 21 or 22,
00748 !C             22 + 10*NEQ                    for MF = 23,
00749 !C             22 + 10*NEQ + (2*ML + MU)*NEQ  for MF = 24 or 25.
00750 !C
00751 !C             The first 20 words of RWORK are reserved for conditional
00752 !C             and optional inputs and optional outputs.
00753 !C
00754 !C             The following word in RWORK is a conditional input:
00755 !C             RWORK(1) = TCRIT, the critical value of t which the
00756 !C                        solver is not to overshoot.  Required if ITASK
00757 !C                        is 4 or 5, and ignored otherwise.  See ITASK.
00758 !C
00759 !C    LRW      The length of the array RWORK, as declared by the user.
00760 !C             (This will be checked by the solver.)
00761 !C
00762 !C    IWORK    An integer work array.  Its length must be at least
00763 !C             20       if MITER = 0 or 3 (MF = 10, 13, 20, 23), or
00764 !C             20 + NEQ otherwise (MF = 11, 12, 14, 15, 21, 22, 24, 25).
00765 !C             (See the MF description below for MITER.)  The first few
00766 !C             words of IWORK are used for conditional and optional
00767 !C             inputs and optional outputs.
00768 !C
00769 !C             The following two words in IWORK are conditional inputs:
00770 !C             IWORK(1) = ML   These are the lower and upper half-
00771 !C             IWORK(2) = MU   bandwidths, respectively, of the banded
00772 !C                             Jacobian, excluding the main diagonal.
00773 !C                        The band is defined by the matrix locations
00774 !C                        (i,j) with i - ML <= j <= i + MU. ML and MU
00775 !C                        must satisfy 0 <= ML,MU <= NEQ - 1. These are
00776 !C                        required if MITER is 4 or 5, and ignored
00777 !C                        otherwise.  ML and MU may in fact be the band
00778 !C                        parameters for a matrix to which df/dy is only
00779 !C                        approximately equal.
00780 !C
00781 !C    LIW      The length of the array IWORK, as declared by the user.
00782 !C             (This will be checked by the solver.)
00783 !C
00784 !C    Note:  The work arrays must not be altered between calls to DLSODE
00785 !C    for the same problem, except possibly for the conditional and
00786 !C    optional inputs, and except for the last 3*NEQ words of RWORK.
00787 !C    The latter space is used for internal scratch space, and so is
00788 !C    available for use by the user outside DLSODE between calls, if
00789 !C    desired (but not for use by F or JAC).
00790 !C
00791 !C    JAC      The name of the user-supplied routine (MITER = 1 or 4) to
00792 !C             compute the Jacobian matrix, df/dy, as a function of the
00793 !C             scalar t and the vector y.  (See the MF description below
00794 !C             for MITER.)  It is to have the form
00795 !C
00796 !C                SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
00797 !C                DOUBLE PRECISION T, Y(*), PD(NROWPD,*)
00798 !C
00799 !C             where NEQ, T, Y, ML, MU, and NROWPD are input and the
00800 !C             array PD is to be loaded with partial derivatives
00801 !C             (elements of the Jacobian matrix) on output.  PD must be
00802 !C             given a first dimension of NROWPD.  T and Y have the same
00803 !C             meaning as in subroutine F.
00804 !C
00805 !C             In the full matrix case (MITER = 1), ML and MU are
00806 !C             ignored, and the Jacobian is to be loaded into PD in
00807 !C             columnwise manner, with df(i)/dy(j) loaded into PD(i,j).
00808 !C
00809 !C             In the band matrix case (MITER = 4), the elements within
00810 !C             the band are to be loaded into PD in columnwise manner,
00811 !C             with diagonal lines of df/dy loaded into the rows of PD.
00812 !C             Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j).  ML
00813 !C             and MU are the half-bandwidth parameters (see IWORK).
00814 !C             The locations in PD in the two triangular areas which
00815 !C             correspond to nonexistent matrix elements can be ignored
00816 !C             or loaded arbitrarily, as they are overwritten by DLSODE.
00817 !C
00818 !C             JAC need not provide df/dy exactly. A crude approximation
00819 !C             (possibly with a smaller bandwidth) will do.
00820 !C
00821 !C             In either case, PD is preset to zero by the solver, so
00822 !C             that only the nonzero elements need be loaded by JAC.
00823 !C             Each call to JAC is preceded by a call to F with the same
00824 !C             arguments NEQ, T, and Y. Thus to gain some efficiency,
00825 !C             intermediate quantities shared by both calculations may
00826 !C             be saved in a user COMMON block by F and not recomputed
00827 !C             by JAC, if desired.  Also, JAC may alter the Y array, if
00828 !C             desired.  JAC must be declared EXTERNAL in the calling
00829 !C             program.
00830 !C
00831 !C             Subroutine JAC may access user-defined quantities in
00832 !C             NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array
00833 !C             (dimensioned in JAC) and/or Y has length exceeding
00834 !C             NEQ(1).  See the descriptions of NEQ and Y above.
00835 !C
00836 !C    MF       The method flag.  Used only for input.  The legal values
00837 !C             of MF are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24,
00838 !C             and 25.  MF has decimal digits METH and MITER:
00839 !C                MF = 10*METH + MITER .
00840 !C
00841 !C             METH indicates the basic linear multistep method:
00842 !C             1   Implicit Adams method.
00843 !C             2   Method based on backward differentiation formulas
00844 !C                 (BDF's).
00845 !C
00846 !C             MITER indicates the corrector iteration method:
00847 !C             0   Functional iteration (no Jacobian matrix is
00848 !C                 involved).
00849 !C             1   Chord iteration with a user-supplied full (NEQ by
00850 !C                 NEQ) Jacobian.
00851 !C             2   Chord iteration with an internally generated
00852 !C                 (difference quotient) full Jacobian (using NEQ
00853 !C                 extra calls to F per df/dy value).
00854 !C             3   Chord iteration with an internally generated
00855 !C                 diagonal Jacobian approximation (using one extra call
00856 !C                 to F per df/dy evaluation).
00857 !C             4   Chord iteration with a user-supplied banded Jacobian.
00858 !C             5   Chord iteration with an internally generated banded
00859 !C                 Jacobian (using ML + MU + 1 extra calls to F per
00860 !C                 df/dy evaluation).
00861 !C
00862 !C             If MITER = 1 or 4, the user must supply a subroutine JAC
00863 !C             (the name is arbitrary) as described above under JAC.
00864 !C             For other values of MITER, a dummy argument can be used.
00865 !C
00866 !C    Optional Inputs
00867 !C    ---------------
00868 !C    The following is a list of the optional inputs provided for in the
00869 !C    call sequence.  (See also Part 2.)  For each such input variable,
00870 !C    this table lists its name as used in this documentation, its
00871 !C    location in the call sequence, its meaning, and the default value.
00872 !C    The use of any of these inputs requires IOPT = 1, and in that case
00873 !C    all of these inputs are examined.  A value of zero for any of
00874 !C    these optional inputs will cause the default value to be used.
00875 !C    Thus to use a subset of the optional inputs, simply preload
00876 !C    locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively,
00877 !C    and then set those of interest to nonzero values.
00878 !C
00879 !C    Name    Location   Meaning and default value
00880 !C    ------  ---------  -----------------------------------------------
00881 !C    H0      RWORK(5)   Step size to be attempted on the first step.
00882 !C                       The default value is determined by the solver.
00883 !C    HMAX    RWORK(6)   Maximum absolute step size allowed.  The
00884 !C                       default value is infinite.
00885 !C    HMIN    RWORK(7)   Minimum absolute step size allowed.  The
00886 !C                       default value is 0.  (This lower bound is not
00887 !C                       enforced on the final step before reaching
00888 !C                       TCRIT when ITASK = 4 or 5.)
00889 !C    MAXORD  IWORK(5)   Maximum order to be allowed.  The default value
00890 !C                       is 12 if METH = 1, and 5 if METH = 2. (See the
00891 !C                       MF description above for METH.)  If MAXORD
00892 !C                       exceeds the default value, it will be reduced
00893 !C                       to the default value.  If MAXORD is changed
00894 !C                       during the problem, it may cause the current
00895 !C                       order to be reduced.
00896 !C    MXSTEP  IWORK(6)   Maximum number of (internally defined) steps
00897 !C                       allowed during one call to the solver.  The
00898 !C                       default value is 500.
00899 !C    MXHNIL  IWORK(7)   Maximum number of messages printed (per
00900 !C                       problem) warning that T + H = T on a step
00901 !C                       (H = step size).  This must be positive to
00902 !C                       result in a nondefault value.  The default
00903 !C                       value is 10.
00904 !C
00905 !C    Optional Outputs
00906 !C    ----------------
00907 !C    As optional additional output from DLSODE, the variables listed
00908 !C    below are quantities related to the performance of DLSODE which
00909 !C    are available to the user.  These are communicated by way of the
00910 !C    work arrays, but also have internal mnemonic names as shown.
00911 !C    Except where stated otherwise, all of these outputs are defined on
00912 !C    any successful return from DLSODE, and on any return with ISTATE =
00913 !C    -1, -2, -4, -5, or -6.  On an illegal input return (ISTATE = -3),
00914 !C    they will be unchanged from their existing values (if any), except
00915 !C    possibly for TOLSF, LENRW, and LENIW.  On any error return,
00916 !C    outputs relevant to the error will be defined, as noted below.
00917 !C
00918 !C    Name   Location   Meaning
00919 !C    -----  ---------  ------------------------------------------------
00920 !C    HU     RWORK(11)  Step size in t last used (successfully).
00921 !C    HCUR   RWORK(12)  Step size to be attempted on the next step.
00922 !C    TCUR   RWORK(13)  Current value of the independent variable which
00923 !C                      the solver has actually reached, i.e., the
00924 !C                      current internal mesh point in t. On output,
00925 !C                      TCUR will always be at least as far as the
00926 !C                      argument T, but may be farther (if interpolation
00927 !C                      was done).
00928 !C    TOLSF  RWORK(14)  Tolerance scale factor, greater than 1.0,
00929 !C                      computed when a request for too much accuracy
00930 !C                      was detected (ISTATE = -3 if detected at the
00931 !C                      start of the problem, ISTATE = -2 otherwise).
00932 !C                      If ITOL is left unaltered but RTOL and ATOL are
00933 !C                      uniformly scaled up by a factor of TOLSF for the
00934 !C                      next call, then the solver is deemed likely to
00935 !C                      succeed.  (The user may also ignore TOLSF and
00936 !C                      alter the tolerance parameters in any other way
00937 !C                      appropriate.)
00938 !C    NST    IWORK(11)  Number of steps taken for the problem so far.
00939 !C    NFE    IWORK(12)  Number of F evaluations for the problem so far.
00940 !C    NJE    IWORK(13)  Number of Jacobian evaluations (and of matrix LU
00941 !C                      decompositions) for the problem so far.
00942 !C    NQU    IWORK(14)  Method order last used (successfully).
00943 !C    NQCUR  IWORK(15)  Order to be attempted on the next step.
00944 !C    IMXER  IWORK(16)  Index of the component of largest magnitude in
00945 !C                      the weighted local error vector ( e(i)/EWT(i) ),
00946 !C                      on an error return with ISTATE = -4 or -5.
00947 !C    LENRW  IWORK(17)  Length of RWORK actually required.  This is
00948 !C                      defined on normal returns and on an illegal
00949 !C                      input return for insufficient storage.
00950 !C    LENIW  IWORK(18)  Length of IWORK actually required.  This is
00951 !C                      defined on normal returns and on an illegal
00952 !C                      input return for insufficient storage.
00953 !C
00954 !C    The following two arrays are segments of the RWORK array which may
00955 !C    also be of interest to the user as optional outputs.  For each
00956 !C    array, the table below gives its internal name, its base address
00957 !C    in RWORK, and its description.
00958 !C
00959 !C    Name  Base address  Description
00960 !C    ----  ------------  ----------------------------------------------
00961 !C    YH    21            The Nordsieck history array, of size NYH by
00962 !C                        (NQCUR + 1), where NYH is the initial value of
00963 !C                        NEQ.  For j = 0,1,...,NQCUR, column j + 1 of
00964 !C                        YH contains HCUR**j/factorial(j) times the jth
00965 !C                        derivative of the interpolating polynomial
00966 !C                        currently representing the solution, evaluated
00967 !C                        at t = TCUR.
00968 !C    ACOR  LENRW-NEQ+1   Array of size NEQ used for the accumulated
00969 !C                        corrections on each step, scaled on output to
00970 !C                        represent the estimated local error in Y on
00971 !C                        the last step.  This is the vector e in the
00972 !C                        description of the error control.  It is
00973 !C                        defined only on successful return from DLSODE.
00974 !C
00975 !C
00976 !C                   Part 2.  Other Callable Routines
00977 !C                   --------------------------------
00978 !C
00979 !C    The following are optional calls which the user may make to gain
00980 !C    additional capabilities in conjunction with DLSODE.
00981 !C
00982 !C    Form of call              Function
00983 !C    ------------------------  ----------------------------------------
00984 !C    CALL XSETUN(LUN)          Set the logical unit number, LUN, for
00985 !C                              output of messages from DLSODE, if the
00986 !C                              default is not desired.  The default
00987 !C                              value of LUN is 6. This call may be made
00988 !C                              at any time and will take effect
00989 !C                              immediately.
00990 !C    CALL XSETF(MFLAG)         Set a flag to control the printing of
00991 !C                              messages by DLSODE.  MFLAG = 0 means do
00992 !C                              not print.  (Danger:  this risks losing
00993 !C                              valuable information.)  MFLAG = 1 means
00994 !C                              print (the default).  This call may be
00995 !C                              made at any time and will take effect
00996 !C                              immediately.
00997 !C    CALL DSRCOM(RSAV,ISAV,JOB)  Saves and restores the contents of the
00998 !C                              internal COMMON blocks used by DLSODE
00999 !C                              (see Part 3 below).  RSAV must be a
01000 !C                              real array of length 218 or more, and
01001 !C                              ISAV must be an integer array of length
01002 !C                              37 or more.  JOB = 1 means save COMMON
01003 !C                              into RSAV/ISAV.  JOB = 2 means restore
01004 !C                              COMMON from same.  DSRCOM is useful if
01005 !C                              one is interrupting a run and restarting
01006 !C                              later, or alternating between two or
01007 !C                              more problems solved with DLSODE.
01008 !C    CALL DINTDY(,,,,,)        Provide derivatives of y, of various
01009 !C    (see below)               orders, at a specified point t, if
01010 !C                              desired.  It may be called only after a
01011 !C                              successful return from DLSODE.  Detailed
01012 !C                              instructions follow.
01013 !C
01014 !C    Detailed instructions for using DINTDY
01015 !C    --------------------------------------
01016 !C    The form of the CALL is:
01017 !C
01018 !C          CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
01019 !C
01020 !C    The input parameters are:
01021 !C
01022 !C    T          Value of independent variable where answers are
01023 !C               desired (normally the same as the T last returned by
01024 !C               DLSODE).  For valid results, T must lie between
01025 !C               TCUR - HU and TCUR.  (See "Optional Outputs" above
01026 !C               for TCUR and HU.)
01027 !C    K          Integer order of the derivative desired.  K must
01028 !C               satisfy 0 <= K <= NQCUR, where NQCUR is the current
01029 !C               order (see "Optional Outputs").  The capability
01030 !C               corresponding to K = 0, i.e., computing y(t), is
01031 !C               already provided by DLSODE directly.  Since
01032 !C               NQCUR >= 1, the first derivative dy/dt is always
01033 !C               available with DINTDY.
01034 !C    RWORK(21)  The base address of the history array YH.
01035 !C    NYH        Column length of YH, equal to the initial value of NEQ.
01036 !C
01037 !C    The output parameters are:
01038 !C
01039 !C    DKY        Real array of length NEQ containing the computed value
01040 !C               of the Kth derivative of y(t).
01041 !C    IFLAG      Integer flag, returned as 0 if K and T were legal,
01042 !C               -1 if K was illegal, and -2 if T was illegal.
01043 !C               On an error return, a message is also written.
01044 !C
01045 !C
01046 !C                         Part 3.  Common Blocks
01047 !C                         ----------------------
01048 !C
01049 !C    If DLSODE is to be used in an overlay situation, the user must
01050 !C    declare, in the primary overlay, the variables in:
01051 !C    (1) the call sequence to DLSODE,
01052 !C    (2) the internal COMMON block /DLS001/, of length 255
01053 !C        (218 double precision words followed by 37 integer words).
01054 !C
01055 !C    If DLSODE is used on a system in which the contents of internal
01056 !C    COMMON blocks are not preserved between calls, the user should
01057 !C    declare the above COMMON block in his main program to insure that
01058 !C    its contents are preserved.
01059 !C
01060 !C    If the solution of a given problem by DLSODE is to be interrupted
01061 !C    and then later continued, as when restarting an interrupted run or
01062 !C    alternating between two or more problems, the user should save,
01063 !C    following the return from the last DLSODE call prior to the
01064 !C    interruption, the contents of the call sequence variables and the
01065 !C    internal COMMON block, and later restore these values before the
01066 !C    next DLSODE call for that problem.   In addition, if XSETUN and/or
01067 !C    XSETF was called for non-default handling of error messages, then
01068 !C    these calls must be repeated.  To save and restore the COMMON
01069 !C    block, use subroutine DSRCOM (see Part 2 above).
01070 !C
01071 !C
01072 !C             Part 4.  Optionally Replaceable Solver Routines
01073 !C             -----------------------------------------------
01074 !C
01075 !C    Below are descriptions of two routines in the DLSODE package which
01076 !C    relate to the measurement of errors.  Either routine can be
01077 !C    replaced by a user-supplied version, if desired.  However, since
01078 !C    such a replacement may have a major impact on performance, it
01079 !C    should be done only when absolutely necessary, and only with great
01080 !C    caution.  (Note:  The means by which the package version of a
01081 !C    routine is superseded by the user's version may be system-
01082 !C    dependent.)
01083 !C
01084 !C    DEWSET
01085 !C    ------
01086 !C    The following subroutine is called just before each internal
01087 !C    integration step, and sets the array of error weights, EWT, as
01088 !C    described under ITOL/RTOL/ATOL above:
01089 !C
01090 !C          SUBROUTINE DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT)
01091 !C
01092 !C    where NEQ, ITOL, RTOL, and ATOL are as in the DLSODE call
01093 !C    sequence, YCUR contains the current dependent variable vector,
01094 !C    and EWT is the array of weights set by DEWSET.
01095 !C
01096 !C    If the user supplies this subroutine, it must return in EWT(i)
01097 !C    (i = 1,...,NEQ) a positive quantity suitable for comparing errors
01098 !C    in Y(i) to.  The EWT array returned by DEWSET is passed to the
01099 !C    DVNORM routine (see below), and also used by DLSODE in the
01100 !C    computation of the optional output IMXER, the diagonal Jacobian
01101 !C    approximation, and the increments for difference quotient
01102 !C    Jacobians.
01103 !C
01104 !C    In the user-supplied version of DEWSET, it may be desirable to use
01105 !C    the current values of derivatives of y. Derivatives up to order NQ
01106 !C    are available from the history array YH, described above under
01107 !C    optional outputs.  In DEWSET, YH is identical to the YCUR array,
01108 !C    extended to NQ + 1 columns with a column length of NYH and scale
01109 !C    factors of H**j/factorial(j).  On the first call for the problem,
01110 !C    given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
01111 !C    NYH is the initial value of NEQ.  The quantities NQ, H, and NST
01112 !C    can be obtained by including in SEWSET the statements:
01113 !C          DOUBLE PRECISION RLS
01114 !C          COMMON /DLS001/ RLS(218),ILS(37)
01115 !C          NQ = ILS(33)
01116 !C          NST = ILS(34)
01117 !C          H = RLS(212)
01118 !C    Thus, for example, the current value of dy/dt can be obtained as
01119 !C    YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is unnecessary
01120 !C    when NST = 0).
01121 !C
01122 !C    DVNORM
01123 !C    ------
01124 !C    DVNORM is a real function routine which computes the weighted
01125 !C    root-mean-square norm of a vector v:
01126 !C
01127 !C       d = DVNORM (n, v, w)
01128 !C
01129 !C    where:
01130 !C    n = the length of the vector,
01131 !C    v = real array of length n containing the vector,
01132 !C    w = real array of length n containing weights,
01133 !C    d = SQRT( (1/n) * sum(v(i)*w(i))**2 ).
01134 !C
01135 !C    DVNORM is called with n = NEQ and with w(i) = 1.0/EWT(i), where
01136 !C    EWT is as set by subroutine DEWSET.
01137 !C
01138 !C    If the user supplies this function, it should return a nonnegative
01139 !C    value of DVNORM suitable for use in the error control in DLSODE.
01140 !C    None of the arguments should be altered by DVNORM.  For example, a
01141 !C    user-supplied DVNORM routine might:
01142 !C    - Substitute a max-norm of (v(i)*w(i)) for the rms-norm, or
01143 !C    - Ignore some components of v in the norm, with the effect of
01144 !C      suppressing the error control on those components of Y.
01145 !C ---------------------------------------------------------------------
01146 !C**ROUTINES CALLED  DEWSET, DINTDY, DUMACH, DSTODE, DVNORM, XERRWD
01147 !C**COMMON BLOCKS    DLS001
01148 !C**REVISION HISTORY  (YYYYMMDD)
01149 !C19791129  DATE WRITTEN
01150 !C19791213  Minor changes to declarations; DELP init. in STODE.
01151 !C19800118  Treat NEQ as array; integer declarations added throughout;
01152 !C          minor changes to prologue.
01153 !C19800306  Corrected TESCO(1,NQP1) setting in CFODE.
01154 !C19800519  Corrected access of YH on forced order reduction;
01155 !C          numerous corrections to prologues and other comments.
01156 !C19800617  In main driver, added loading of SQRT(UROUND) in RWORK;
01157 !C          minor corrections to main prologue.
01158 !C19800923  Added zero initialization of HU and NQU.
01159 !C19801218  Revised XERRWD routine; minor corrections to main prologue.
01160 !C19810401  Minor changes to comments and an error message.
01161 !C19810814  Numerous revisions: replaced EWT by 1/EWT; used flags
01162 !C          JCUR, ICF, IERPJ, IERSL between STODE and subordinates;
01163 !C          added tuning parameters CCMAX, MAXCOR, MSBP, MXNCF;
01164 !C          reorganized returns from STODE; reorganized type decls.;
01165 !C          fixed message length in XERRWD; changed default LUNIT to 6;
01166 !C          changed Common lengths; changed comments throughout.
01167 !C19870330  Major update by ACH: corrected comments throughout;
01168 !C          removed TRET from Common; rewrote EWSET with 4 loops;
01169 !C          fixed t test in INTDY; added Cray directives in STODE;
01170 !C          in STODE, fixed DELP init. and logic around PJAC call;
01171 !C          combined routines to save/restore Common;
01172 !C          passed LEVEL = 0 in error message calls (except run abort).
01173 !C19890426  Modified prologue to SLATEC/LDOC format.  (FNF)
01174 !C19890501  Many improvements to prologue.  (FNF)
01175 !C19890503  A few final corrections to prologue.  (FNF)
01176 !C19890504  Minor cosmetic changes.  (FNF)
01177 !C19890510  Corrected description of Y in Arguments section.  (FNF)
01178 !C19890517  Minor corrections to prologue.  (FNF)
01179 !C19920514  Updated with prologue edited 891025 by G. Shaw for manual.
01180 !C19920515  Converted source lines to upper case.  (FNF)
01181 !C19920603  Revised XERRWD calls using mixed upper-lower case.  (ACH)
01182 !C19920616  Revised prologue comment regarding CFT.  (ACH)
01183 !C19921116  Revised prologue comments regarding Common.  (ACH).
01184 !C19930326  Added comment about non-reentrancy.  (FNF)
01185 !C19930723  Changed D1MACH to DUMACH. (FNF)
01186 !C19930801  Removed ILLIN and NTREP from Common (affects driver logic);
01187 !C          minor changes to prologue and internal comments;
01188 !C          changed Hollerith strings to quoted strings;
01189 !C          changed internal comments to mixed case;
01190 !C          replaced XERRWD with new version using character type;
01191 !C          changed dummy dimensions from 1 to *. (ACH)
01192 !C19930809  Changed to generic intrinsic names; changed names of
01193 !C          subprograms and Common blocks to DLSODE etc. (ACH)
01194 !C19930929  Eliminated use of REAL intrinsic; other minor changes. (ACH)
01195 !C20010412  Removed all 'own' variables from Common block /DLS001/
01196 !C          (affects declarations in 6 routines). (ACH)
01197 !C20010509  Minor corrections to prologue. (ACH)
01198 !C20031105  Restored 'own' variables to Common block /DLS001/, to
01199 !C          enable interrupt/restart feature. (ACH)
01200 !C20031112  Added SAVE statements for data-loaded constants.
01201 !C
01202 !C**END PROLOGUE  DLSODE
01203 !C
01204 !CInternal Notes:
01205 !C
01206 !COther Routines in the DLSODE Package.
01207 !C
01208 !CIn addition to Subroutine DLSODE, the DLSODE package includes the
01209 !Cfollowing subroutines and function routines:
01210 !C DINTDY   computes an interpolated value of the y vector at t = TOUT.
01211 !C DSTODE   is the core integrator, which does one step of the
01212 !C          integration and the associated error control.
01213 !C DCFODE   sets all method coefficients and test constants.
01214 !C DPREPJ   computes and preprocesses the Jacobian matrix J = df/dy
01215 !C          and the Newton iteration matrix P = I - h*l0*J.
01216 !C DSOLSY   manages solution of linear system in chord iteration.
01217 !C DEWSET   sets the error weight vector EWT before each step.
01218 !C DVNORM   computes the weighted R.M.S. norm of a vector.
01219 !C DSRCOM   is a user-callable routine to save and restore
01220 !C          the contents of the internal Common block.
01221 !C DGEFA and DGESL   are routines from LINPACK for solving full
01222 !C          systems of linear algebraic equations.
01223 !C DGBFA and DGBSL   are routines from LINPACK for solving banded
01224 !C          linear systems.
01225 !C DUMACH   computes the unit roundoff in a machine-independent manner.
01226 !C XERRWD, XSETUN, XSETF, IXSAV, IUMACH   handle the printing of all
01227 !C          error messages and warnings.  XERRWD is machine-dependent.
01228 !CNote: DVNORM, DUMACH, IXSAV, and IUMACH are function routines.
01229 !CAll the others are subroutines.
01230 !C
01231 !C*End
01232 !C
01233 !C Declare externals.
01234       EXTERNAL DPREPJ, DSOLSY
01235       DOUBLE PRECISION DUMACH, DVNORM
01236 !C
01237 !C Declare all other variables.
01238       INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH, IOWNS, 
01239          ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 
01240          LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 
01241          MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
01242       INTEGER I, I1, I2, IFLAG, IMXER, KGO, LF0, 
01243          LENIW, LENRW, LENWM, ML, MORD, MU, MXHNL0, MXSTP0
01244       DOUBLE PRECISION ROWNS, 
01245          CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
01246       DOUBLE PRECISION ATOLI, AYI, BIG, EWTI, H0, HMAX, HMX, RH, RTOLI, 
01247          TCRIT, TDIST, TNEXT, TOL, TOLSF, TP, SIZE, SUM, W0
01248       DIMENSION MORD(2)
01249       LOGICAL IHIT
01250       CHARACTER*80 MSG
01251       SAVE MORD, MXSTP0, MXHNL0
01252 !C----------------------------------------------------------------------
01253 !CThe following internal Common block contains
01254 !C(a) variables which are local to any subroutine but whose values must
01255 !C    be preserved between calls to the routine ("own" variables), and
01256 !C(b) variables which are communicated between subroutines.
01257 !CThe block DLS001 is declared in subroutines DLSODE, DINTDY, DSTODE,
01258 !CDPREPJ, and DSOLSY.
01259 !CGroups of variables are replaced by dummy arrays in the Common
01260 !Cdeclarations in routines where those variables are not used.
01261 !C----------------------------------------------------------------------
01262       COMMON /DLS001/ ROWNS(209), &
01263          CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
01264          INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH, IOWNS(6), &
01265          ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
01266          LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
01267          MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
01268 !C
01269       DATA  MORD(1),MORD(2)/12,5/, MXSTP0/10000/, MXHNL0/10/
01270 !C----------------------------------------------------------------------
01271 !CBlock A.
01272 !CThis code block is executed on every call.
01273 !CIt tests ISTATE and ITASK for legality and branches appropriately.
01274 !CIf ISTATE .GT. 1 but the flag INIT shows that initialization has
01275 !Cnot yet been done, an error return occurs.
01276 !CIf ISTATE = 1 and TOUT = T, return immediately.
01277 !C----------------------------------------------------------------------
01278 !C
01279 !C**FIRST EXECUTABLE STATEMENT  DLSODE
01280       IF (ISTATE .LT. 1 .OR. ISTATE .GT. 3) GO TO 601
01281       IF (ITASK .LT. 1 .OR. ITASK .GT. 5) GO TO 602
01282       IF (ISTATE .EQ. 1) GO TO 10
01283       IF (INIT .EQ. 0) GO TO 603
01284       IF (ISTATE .EQ. 2) GO TO 200
01285       GO TO 20
01286  10   INIT = 0
01287       IF (TOUT .EQ. T) RETURN
01288 !C----------------------------------------------------------------------
01289 !CBlock B.
01290 !CThe next code block is executed for the initial call (ISTATE = 1),
01291 !Cor for a continuation call with parameter changes (ISTATE = 3).
01292 !CIt contains checking of all inputs and various initializations.
01293 !C
01294 !CFirst check legality of the non-optional inputs NEQ, ITOL, IOPT,
01295 !CMF, ML, and MU.
01296 !C----------------------------------------------------------------------
01297  20   IF (NEQ(1) .LE. 0) GO TO 604
01298       IF (ISTATE .EQ. 1) GO TO 25
01299       IF (NEQ(1) .GT. N) GO TO 605
01300  25   N = NEQ(1)
01301       IF (ITOL .LT. 1 .OR. ITOL .GT. 4) GO TO 606
01302       IF (IOPT .LT. 0 .OR. IOPT .GT. 1) GO TO 607
01303       METH = MF/10
01304       MITER = MF - 10*METH
01305       IF (METH .LT. 1 .OR. METH .GT. 2) GO TO 608
01306       IF (MITER .LT. 0 .OR. MITER .GT. 5) GO TO 608
01307       IF (MITER .LE. 3) GO TO 30
01308       ML = IWORK(1)
01309       MU = IWORK(2)
01310       IF (ML .LT. 0 .OR. ML .GE. N) GO TO 609
01311       IF (MU .LT. 0 .OR. MU .GE. N) GO TO 610
01312  30   CONTINUE
01313 !CNext process and check the optional inputs. --------------------------
01314       IF (IOPT .EQ. 1) GO TO 40
01315       MAXORD = MORD(METH)
01316       MXSTEP = MXSTP0
01317       MXHNIL = MXHNL0
01318       IF (ISTATE .EQ. 1) H0 = 0.0D0
01319       HMXI = 0.0D0
01320       HMIN = 0.0D0
01321       GO TO 60
01322  40   MAXORD = IWORK(5)
01323       IF (MAXORD .LT. 0) GO TO 611
01324       IF (MAXORD .EQ. 0) MAXORD = 100
01325       MAXORD = MIN(MAXORD,MORD(METH))
01326       MXSTEP = IWORK(6)
01327       IF (MXSTEP .LT. 0) GO TO 612
01328       IF (MXSTEP .EQ. 0) MXSTEP = MXSTP0
01329       MXHNIL = IWORK(7)
01330       IF (MXHNIL .LT. 0) GO TO 613
01331       IF (MXHNIL .EQ. 0) MXHNIL = MXHNL0
01332       IF (ISTATE .NE. 1) GO TO 50
01333       H0 = RWORK(5)
01334       IF ((TOUT - T)*H0 .LT. 0.0D0) GO TO 614
01335  50   HMAX = RWORK(6)
01336       IF (HMAX .LT. 0.0D0) GO TO 615
01337       HMXI = 0.0D0
01338       IF (HMAX .GT. 0.0D0) HMXI = 1.0D0/HMAX
01339       HMIN = RWORK(7)
01340       IF (HMIN .LT. 0.0D0) GO TO 616
01341 !C----------------------------------------------------------------------
01342 !CSet work array pointers and check lengths LRW and LIW.
01343 !CPointers to segments of RWORK and IWORK are named by prefixing L to
01344 !Cthe name of the segment.  E.g., the segment YH starts at RWORK(LYH).
01345 !CSegments of RWORK (in order) are denoted  YH, WM, EWT, SAVF, ACOR.
01346 !C----------------------------------------------------------------------
01347  60   LYH = 21
01348       IF (ISTATE .EQ. 1) NYH = N
01349       LWM = LYH + (MAXORD + 1)*NYH
01350       IF (MITER .EQ. 0) LENWM = 0
01351       IF (MITER .EQ. 1 .OR. MITER .EQ. 2) LENWM = N*N + 2
01352       IF (MITER .EQ. 3) LENWM = N + 2
01353       IF (MITER .GE. 4) LENWM = (2*ML + MU + 1)*N + 2
01354       LEWT = LWM + LENWM
01355       LSAVF = LEWT + N
01356       LACOR = LSAVF + N
01357       LENRW = LACOR + N - 1
01358       IWORK(17) = LENRW
01359       LIWM = 1
01360       LENIW = 20 + N
01361       IF (MITER .EQ. 0 .OR. MITER .EQ. 3) LENIW = 20
01362       IWORK(18) = LENIW
01363       IF (LENRW .GT. LRW) GO TO 617
01364       IF (LENIW .GT. LIW) GO TO 618
01365 !CCheck RTOL and ATOL for legality. ------------------------------------
01366       RTOLI = RTOL(1)
01367       ATOLI = ATOL(1)
01368       DO 70 I = 1,N
01369         IF (ITOL .GE. 3) RTOLI = RTOL(I)
01370         IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
01371         IF (RTOLI .LT. 0.0D0) GO TO 619
01372         IF (ATOLI .LT. 0.0D0) GO TO 620
01373  70     CONTINUE
01374       IF (ISTATE .EQ. 1) GO TO 100
01375 !CIf ISTATE = 3, set flag to signal parameter changes to DSTODE. -------
01376       JSTART = -1
01377       IF (NQ .LE. MAXORD) GO TO 90
01378 !CMAXORD was reduced below NQ.  Copy YH(*,MAXORD+2) into SAVF. ---------
01379       DO 80 I = 1,N
01380  80     RWORK(I+LSAVF-1) = RWORK(I+LWM-1)
01381 !CReload WM(1) = RWORK(LWM), since LWM may have changed. ---------------
01382  90   IF (MITER .GT. 0) RWORK(LWM) = DSQRT(UROUND)
01383       IF (N .EQ. NYH) GO TO 200
01384 !CNEQ was reduced.  Zero part of YH to avoid undefined references. -----
01385       I1 = LYH + L*NYH
01386       I2 = LYH + (MAXORD + 1)*NYH - 1
01387       IF (I1 .GT. I2) GO TO 200
01388       DO 95 I = I1,I2
01389  95     RWORK(I) = 0.0D0
01390       GO TO 200
01391 !C----------------------------------------------------------------------
01392 !CBlock C.
01393 !CThe next block is for the initial call only (ISTATE = 1).
01394 !CIt contains all remaining initializations, the initial call to F,
01395 !Cand the calculation of the initial step size.
01396 !CThe error weights in EWT are inverted after being loaded.
01397 !C----------------------------------------------------------------------
01398  100  UROUND = DUMACH()
01399       TN = T
01400       IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 110
01401       TCRIT = RWORK(1)
01402       IF ((TCRIT - TOUT)*(TOUT - T) .LT. 0.0D0) GO TO 625
01403       IF (H0 .NE. 0.0D0 .AND. (T + H0 - TCRIT)*H0 .GT. 0.0D0) &
01404          H0 = TCRIT - T
01405  110  JSTART = 0
01406       IF (MITER .GT. 0) RWORK(LWM) = DSQRT(UROUND)
01407       NHNIL = 0
01408       NST = 0
01409       NJE = 0
01410       NSLAST = 0
01411       HU = 0.0D0
01412       NQU = 0
01413       CCMAX = 0.3D0
01414       MAXCOR = 3
01415       MSBP = 20
01416       MXNCF = 10
01417 !CInitial call to F.  (LF0 points to YH(*,2).) -------------------------
01418       LF0 = LYH + NYH
01419       CALL F (NEQ, T, Y, RWORK(LF0))
01420       NFE = 1
01421 !CLoad the initial value vector in YH. ---------------------------------
01422       DO 115 I = 1,N
01423  115    RWORK(I+LYH-1) = Y(I)
01424 !CLoad and invert the EWT array.  (H is temporarily set to 1.0.) -------
01425       NQ = 1
01426       H = 1.0D0
01427       CALL DEWSET (N, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT))
01428       DO 120 I = 1,N
01429         IF (RWORK(I+LEWT-1) .LE. 0.0D0) GO TO 621
01430  120    RWORK(I+LEWT-1) = 1.0D0/RWORK(I+LEWT-1)
01431 !C----------------------------------------------------------------------
01432 !CThe coding below computes the step size, H0, to be attempted on the
01433 !Cfirst step, unless the user has supplied a value for this.
01434 !CFirst check that TOUT - T differs significantly from zero.
01435 !CA scalar tolerance quantity TOL is computed, as MAX(RTOL(I))
01436 !Cif this is positive, or MAX(ATOL(I)/ABS(Y(I))) otherwise, adjusted
01437 !Cso as to be between 100*UROUND and 1.0E-3.
01438 !CThen the computed value H0 is given by..
01439 !C                                     NEQ
01440 !C  H0**2 = TOL / ( w0**-2 + (1/NEQ) * SUM ( f(i)/ywt(i) )**2  )
01441 !C                                      1
01442 !Cwhere   w0     = MAX ( ABS(T), ABS(TOUT) ),
01443 !C        f(i)   = i-th component of initial value of f,
01444 !C        ywt(i) = EWT(i)/TOL  (a weight for y(i)).
01445 !CThe sign of H0 is inferred from the initial values of TOUT and T.
01446 !C----------------------------------------------------------------------
01447       IF (H0 .NE. 0.0D0) GO TO 180
01448       TDIST = DABS(TOUT - T)
01449       W0 = DMAX1(DABS(T),DABS(TOUT))
01450       IF (TDIST .LT. 2.0D0*UROUND*W0) GO TO 622
01451       TOL = RTOL(1)
01452       IF (ITOL .LE. 2) GO TO 140
01453       DO 130 I = 1,N
01454  130    TOL = DMAX1(TOL,RTOL(I))
01455  140  IF (TOL .GT. 0.0D0) GO TO 160
01456       ATOLI = ATOL(1)
01457       DO 150 I = 1,N
01458         IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
01459         AYI = DABS(Y(I))
01460         IF (AYI .NE. 0.0D0) TOL = DMAX1(TOL,ATOLI/AYI)
01461  150    CONTINUE
01462  160  TOL = DMAX1(TOL,100.0D0*UROUND)
01463       TOL = DMIN1(TOL,0.001D0)
01464       SUM = DVNORM (N, RWORK(LF0), RWORK(LEWT))
01465       SUM = 1.0D0/(TOL*W0*W0) + TOL*SUM**2
01466       H0 = 1.0D0/DSQRT(SUM)
01467       H0 = DMIN1(H0,TDIST)
01468       H0 = SIGN(H0,TOUT-T)
01469 !CAdjust H0 if necessary to meet HMAX bound. ---------------------------
01470  180  RH = DABS(H0)*HMXI
01471       IF (RH .GT. 1.0D0) H0 = H0/RH
01472 !CLoad H with H0 and scale YH(*,2) by H0. ------------------------------
01473       H = H0
01474       DO 190 I = 1,N
01475  190    RWORK(I+LF0-1) = H0*RWORK(I+LF0-1)
01476       GO TO 270
01477 !C----------------------------------------------------------------------
01478 !CBlock D.
01479 !CThe next code block is for continuation calls only (ISTATE = 2 or 3)
01480 !Cand is to check stop conditions before taking a step.
01481 !C----------------------------------------------------------------------
01482  200  NSLAST = NST
01483       GO TO (210, 250, 220, 230, 240), ITASK
01484  210  IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
01485       CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
01486       IF (IFLAG .NE. 0) GO TO 627
01487       T = TOUT
01488       GO TO 420
01489  220  TP = TN - HU*(1.0D0 + 100.0D0*UROUND)
01490       IF ((TP - TOUT)*H .GT. 0.0D0) GO TO 623
01491       IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
01492       GO TO 400
01493  230  TCRIT = RWORK(1)
01494       IF ((TN - TCRIT)*H .GT. 0.0D0) GO TO 624
01495       IF ((TCRIT - TOUT)*H .LT. 0.0D0) GO TO 625
01496       IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 245
01497       CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
01498       IF (IFLAG .NE. 0) GO TO 627
01499       T = TOUT
01500       GO TO 420
01501  240  TCRIT = RWORK(1)
01502       IF ((TN - TCRIT)*H .GT. 0.0D0) GO TO 624
01503  245  HMX = DABS(TN) + DABS(H)
01504       IHIT = DABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX
01505       IF (IHIT) GO TO 400
01506       TNEXT = TN + H*(1.0D0 + 4.0D0*UROUND)
01507       IF ((TNEXT - TCRIT)*H .LE. 0.0D0) GO TO 250
01508       H = (TCRIT - TN)*(1.0D0 - 4.0D0*UROUND)
01509       IF (ISTATE .EQ. 2) JSTART = -2
01510 !C----------------------------------------------------------------------
01511 !CBlock E.
01512 !CThe next block is normally executed for all calls and contains
01513 !Cthe call to the one-step core integrator DSTODE.
01514 !C
01515 !CThis is a looping point for the integration steps.
01516 !C
01517 !CFirst check for too many steps being taken, update EWT (if not at
01518 !Cstart of problem), check for too much accuracy being requested, and
01519 !Ccheck for H below the roundoff level in T.
01520 !C----------------------------------------------------------------------
01521  250  CONTINUE
01522       IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500
01523       CALL DEWSET (N, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT))
01524       DO 260 I = 1,N
01525         IF (RWORK(I+LEWT-1) .LE. 0.0D0) GO TO 510
01526  260    RWORK(I+LEWT-1) = 1.0D0/RWORK(I+LEWT-1)
01527  270  TOLSF = UROUND*DVNORM (N, RWORK(LYH), RWORK(LEWT))
01528       IF (TOLSF .LE. 1.0D0) GO TO 280
01529       TOLSF = TOLSF*2.0D0
01530       IF (NST .EQ. 0) GO TO 626
01531       GO TO 520
01532  280  IF ((TN + H) .NE. TN) GO TO 290
01533       NHNIL = NHNIL + 1
01534       IF (NHNIL .GT. MXHNIL) GO TO 290
01535       MSG = 'DLSODE-  Warning..internal T (=R1) and H (=R2) are'
01536       CALL XERRWD (MSG, 50, 101, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01537       MSG='      such that in the machine, T + H = T on the next step  '
01538       CALL XERRWD (MSG, 60, 101, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01539       MSG = '      (H = step size). Solver will continue anyway'
01540       CALL XERRWD (MSG, 50, 101, 0, 0, 0, 0, 2, TN, H)
01541       IF (NHNIL .LT. MXHNIL) GO TO 290
01542       MSG = 'DLSODE-  Above warning has been issued I1 times.  '
01543       CALL XERRWD (MSG, 50, 102, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01544       MSG = '      It will not be issued again for this problem'
01545       CALL XERRWD (MSG, 50, 102, 0, 1, MXHNIL, 0, 0, 0.0D0, 0.0D0)
01546  290  CONTINUE
01547 !C----------------------------------------------------------------------
01548 !C CALL DSTODE(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,DPREPJ,DSOLSY)
01549 !C----------------------------------------------------------------------
01550       CALL DSTODE (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), RWORK(LEWT), &
01551          RWORK(LSAVF), RWORK(LACOR), RWORK(LWM), IWORK(LIWM), &
01552          F, JAC, DPREPJ, DSOLSY)
01553       KGO = 1 - KFLAG
01554       GO TO (300, 530, 540), KGO
01555 !C----------------------------------------------------------------------
01556 !CBlock F.
01557 !CThe following block handles the case of a successful return from the
01558 !Ccore integrator (KFLAG = 0).  Test for stop conditions.
01559 !C----------------------------------------------------------------------
01560  300  INIT = 1
01561       GO TO (310, 400, 330, 340, 350), ITASK
01562 !CITASK = 1.  If TOUT has been reached, interpolate. -------------------
01563  310  IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
01564       CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
01565       T = TOUT
01566       GO TO 420
01567 !CITASK = 3.  Jump to exit if TOUT was reached. ------------------------
01568  330  IF ((TN - TOUT)*H .GE. 0.0D0) GO TO 400
01569       GO TO 250
01570 !CITASK = 4.  See if TOUT or TCRIT was reached.  Adjust H if necessary.
01571  340  IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 345
01572       CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
01573       T = TOUT
01574       GO TO 420
01575  345  HMX = DABS(TN) + ABS(H)
01576       IHIT = DABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX
01577       IF (IHIT) GO TO 400
01578       TNEXT = TN + H*(1.0D0 + 4.0D0*UROUND)
01579       IF ((TNEXT - TCRIT)*H .LE. 0.0D0) GO TO 250
01580       H = (TCRIT - TN)*(1.0D0 - 4.0D0*UROUND)
01581       JSTART = -2
01582       GO TO 250
01583 !CITASK = 5.  See if TCRIT was reached and jump to exit. ---------------
01584  350  HMX = DABS(TN) + ABS(H)
01585       IHIT = DABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX
01586 !C----------------------------------------------------------------------
01587 !CBlock G.
01588 !CThe following block handles all successful returns from DLSODE.
01589 !CIf ITASK .NE. 1, Y is loaded from YH and T is set accordingly.
01590 !CISTATE is set to 2, and the optional outputs are loaded into the
01591 !Cwork arrays before returning.
01592 !C----------------------------------------------------------------------
01593  400  DO 410 I = 1,N
01594  410    Y(I) = RWORK(I+LYH-1)
01595       T = TN
01596       IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 420
01597       IF (IHIT) T = TCRIT
01598  420  ISTATE = 2
01599       RWORK(11) = HU
01600       RWORK(12) = H
01601       RWORK(13) = TN
01602       IWORK(11) = NST
01603       IWORK(12) = NFE
01604       IWORK(13) = NJE
01605       IWORK(14) = NQU
01606       IWORK(15) = NQ
01607       RETURN
01608 !C----------------------------------------------------------------------
01609 !CBlock H.
01610 !CThe following block handles all unsuccessful returns other than
01611 !Cthose for illegal input.  First the error message routine is called.
01612 !CIf there was an error test or convergence test failure, IMXER is set.
01613 !CThen Y is loaded from YH and T is set to TN.  The optional outputs
01614 !Care loaded into the work arrays before returning.
01615 !C----------------------------------------------------------------------
01616 !CThe maximum number of steps was taken before reaching TOUT. ----------
01617  500  MSG = 'DLSODE-  At current T (=R1), MXSTEP (=I1) steps   '
01618       CALL XERRWD (MSG, 50, 201, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01619       MSG = '      taken on this call before reaching TOUT     '
01620       CALL XERRWD (MSG, 50, 201, 0, 1, MXSTEP, 0, 1, TN, 0.0D0)
01621       ISTATE = -1
01622       GO TO 580
01623 !CEWT(I) .LE. 0.0 for some I (not at start of problem). ----------------
01624  510  EWTI = RWORK(LEWT+I-1)
01625       MSG = 'DLSODE-  At T (=R1), EWT(I1) has become R2 .LE. 0.'
01626       CALL XERRWD (MSG, 50, 202, 0, 1, I, 0, 2, TN, EWTI)
01627       ISTATE = -6
01628       GO TO 580
01629 !CToo much accuracy requested for machine precision. -------------------
01630  520  MSG = 'DLSODE-  At T (=R1), too much accuracy requested  '
01631       CALL XERRWD (MSG, 50, 203, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01632       MSG = '      for precision of machine..  see TOLSF (=R2) '
01633       CALL XERRWD (MSG, 50, 203, 0, 0, 0, 0, 2, TN, TOLSF)
01634       RWORK(14) = TOLSF
01635       ISTATE = -2
01636       GO TO 580
01637 !CKFLAG = -1.  Error test failed repeatedly or with ABS(H) = HMIN. -----
01638  530  MSG = 'DLSODE-  At T(=R1) and step size H(=R2), the error'
01639       CALL XERRWD (MSG, 50, 204, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01640       MSG = '      test failed repeatedly or with ABS(H) = HMIN'
01641       CALL XERRWD (MSG, 50, 204, 0, 0, 0, 0, 2, TN, H)
01642       ISTATE = -4
01643       GO TO 560
01644 !CKFLAG = -2.  Convergence failed repeatedly or with ABS(H) = HMIN. ----
01645  540  MSG = 'DLSODE-  At T (=R1) and step size H (=R2), the    '
01646       CALL XERRWD (MSG, 50, 205, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01647       MSG = '      corrector convergence failed repeatedly     '
01648       CALL XERRWD (MSG, 50, 205, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01649       MSG = '      or with ABS(H) = HMIN   '
01650       CALL XERRWD (MSG, 30, 205, 0, 0, 0, 0, 2, TN, H)
01651       ISTATE = -5
01652 !CCompute IMXER if relevant. -------------------------------------------
01653  560  BIG = 0.0D0
01654       IMXER = 1
01655       DO 570 I = 1,N
01656         SIZE = DABS(RWORK(I+LACOR-1)*RWORK(I+LEWT-1))
01657         IF (BIG .GE. SIZE) GO TO 570
01658         BIG = SIZE
01659         IMXER = I
01660  570    CONTINUE
01661       IWORK(16) = IMXER
01662 !CSet Y vector, T, and optional outputs. -------------------------------
01663  580  DO 590 I = 1,N
01664  590    Y(I) = RWORK(I+LYH-1)
01665       T = TN
01666       RWORK(11) = HU
01667       RWORK(12) = H
01668       RWORK(13) = TN
01669       IWORK(11) = NST
01670       IWORK(12) = NFE
01671       IWORK(13) = NJE
01672       IWORK(14) = NQU
01673       IWORK(15) = NQ
01674       RETURN
01675 !C----------------------------------------------------------------------
01676 !CBlock I.
01677 !CThe following block handles all error returns due to illegal input
01678 !C(ISTATE = -3), as detected before calling the core integrator.
01679 !CFirst the error message routine is called.  If the illegal input
01680 !Cis a negative ISTATE, the run is aborted (apparent infinite loop).
01681 !C----------------------------------------------------------------------
01682  601  MSG = 'DLSODE-  ISTATE (=I1) illegal '
01683       CALL XERRWD (MSG, 30, 1, 0, 1, ISTATE, 0, 0, 0.0D0, 0.0D0)
01684       IF (ISTATE .LT. 0) GO TO 800
01685       GO TO 700
01686  602  MSG = 'DLSODE-  ITASK (=I1) illegal  '
01687       CALL XERRWD (MSG, 30, 2, 0, 1, ITASK, 0, 0, 0.0D0, 0.0D0)
01688       GO TO 700
01689  603  MSG = 'DLSODE-  ISTATE .GT. 1 but DLSODE not initialized '
01690       CALL XERRWD (MSG, 50, 3, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01691       GO TO 700
01692  604  MSG = 'DLSODE-  NEQ (=I1) .LT. 1     '
01693       CALL XERRWD (MSG, 30, 4, 0, 1, NEQ(1), 0, 0, 0.0D0, 0.0D0)
01694       GO TO 700
01695  605  MSG = 'DLSODE-  ISTATE = 3 and NEQ increased (I1 to I2)  '
01696       CALL XERRWD (MSG, 50, 5, 0, 2, N, NEQ(1), 0, 0.0D0, 0.0D0)
01697       GO TO 700
01698  606  MSG = 'DLSODE-  ITOL (=I1) illegal   '
01699       CALL XERRWD (MSG, 30, 6, 0, 1, ITOL, 0, 0, 0.0D0, 0.0D0)
01700       GO TO 700
01701  607  MSG = 'DLSODE-  IOPT (=I1) illegal   '
01702       CALL XERRWD (MSG, 30, 7, 0, 1, IOPT, 0, 0, 0.0D0, 0.0D0)
01703       GO TO 700
01704  608  MSG = 'DLSODE-  MF (=I1) illegal     '
01705       CALL XERRWD (MSG, 30, 8, 0, 1, MF, 0, 0, 0.0D0, 0.0D0)
01706       GO TO 700
01707  609  MSG = 'DLSODE-  ML (=I1) illegal.. .LT.0 or .GE.NEQ (=I2)'
01708       CALL XERRWD (MSG, 50, 9, 0, 2, ML, NEQ(1), 0, 0.0D0, 0.0D0)
01709       GO TO 700
01710  610  MSG = 'DLSODE-  MU (=I1) illegal.. .LT.0 or .GE.NEQ (=I2)'
01711       CALL XERRWD (MSG, 50, 10, 0, 2, MU, NEQ(1), 0, 0.0D0, 0.0D0)
01712       GO TO 700
01713  611  MSG = 'DLSODE-  MAXORD (=I1) .LT. 0  '
01714       CALL XERRWD (MSG, 30, 11, 0, 1, MAXORD, 0, 0, 0.0D0, 0.0D0)
01715       GO TO 700
01716  612  MSG = 'DLSODE-  MXSTEP (=I1) .LT. 0  '
01717       CALL XERRWD (MSG, 30, 12, 0, 1, MXSTEP, 0, 0, 0.0D0, 0.0D0)
01718       GO TO 700
01719  613  MSG = 'DLSODE-  MXHNIL (=I1) .LT. 0  '
01720       CALL XERRWD (MSG, 30, 13, 0, 1, MXHNIL, 0, 0, 0.0D0, 0.0D0)
01721       GO TO 700
01722  614  MSG = 'DLSODE-  TOUT (=R1) behind T (=R2)      '
01723       CALL XERRWD (MSG, 40, 14, 0, 0, 0, 0, 2, TOUT, T)
01724       MSG = '      Integration direction is given by H0 (=R1)  '
01725       CALL XERRWD (MSG, 50, 14, 0, 0, 0, 0, 1, H0, 0.0D0)
01726       GO TO 700
01727  615  MSG = 'DLSODE-  HMAX (=R1) .LT. 0.0  '
01728       CALL XERRWD (MSG, 30, 15, 0, 0, 0, 0, 1, HMAX, 0.0D0)
01729       GO TO 700
01730  616  MSG = 'DLSODE-  HMIN (=R1) .LT. 0.0  '
01731       CALL XERRWD (MSG, 30, 16, 0, 0, 0, 0, 1, HMIN, 0.0D0)
01732       GO TO 700
01733  617  CONTINUE
01734       MSG='DLSODE-  RWORK length needed, LENRW (=I1), exceeds LRW (=I2)'
01735       CALL XERRWD (MSG, 60, 17, 0, 2, LENRW, LRW, 0, 0.0D0, 0.0D0)
01736       GO TO 700
01737  618  CONTINUE
01738       MSG='DLSODE-  IWORK length needed, LENIW (=I1), exceeds LIW (=I2)'
01739       CALL XERRWD (MSG, 60, 18, 0, 2, LENIW, LIW, 0, 0.0D0, 0.0D0)
01740       GO TO 700
01741  619  MSG = 'DLSODE-  RTOL(I1) is R1 .LT. 0.0        '
01742       CALL XERRWD (MSG, 40, 19, 0, 1, I, 0, 1, RTOLI, 0.0D0)
01743       GO TO 700
01744  620  MSG = 'DLSODE-  ATOL(I1) is R1 .LT. 0.0        '
01745       CALL XERRWD (MSG, 40, 20, 0, 1, I, 0, 1, ATOLI, 0.0D0)
01746       GO TO 700
01747  621  EWTI = RWORK(LEWT+I-1)
01748       MSG = 'DLSODE-  EWT(I1) is R1 .LE. 0.0         '
01749       CALL XERRWD (MSG, 40, 21, 0, 1, I, 0, 1, EWTI, 0.0D0)
01750       GO TO 700
01751  622  CONTINUE
01752       MSG='DLSODE-  TOUT (=R1) too close to T(=R2) to start integration'
01753       CALL XERRWD (MSG, 60, 22, 0, 0, 0, 0, 2, TOUT, T)
01754       GO TO 700
01755  623  CONTINUE
01756       MSG='DLSODE-  ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2)  '
01757       CALL XERRWD (MSG, 60, 23, 0, 1, ITASK, 0, 2, TOUT, TP)
01758       GO TO 700
01759  624  CONTINUE
01760       MSG='DLSODE-  ITASK = 4 OR 5 and TCRIT (=R1) behind TCUR (=R2)   '
01761       CALL XERRWD (MSG, 60, 24, 0, 0, 0, 0, 2, TCRIT, TN)
01762       GO TO 700
01763  625  CONTINUE
01764       MSG='DLSODE-  ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2)   '
01765       CALL XERRWD (MSG, 60, 25, 0, 0, 0, 0, 2, TCRIT, TOUT)
01766       GO TO 700
01767  626  MSG = 'DLSODE-  At start of problem, too much accuracy   '
01768       CALL XERRWD (MSG, 50, 26, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01769       MSG='      requested for precision of machine..  See TOLSF (=R1) '
01770       CALL XERRWD (MSG, 60, 26, 0, 0, 0, 0, 1, TOLSF, 0.0D0)
01771       RWORK(14) = TOLSF
01772       GO TO 700
01773  627  MSG = 'DLSODE-  Trouble in DINTDY.  ITASK = I1, TOUT = R1'
01774       CALL XERRWD (MSG, 50, 27, 0, 1, ITASK, 0, 1, TOUT, 0.0D0)
01775 !C
01776  700  ISTATE = -3
01777       RETURN
01778 !C
01779  800  MSG = 'DLSODE-  Run aborted.. apparent infinite loop     '
01780       CALL XERRWD (MSG, 50, 303, 2, 0, 0, 0, 0, 0.0D0, 0.0D0)
01781       RETURN
01782 !C---------------------- END OF SUBROUTINE DLSODE ----------------------
01783       END
01784 !*DECK DUMACH
01785       DOUBLE PRECISION FUNCTION DUMACH ()
01786 !C**BEGIN PROLOGUE  DUMACH
01787 !C**PURPOSE  Compute the unit roundoff of the machine.
01788 !C**CATEGORY  R1
01789 !C**TYPE      DOUBLE PRECISION (RUMACH-S, DUMACH-D)
01790 !C**KEYWORDS  MACHINE CONSTANTS
01791 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
01792 !C**DESCRIPTION
01793 !C*Usage:
01794 !C       DOUBLE PRECISION  A, DUMACH
01795 !C       A = DUMACH()
01796 !C
01797 !C*Function Return Values:
01798 !C    A : the unit roundoff of the machine.
01799 !C
01800 !C*Description:
01801 !C    The unit roundoff is defined as the smallest positive machine
01802 !C    number u such that  1.0 + u .ne. 1.0.  This is computed by DUMACH
01803 !C    in a machine-independent manner.
01804 !C
01805 !C**REFERENCES  (NONE)
01806 !C**ROUTINES CALLED  DUMSUM
01807 !C**REVISION HISTORY  (YYYYMMDD)
01808 !C  19930216  DATE WRITTEN
01809 !C  19930818  Added SLATEC-format prologue.  (FNF)
01810 !C  20030707  Added DUMSUM to force normal storage of COMP.  (ACH)
01811 !C**END PROLOGUE  DUMACH
01812 !C
01813       DOUBLE PRECISION U, COMP
01814 !C**FIRST EXECUTABLE STATEMENT  DUMACH
01815       U = 1.0D0
01816  10   U = U*0.5D0
01817       CALL DUMSUM(1.0D0, U, COMP)
01818       IF (COMP .NE. 1.0D0) GO TO 10
01819       DUMACH = U*2.0D0
01820       RETURN
01821 !C---------------------- End of Function DUMACH ------------------------
01822       END
01823       SUBROUTINE DUMSUM(A,B,C)
01824 !C    Routine to force normal storing of A + B, for DUMACH.
01825       DOUBLE PRECISION A, B, C
01826       C = A + B
01827       RETURN
01828       END
01829 !*DECK DCFODE
01830       SUBROUTINE DCFODE (METH, ELCO, TESCO)
01831 !C**BEGIN PROLOGUE  DCFODE
01832 !C**SUBSIDIARY
01833 !C**PURPOSE  Set ODE integrator coefficients.
01834 !C**TYPE      DOUBLE PRECISION (SCFODE-S, DCFODE-D)
01835 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
01836 !C**DESCRIPTION
01837 !C
01838 !C DCFODE is called by the integrator routine to set coefficients
01839 !C needed there.  The coefficients for the current method, as
01840 !C given by the value of METH, are set for all orders and saved.
01841 !C The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2.
01842 !C (A smaller value of the maximum order is also allowed.)
01843 !C DCFODE is called once at the beginning of the problem,
01844 !C and is not called again unless and until METH is changed.
01845 !C
01846 !C The ELCO array contains the basic method coefficients.
01847 !C The coefficients el(i), 1 .le. i .le. nq+1, for the method of
01848 !C order nq are stored in ELCO(i,nq).  They are given by a genetrating
01849 !C polynomial, i.e.,
01850 !C     l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq.
01851 !C For the implicit Adams methods, l(x) is given by
01852 !C     dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1),    l(-1) = 0.
01853 !C For the BDF methods, l(x) is given by
01854 !C     l(x) = (x+1)*(x+2)* ... *(x+nq)/K,
01855 !C where         K = factorial(nq)*(1 + 1/2 + ... + 1/nq).
01856 !C
01857 !C The TESCO array contains test constants used for the
01858 !C local error test and the selection of step size and/or order.
01859 !C At order nq, TESCO(k,nq) is used for the selection of step
01860 !C size at order nq - 1 if k = 1, at order nq if k = 2, and at order
01861 !C nq + 1 if k = 3.
01862 !C
01863 !C**SEE ALSO  DLSODE
01864 !C**ROUTINES CALLED  (NONE)
01865 !C**REVISION HISTORY  (YYMMDD)
01866 !C  791129  DATE WRITTEN
01867 !C  890501  Modified prologue to SLATEC/LDOC format.  (FNF)
01868 !C  890503  Minor cosmetic changes.  (FNF)
01869 !C  930809  Renamed to allow single/double precision versions. (ACH)
01870 !C**END PROLOGUE  DCFODE
01871 !C*End
01872       INTEGER METH
01873       INTEGER I, IB, NQ, NQM1, NQP1
01874       DOUBLE PRECISION ELCO, TESCO
01875       DOUBLE PRECISION AGAMQ, FNQ, FNQM1, PC, PINT, RAGQ, 
01876          RQFAC, RQ1FAC, TSIGN, XPIN
01877       DIMENSION ELCO(13,12), TESCO(3,12)
01878       DIMENSION PC(12)
01879 !C
01880 !C**FIRST EXECUTABLE STATEMENT  DCFODE
01881       GO TO (100, 200), METH
01882 !C
01883  100  ELCO(1,1) = 1.0D0
01884       ELCO(2,1) = 1.0D0
01885       TESCO(1,1) = 0.0D0
01886       TESCO(2,1) = 2.0D0
01887       TESCO(1,2) = 1.0D0
01888       TESCO(3,12) = 0.0D0
01889       PC(1) = 1.0D0
01890       RQFAC = 1.0D0
01891       DO 140 NQ = 2,12
01892 !C----------------------------------------------------------------------
01893 !CThe PC array will contain the coefficients of the polynomial
01894 !C    p(x) = (x+1)*(x+2)*...*(x+nq-1).
01895 !CInitially, p(x) = 1.
01896 !C----------------------------------------------------------------------
01897         RQ1FAC = RQFAC
01898         RQFAC = RQFAC/NQ
01899         NQM1 = NQ - 1
01900         FNQM1 = NQM1
01901         NQP1 = NQ + 1
01902 !CForm coefficients of p(x)*(x+nq-1). ----------------------------------
01903         PC(NQ) = 0.0D0
01904         DO 110 IB = 1,NQM1
01905           I = NQP1 - IB
01906  110      PC(I) = PC(I-1) + FNQM1*PC(I)
01907         PC(1) = FNQM1*PC(1)
01908 !CCompute integral, -1 to 0, of p(x) and x*p(x). -----------------------
01909         PINT = PC(1)
01910         XPIN = PC(1)/2.0D0
01911         TSIGN = 1.0D0
01912         DO 120 I = 2,NQ
01913           TSIGN = -TSIGN
01914           PINT = PINT + TSIGN*PC(I)/I
01915  120      XPIN = XPIN + TSIGN*PC(I)/(I+1)
01916 !CStore coefficients in ELCO and TESCO. --------------------------------
01917         ELCO(1,NQ) = PINT*RQ1FAC
01918         ELCO(2,NQ) = 1.0D0
01919         DO 130 I = 2,NQ
01920  130      ELCO(I+1,NQ) = RQ1FAC*PC(I)/I
01921         AGAMQ = RQFAC*XPIN
01922         RAGQ = 1.0D0/AGAMQ
01923         TESCO(2,NQ) = RAGQ
01924         IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/NQP1
01925         TESCO(3,NQM1) = RAGQ
01926  140    CONTINUE
01927       RETURN
01928 !C
01929  200  PC(1) = 1.0D0
01930       RQ1FAC = 1.0D0
01931       DO 230 NQ = 1,5
01932 !C----------------------------------------------------------------------
01933 !CThe PC array will contain the coefficients of the polynomial
01934 !C    p(x) = (x+1)*(x+2)*...*(x+nq).
01935 !CInitially, p(x) = 1.
01936 !C----------------------------------------------------------------------
01937         FNQ = NQ
01938         NQP1 = NQ + 1
01939 !CForm coefficients of p(x)*(x+nq). ------------------------------------
01940         PC(NQP1) = 0.0D0
01941         DO 210 IB = 1,NQ
01942           I = NQ + 2 - IB
01943  210      PC(I) = PC(I-1) + FNQ*PC(I)
01944         PC(1) = FNQ*PC(1)
01945 !CStore coefficients in ELCO and TESCO. --------------------------------
01946         DO 220 I = 1,NQP1
01947  220      ELCO(I,NQ) = PC(I)/PC(2)
01948         ELCO(2,NQ) = 1.0D0
01949         TESCO(1,NQ) = RQ1FAC
01950         TESCO(2,NQ) = NQP1/ELCO(1,NQ)
01951         TESCO(3,NQ) = (NQ+2)/ELCO(1,NQ)
01952         RQ1FAC = RQ1FAC/FNQ
01953  230    CONTINUE
01954       RETURN
01955 !C---------------------- END OF SUBROUTINE DCFODE ----------------------
01956       END
01957 !*DECK DINTDY
01958       SUBROUTINE DINTDY (T, K, YH, NYH, DKY, IFLAG)
01959 !C**BEGIN PROLOGUE  DINTDY
01960 !C**SUBSIDIARY
01961 !C**PURPOSE  Interpolate solution derivatives.
01962 !C**TYPE      DOUBLE PRECISION (SINTDY-S, DINTDY-D)
01963 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
01964 !C**DESCRIPTION
01965 !C
01966 !C DINTDY computes interpolated values of the K-th derivative of the
01967 !C dependent variable vector y, and stores it in DKY.  This routine
01968 !C is called within the package with K = 0 and T = TOUT, but may
01969 !C also be called by the user for any K up to the current order.
01970 !C (See detailed instructions in the usage documentation.)
01971 !C
01972 !C The computed values in DKY are gotten by interpolation using the
01973 !C Nordsieck history array YH.  This array corresponds uniquely to a
01974 !C vector-valued polynomial of degree NQCUR or less, and DKY is set
01975 !C to the K-th derivative of this polynomial at T.
01976 !C The formula for DKY is:
01977 !C              q
01978 !C  DKY(i)  =  sum  c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1)
01979 !C             j=K
01980 !C where  c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR.
01981 !C The quantities  nq = NQCUR, l = nq+1, N = NEQ, tn, and h are
01982 !C communicated by COMMON.  The above sum is done in reverse order.
01983 !C IFLAG is returned negative if either K or T is out of bounds.
01984 !C
01985 !C**SEE ALSO  DLSODE
01986 !C**ROUTINES CALLED  XERRWD
01987 !C**COMMON BLOCKS    DLS001
01988 !C**REVISION HISTORY  (YYMMDD)
01989 !C  791129  DATE WRITTEN
01990 !C  890501  Modified prologue to SLATEC/LDOC format.  (FNF)
01991 !C  890503  Minor cosmetic changes.  (FNF)
01992 !C  930809  Renamed to allow single/double precision versions. (ACH)
01993 !C  010418  Reduced size of Common block /DLS001/. (ACH)
01994 !C  031105  Restored 'own' variables to Common block /DLS001/, to
01995 !C          enable interrupt/restart feature. (ACH)
01996 !C  050427  Corrected roundoff decrement in TP. (ACH)
01997 !C**END PROLOGUE  DINTDY
01998 !C*End
01999       INTEGER K, NYH, IFLAG
02000       DOUBLE PRECISION T, YH, DKY
02001       DIMENSION YH(NYH,*), DKY(*)
02002       INTEGER IOWND, IOWNS, 
02003          ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 
02004          LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 
02005          MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
02006       DOUBLE PRECISION ROWNS, 
02007         CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
02008       COMMON /DLS001/ ROWNS(209), &
02009         CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
02010         IOWND(6), IOWNS(6), &
02011         ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
02012         LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
02013         MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
02014       INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
02015       DOUBLE PRECISION C, R, S, TP
02016       CHARACTER*80 MSG
02017 !C
02018 !C**FIRST EXECUTABLE STATEMENT  DINTDY
02019       IFLAG = 0
02020       IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80
02021       TP = TN - HU -  100.0D0*UROUND*SIGN(DABS(TN) + DABS(HU), HU)
02022       IF ((T-TP)*(T-TN) .GT. 0.0D0) GO TO 90
02023 !C
02024       S = (T - TN)/H
02025       IC = 1
02026       IF (K .EQ. 0) GO TO 15
02027       JJ1 = L - K
02028       DO 10 JJ = JJ1,NQ
02029  10     IC = IC*JJ
02030  15   C = IC
02031       DO 20 I = 1,N
02032  20     DKY(I) = C*YH(I,L)
02033       IF (K .EQ. NQ) GO TO 55
02034       JB2 = NQ - K
02035       DO 50 JB = 1,JB2
02036         J = NQ - JB
02037         JP1 = J + 1
02038         IC = 1
02039         IF (K .EQ. 0) GO TO 35
02040         JJ1 = JP1 - K
02041         DO 30 JJ = JJ1,J
02042  30       IC = IC*JJ
02043  35     C = IC
02044         DO 40 I = 1,N
02045  40       DKY(I) = C*YH(I,JP1) + S*DKY(I)
02046  50     CONTINUE
02047       IF (K .EQ. 0) RETURN
02048  55   R = H**(-K)
02049       DO 60 I = 1,N
02050  60     DKY(I) = R*DKY(I)
02051       RETURN
02052 !C
02053  80   MSG = 'DINTDY-  K (=I1) illegal      '
02054       CALL XERRWD (MSG, 30, 51, 0, 1, K, 0, 0, 0.0D0, 0.0D0)
02055       IFLAG = -1
02056       RETURN
02057  90   MSG = 'DINTDY-  T (=R1) illegal      '
02058       CALL XERRWD (MSG, 30, 52, 0, 0, 0, 0, 1, T, 0.0D0)
02059       MSG='      T not in interval TCUR - HU (= R1) to TCUR (=R2)      '
02060       CALL XERRWD (MSG, 60, 52, 0, 0, 0, 0, 2, TP, TN)
02061       IFLAG = -2
02062       RETURN
02063 !C---------------------- END OF SUBROUTINE DINTDY ----------------------
02064       END
02065 !*DECK DPREPJ
02066       SUBROUTINE DPREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, &
02067          F, JAC)
02068 !C**BEGIN PROLOGUE  DPREPJ
02069 !C**SUBSIDIARY
02070 !C**PURPOSE  Compute and process Newton iteration matrix.
02071 !C**TYPE      DOUBLE PRECISION (SPREPJ-S, DPREPJ-D)
02072 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
02073 !C**DESCRIPTION
02074 !C
02075 !C DPREPJ is called by DSTODE to compute and process the matrix
02076 !C P = I - h*el(1)*J , where J is an approximation to the Jacobian.
02077 !C Here J is computed by the user-supplied routine JAC if
02078 !C MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5.
02079 !C If MITER = 3, a diagonal approximation to J is used.
02080 !C J is stored in WM and replaced by P.  If MITER .ne. 3, P is then
02081 !C subjected to LU decomposition in preparation for later solution
02082 !C of linear systems with P as coefficient matrix.  This is done
02083 !C by DGEFA if MITER = 1 or 2, and by DGBFA if MITER = 4 or 5.
02084 !C
02085 !C In addition to variables described in DSTODE and DLSODE prologues,
02086 !C communication with DPREPJ uses the following:
02087 !C Y     = array containing predicted values on entry.
02088 !C FTEM  = work array of length N (ACOR in DSTODE).
02089 !C SAVF  = array containing f evaluated at predicted y.
02090 !C WM    = real work space for matrices.  On output it contains the
02091 !C         inverse diagonal matrix if MITER = 3 and the LU decomposition
02092 !C         of P if MITER is 1, 2 , 4, or 5.
02093 !C         Storage of matrix elements starts at WM(3).
02094 !C         WM also contains the following matrix-related data:
02095 !C         WM(1) = SQRT(UROUND), used in numerical Jacobian increments.
02096 !C         WM(2) = H*EL0, saved for later use if MITER = 3.
02097 !C IWM   = integer work space containing pivot information, starting at
02098 !C         IWM(21), if MITER is 1, 2, 4, or 5.  IWM also contains band
02099 !C         parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
02100 !C EL0   = EL(1) (input).
02101 !C IERPJ = output error flag,  = 0 if no trouble, .gt. 0 if
02102 !C         P matrix found to be singular.
02103 !C JCUR  = output flag = 1 to indicate that the Jacobian matrix
02104 !C         (or approximation) is now current.
02105 !C This routine also uses the COMMON variables EL0, H, TN, UROUND,
02106 !C MITER, N, NFE, and NJE.
02107 !C
02108 !C**SEE ALSO  DLSODE
02109 !C**ROUTINES CALLED  DGBFA, DGEFA, DVNORM
02110 !C**COMMON BLOCKS    DLS001
02111 !C**REVISION HISTORY  (YYMMDD)
02112 !C  791129  DATE WRITTEN
02113 !C  890501  Modified prologue to SLATEC/LDOC format.  (FNF)
02114 !C  890504  Minor cosmetic changes.  (FNF)
02115 !C  930809  Renamed to allow single/double precision versions. (ACH)
02116 !C  010418  Reduced size of Common block /DLS001/. (ACH)
02117 !C  031105  Restored 'own' variables to Common block /DLS001/, to
02118 !C          enable interrupt/restart feature. (ACH)
02119 !C**END PROLOGUE  DPREPJ
02120 !C*End
02121       EXTERNAL F, JAC
02122       INTEGER NEQ, NYH, IWM
02123       DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WM
02124       DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*), &
02125          WM(*), IWM(*)
02126       INTEGER IOWND, IOWNS, 
02127          ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 
02128          LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 
02129          MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
02130       DOUBLE PRECISION ROWNS, 
02131          CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
02132       COMMON /DLS001/ ROWNS(209), &
02133          CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
02134          IOWND(6), IOWNS(6), &
02135          ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
02136          LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
02137         MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
02138       INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP, 
02139          MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1
02140       DOUBLE PRECISION CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ, 
02141          DVNORM
02142 !C
02143 !C**FIRST EXECUTABLE STATEMENT  DPREPJ
02144       NJE = NJE + 1
02145       IERPJ = 0
02146       JCUR = 1
02147       HL0 = H*EL0
02148       GO TO (100, 200, 300, 400, 500), MITER
02149 !CIf MITER = 1, call JAC and multiply by scalar. -----------------------
02150  100  LENP = N*N
02151       DO 110 I = 1,LENP
02152  110    WM(I+2) = 0.0D0
02153       CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N)
02154       CON = -HL0
02155       DO 120 I = 1,LENP
02156  120    WM(I+2) = WM(I+2)*CON
02157       GO TO 240
02158 !CIf MITER = 2, make N calls to F to approximate J. --------------------
02159  200  FAC = DVNORM (N, SAVF, EWT)
02160       R0 = 1000.0D0*DABS(H)*UROUND*N*FAC
02161       IF (R0 .EQ. 0.0D0) R0 = 1.0D0
02162       SRUR = WM(1)
02163       J1 = 2
02164       DO 230 J = 1,N
02165         YJ = Y(J)
02166         R = DMAX1(SRUR*DABS(YJ),R0/EWT(J))
02167         Y(J) = Y(J) + R
02168         FAC = -HL0/R
02169         CALL F (NEQ, TN, Y, FTEM)
02170         DO 220 I = 1,N
02171  220      WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
02172         Y(J) = YJ
02173         J1 = J1 + N
02174  230    CONTINUE
02175       NFE = NFE + N
02176 !CAdd identity matrix. -------------------------------------------------
02177  240  J = 3
02178       NP1 = N + 1
02179       DO 250 I = 1,N
02180         WM(J) = WM(J) + 1.0D0
02181  250    J = J + NP1
02182 !CDo LU decomposition on P. --------------------------------------------
02183       CALL DGEFA (WM(3), N, N, IWM(21), IER)
02184       IF (IER .NE. 0) IERPJ = 1
02185       RETURN
02186 !CIf MITER = 3, construct a diagonal approximation to J and P. ---------
02187  300  WM(2) = HL0
02188       R = EL0*0.1D0
02189       DO 310 I = 1,N
02190  310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
02191       CALL F (NEQ, TN, Y, WM(3))
02192       NFE = NFE + 1
02193       DO 320 I = 1,N
02194         R0 = H*SAVF(I) - YH(I,2)
02195         DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I))
02196         WM(I+2) = 1.0D0
02197         IF (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320
02198         IF (DABS(DI) .EQ. 0.0D0) GO TO 330
02199         WM(I+2) = 0.1D0*R0/DI
02200  320    CONTINUE
02201       RETURN
02202  330  IERPJ = 1
02203       RETURN
02204 !CIf MITER = 4, call JAC and multiply by scalar. -----------------------
02205  400  ML = IWM(1)
02206       MU = IWM(2)
02207       ML3 = ML + 3
02208       MBAND = ML + MU + 1
02209       MEBAND = MBAND + ML
02210       LENP = MEBAND*N
02211       DO 410 I = 1,LENP
02212  410    WM(I+2) = 0.0D0
02213       CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND)
02214       CON = -HL0
02215       DO 420 I = 1,LENP
02216  420    WM(I+2) = WM(I+2)*CON
02217       GO TO 570
02218 !CIf MITER = 5, make MBAND calls to F to approximate J. ----------------
02219  500  ML = IWM(1)
02220       MU = IWM(2)
02221       MBAND = ML + MU + 1
02222       MBA = MIN(MBAND,N)
02223       MEBAND = MBAND + ML
02224       MEB1 = MEBAND - 1
02225       SRUR = WM(1)
02226       FAC = DVNORM (N, SAVF, EWT)
02227       R0 = 1000.0D0*DABS(H)*UROUND*N*FAC
02228       IF (R0 .EQ. 0.0D0) R0 = 1.0D0
02229       DO 560 J = 1,MBA
02230         DO 530 I = J,N,MBAND
02231           YI = Y(I)
02232           R = DMAX1(SRUR*DABS(YI),R0/EWT(I))
02233  530      Y(I) = Y(I) + R
02234         CALL F (NEQ, TN, Y, FTEM)
02235         DO 550 JJ = J,N,MBAND
02236           Y(JJ) = YH(JJ,1)
02237           YJJ = Y(JJ)
02238           R = DMAX1(SRUR*DABS(YJJ),R0/EWT(JJ))
02239           FAC = -HL0/R
02240           I1 = MAX(JJ-MU,1)
02241           I2 = MIN(JJ+ML,N)
02242           II = JJ*MEB1 - ML + 2
02243           DO 540 I = I1,I2
02244  540        WM(II+I) = (FTEM(I) - SAVF(I))*FAC
02245  550      CONTINUE
02246  560    CONTINUE
02247       NFE = NFE + MBA
02248 !CAdd identity matrix. -------------------------------------------------
02249  570  II = MBAND + 2
02250       DO 580 I = 1,N
02251         WM(II) = WM(II) + 1.0D0
02252  580    II = II + MEBAND
02253 !CDo LU decomposition of P. --------------------------------------------
02254       CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
02255       IF (IER .NE. 0) IERPJ = 1
02256       RETURN
02257 !C---------------------- END OF SUBROUTINE DPREPJ ----------------------
02258       END
02259 !*DECK DSOLSY
02260       SUBROUTINE DSOLSY (WM, IWM, X, TEM)
02261 !C**BEGIN PROLOGUE  DSOLSY
02262 !C**SUBSIDIARY
02263 !C**PURPOSE  ODEPACK linear system solver.
02264 !C**TYPE      DOUBLE PRECISION (SSOLSY-S, DSOLSY-D)
02265 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
02266 !C**DESCRIPTION
02267 !C
02268 !C This routine manages the solution of the linear system arising from
02269 !C a chord iteration.  It is called if MITER .ne. 0.
02270 !C If MITER is 1 or 2, it calls DGESL to accomplish this.
02271 !C If MITER = 3 it updates the coefficient h*EL0 in the diagonal
02272 !C matrix, and then computes the solution.
02273 !C If MITER is 4 or 5, it calls DGBSL.
02274 !C Communication with DSOLSY uses the following variables:
02275 !C WM    = real work space containing the inverse diagonal matrix if
02276 !C         MITER = 3 and the LU decomposition of the matrix otherwise.
02277 !C         Storage of matrix elements starts at WM(3).
02278 !C         WM also contains the following matrix-related data:
02279 !C         WM(1) = SQRT(UROUND) (not used here),
02280 !C         WM(2) = HL0, the previous value of h*EL0, used if MITER = 3.
02281 !C IWM   = integer work space containing pivot information, starting at
02282 !C         IWM(21), if MITER is 1, 2, 4, or 5.  IWM also contains band
02283 !C         parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
02284 !C X     = the right-hand side vector on input, and the solution vector
02285 !C         on output, of length N.
02286 !C TEM   = vector of work space of length N, not used in this version.
02287 !C IERSL = output flag (in COMMON).  IERSL = 0 if no trouble occurred.
02288 !C         IERSL = 1 if a singular matrix arose with MITER = 3.
02289 !C This routine also uses the COMMON variables EL0, H, MITER, and N.
02290 !C
02291 !C**SEE ALSO  DLSODE
02292 !C**ROUTINES CALLED  DGBSL, DGESL
02293 !C**COMMON BLOCKS    DLS001
02294 !C**REVISION HISTORY  (YYMMDD)
02295 !C  791129  DATE WRITTEN
02296 !C  890501  Modified prologue to SLATEC/LDOC format.  (FNF)
02297 !C  890503  Minor cosmetic changes.  (FNF)
02298 !C  930809  Renamed to allow single/double precision versions. (ACH)
02299 !C  010418  Reduced size of Common block /DLS001/. (ACH)
02300 !C  031105  Restored 'own' variables to Common block /DLS001/, to
02301 !C          enable interrupt/restart feature. (ACH)
02302 !C**END PROLOGUE  DSOLSY
02303 !C*End
02304       INTEGER IWM
02305       DOUBLE PRECISION WM, X, TEM
02306       DIMENSION WM(*), IWM(*), X(*), TEM(*)
02307       INTEGER IOWND, IOWNS, 
02308         ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 
02309         LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 
02310         MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
02311       DOUBLE PRECISION ROWNS, 
02312         CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
02313       COMMON /DLS001/ ROWNS(209), &
02314         CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
02315         IOWND(6), IOWNS(6), &
02316         ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
02317         LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
02318         MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
02319       INTEGER I, MEBAND, ML, MU
02320       DOUBLE PRECISION DI, HL0, PHL0, R
02321 !C
02322 !C**FIRST EXECUTABLE STATEMENT  DSOLSY
02323       IERSL = 0
02324       GO TO (100, 100, 300, 400, 400), MITER
02325  100  CALL DGESL (WM(3), N, N, IWM(21), X, 0)
02326       RETURN
02327 !C
02328  300  PHL0 = WM(2)
02329       HL0 = H*EL0
02330       WM(2) = HL0
02331       IF (HL0 .EQ. PHL0) GO TO 330
02332       R = HL0/PHL0
02333       DO 320 I = 1,N
02334         DI = 1.0D0 - R*(1.0D0 - 1.0D0/WM(I+2))
02335         IF (DABS(DI) .EQ. 0.0D0) GO TO 390
02336  320    WM(I+2) = 1.0D0/DI
02337  330  DO 340 I = 1,N
02338  340    X(I) = WM(I+2)*X(I)
02339       RETURN
02340  390  IERSL = 1
02341       RETURN
02342 !C
02343  400  ML = IWM(1)
02344       MU = IWM(2)
02345       MEBAND = 2*ML + MU + 1
02346       CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0)
02347       RETURN
02348 !C---------------------- END OF SUBROUTINE DSOLSY ----------------------
02349       END
02350 !*DECK DSRCOM
02351       SUBROUTINE DSRCOM (RSAV, ISAV, JOB)
02352 !C**BEGIN PROLOGUE  DSRCOM
02353 !C**SUBSIDIARY
02354 !C**PURPOSE  Save/restore ODEPACK COMMON blocks.
02355 !C**TYPE      DOUBLE PRECISION (SSRCOM-S, DSRCOM-D)
02356 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
02357 !C**DESCRIPTION
02358 !C
02359 !C This routine saves or restores (depending on JOB) the contents of
02360 !C the COMMON block DLS001, which is used internally
02361 !C by one or more ODEPACK solvers.
02362 !C
02363 !C RSAV = real array of length 218 or more.
02364 !C ISAV = integer array of length 37 or more.
02365 !C JOB  = flag indicating to save or restore the COMMON blocks:
02366 !C        JOB  = 1 if COMMON is to be saved (written to RSAV/ISAV)
02367 !C        JOB  = 2 if COMMON is to be restored (read from RSAV/ISAV)
02368 !C        A call with JOB = 2 presumes a prior call with JOB = 1.
02369 !C
02370 !C**SEE ALSO  DLSODE
02371 !C**ROUTINES CALLED  (NONE)
02372 !C**COMMON BLOCKS    DLS001
02373 !C**REVISION HISTORY  (YYMMDD)
02374 !C  791129  DATE WRITTEN
02375 !C  890501  Modified prologue to SLATEC/LDOC format.  (FNF)
02376 !C  890503  Minor cosmetic changes.  (FNF)
02377 !C  921116  Deleted treatment of block /EH0001/.  (ACH)
02378 !C  930801  Reduced Common block length by 2.  (ACH)
02379 !C  930809  Renamed to allow single/double precision versions. (ACH)
02380 !C  010418  Reduced Common block length by 209+12. (ACH)
02381 !C  031105  Restored 'own' variables to Common block /DLS001/, to
02382 !C          enable interrupt/restart feature. (ACH)
02383 !C  031112  Added SAVE statement for data-loaded constants.
02384 !C**END PROLOGUE  DSRCOM
02385 !C*End
02386       INTEGER ISAV, JOB
02387       INTEGER ILS
02388       INTEGER I, LENILS, LENRLS
02389       DOUBLE PRECISION RSAV,   RLS
02390       DIMENSION RSAV(*), ISAV(*)
02391       SAVE LENRLS, LENILS
02392       COMMON /DLS001/ RLS(218), ILS(37)
02393       DATA LENRLS/218/, LENILS/37/
02394 !C
02395 !C**FIRST EXECUTABLE STATEMENT  DSRCOM
02396       IF (JOB .EQ. 2) GO TO 100
02397 !C
02398       DO 10 I = 1,LENRLS
02399  10     RSAV(I) = RLS(I)
02400       DO 20 I = 1,LENILS
02401  20     ISAV(I) = ILS(I)
02402       RETURN
02403 !C
02404  100  CONTINUE
02405       DO 110 I = 1,LENRLS
02406  110     RLS(I) = RSAV(I)
02407       DO 120 I = 1,LENILS
02408  120     ILS(I) = ISAV(I)
02409       RETURN
02410 !C---------------------- END OF SUBROUTINE DSRCOM ----------------------
02411       END
02412 !*DECK DSTODE
02413       SUBROUTINE DSTODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, &
02414          WM, IWM, F, JAC, PJAC, SLVS)
02415 !C**BEGIN PROLOGUE  DSTODE
02416 !C**SUBSIDIARY
02417 !C**PURPOSE  Performs one step of an ODEPACK integration.
02418 !C**TYPE      DOUBLE PRECISION (SSTODE-S, DSTODE-D)
02419 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
02420 !C**DESCRIPTION
02421 !C
02422 !C DSTODE performs one step of the integration of an initial value
02423 !C problem for a system of ordinary differential equations.
02424 !C Note:  DSTODE is independent of the value of the iteration method
02425 !C indicator MITER, when this is .ne. 0, and hence is independent
02426 !C of the type of chord method used, or the Jacobian structure.
02427 !C Communication with DSTODE is done with the following variables:
02428 !C
02429 !C NEQ    = integer array containing problem size in NEQ(1), and
02430 !C          passed as the NEQ argument in all calls to F and JAC.
02431 !C Y      = an array of length .ge. N used as the Y argument in
02432 !C          all calls to F and JAC.
02433 !C YH     = an NYH by LMAX array containing the dependent variables
02434 !C          and their approximate scaled derivatives, where
02435 !C          LMAX = MAXORD + 1.  YH(i,j+1) contains the approximate
02436 !C          j-th derivative of y(i), scaled by h**j/factorial(j)
02437 !C          (j = 0,1,...,NQ).  on entry for the first step, the first
02438 !C          two columns of YH must be set from the initial values.
02439 !C NYH    = a constant integer .ge. N, the first dimension of YH.
02440 !C YH1    = a one-dimensional array occupying the same space as YH.
02441 !C EWT    = an array of length N containing multiplicative weights
02442 !C          for local error measurements.  Local errors in Y(i) are
02443 !C          compared to 1.0/EWT(i) in various error tests.
02444 !C SAVF   = an array of working storage, of length N.
02445 !C          Also used for input of YH(*,MAXORD+2) when JSTART = -1
02446 !C          and MAXORD .lt. the current order NQ.
02447 !C ACOR   = a work array of length N, used for the accumulated
02448 !C          corrections.  On a successful return, ACOR(i) contains
02449 !C          the estimated one-step local error in Y(i).
02450 !C WM,IWM = real and integer work arrays associated with matrix
02451 !C          operations in chord iteration (MITER .ne. 0).
02452 !C PJAC   = name of routine to evaluate and preprocess Jacobian matrix
02453 !C          and P = I - h*el0*JAC, if a chord method is being used.
02454 !C SLVS   = name of routine to solve linear system in chord iteration.
02455 !C CCMAX  = maximum relative change in h*el0 before PJAC is called.
02456 !C H      = the step size to be attempted on the next step.
02457 !C          H is altered by the error control algorithm during the
02458 !C          problem.  H can be either positive or negative, but its
02459 !C          sign must remain constant throughout the problem.
02460 !C HMIN   = the minimum absolute value of the step size h to be used.
02461 !C HMXI   = inverse of the maximum absolute value of h to be used.
02462 !C          HMXI = 0.0 is allowed and corresponds to an infinite hmax.
02463 !C          HMIN and HMXI may be changed at any time, but will not
02464 !C          take effect until the next change of h is considered.
02465 !C TN     = the independent variable. TN is updated on each step taken.
02466 !C JSTART = an integer used for input only, with the following
02467 !C          values and meanings:
02468 !C               0  perform the first step.
02469 !C           .gt.0  take a new step continuing from the last.
02470 !C              -1  take the next step with a new value of H, MAXORD,
02471 !C                    N, METH, MITER, and/or matrix parameters.
02472 !C              -2  take the next step with a new value of H,
02473 !C                    but with other inputs unchanged.
02474 !C          On return, JSTART is set to 1 to facilitate continuation.
02475 !C KFLAG  = a completion code with the following meanings:
02476 !C               0  the step was succesful.
02477 !C              -1  the requested error could not be achieved.
02478 !C              -2  corrector convergence could not be achieved.
02479 !C              -3  fatal error in PJAC or SLVS.
02480 !C          A return with KFLAG = -1 or -2 means either
02481 !C          abs(H) = HMIN or 10 consecutive failures occurred.
02482 !C          On a return with KFLAG negative, the values of TN and
02483 !C          the YH array are as of the beginning of the last
02484 !C          step, and H is the last step size attempted.
02485 !C MAXORD = the maximum order of integration method to be allowed.
02486 !C MAXCOR = the maximum number of corrector iterations allowed.
02487 !C MSBP   = maximum number of steps between PJAC calls (MITER .gt. 0).
02488 !C MXNCF  = maximum number of convergence failures allowed.
02489 !C METH/MITER = the method flags.  See description in driver.
02490 !C N      = the number of first-order differential equations.
02491 !C The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD,
02492 !C MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON.
02493 !C
02494 !C**SEE ALSO  DLSODE
02495 !C**ROUTINES CALLED  DCFODE, DVNORM
02496 !C**COMMON BLOCKS    DLS001
02497 !C**REVISION HISTORY  (YYMMDD)
02498 !C  791129  DATE WRITTEN
02499 !C  890501  Modified prologue to SLATEC/LDOC format.  (FNF)
02500 !C  890503  Minor cosmetic changes.  (FNF)
02501 !C  930809  Renamed to allow single/double precision versions. (ACH)
02502 !C  010418  Reduced size of Common block /DLS001/. (ACH)
02503 !C  031105  Restored 'own' variables to Common block /DLS001/, to
02504 !C          enable interrupt/restart feature. (ACH)
02505 !C**END PROLOGUE  DSTODE
02506 !C*End
02507       EXTERNAL F, JAC, PJAC, SLVS
02508       INTEGER NEQ, NYH, IWM
02509       DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM
02510       DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*), &
02511          ACOR(*), WM(*), IWM(*)
02512       INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, 
02513         ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 
02514         LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 
02515         MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
02516       INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
02517       DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, 
02518          CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
02519       DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP, 
02520          R, RH, RHDN, RHSM, RHUP, TOLD, DVNORM
02521       COMMON /DLS001/ CONIT, CRATE, EL(13), ELCO(13,12), &
02522         HOLD, RMAX, TESCO(3,12), &
02523         CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
02524         IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, &
02525         ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
02526         LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
02527         MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
02528 !C
02529 !C**FIRST EXECUTABLE STATEMENT  DSTODE
02530       KFLAG = 0
02531       TOLD = TN
02532       NCF = 0
02533       IERPJ = 0
02534       IERSL = 0
02535       JCUR = 0
02536       ICF = 0
02537       DELP = 0.0D0
02538       IF (JSTART .GT. 0) GO TO 200
02539       IF (JSTART .EQ. -1) GO TO 100
02540       IF (JSTART .EQ. -2) GO TO 160
02541 !C----------------------------------------------------------------------
02542 !COn the first call, the order is set to 1, and other variables are
02543 !Cinitialized.  RMAX is the maximum ratio by which H can be increased
02544 !Cin a single step.  It is initially 1.E4 to compensate for the small
02545 !Cinitial H, but then is normally equal to 10.  If a failure
02546 !Coccurs (in corrector convergence or error test), RMAX is set to 2
02547 !Cfor the next increase.
02548 !C----------------------------------------------------------------------
02549       LMAX = MAXORD + 1
02550       NQ = 1
02551       L = 2
02552       IALTH = 2
02553       RMAX = 10000.0D0
02554       RC = 0.0D0
02555       EL0 = 1.0D0
02556       CRATE = 0.7D0
02557       HOLD = H
02558       MEO = METH
02559       NSLP = 0
02560       IPUP = MITER
02561       IRET = 3
02562       GO TO 140
02563 !C----------------------------------------------------------------------
02564 !CThe following block handles preliminaries needed when JSTART = -1.
02565 !CIPUP is set to MITER to force a matrix update.
02566 !CIf an order increase is about to be considered (IALTH = 1),
02567 !CIALTH is reset to 2 to postpone consideration one more step.
02568 !CIf the caller has changed METH, DCFODE is called to reset
02569 !Cthe coefficients of the method.
02570 !CIf the caller has changed MAXORD to a value less than the current
02571 !Corder NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
02572 !CIf H is to be changed, YH must be rescaled.
02573 !CIf H or METH is being changed, IALTH is reset to L = NQ + 1
02574 !Cto prevent further changes in H for that many steps.
02575 !C----------------------------------------------------------------------
02576  100  IPUP = MITER
02577       LMAX = MAXORD + 1
02578       IF (IALTH .EQ. 1) IALTH = 2
02579       IF (METH .EQ. MEO) GO TO 110
02580       CALL DCFODE (METH, ELCO, TESCO)
02581       MEO = METH
02582       IF (NQ .GT. MAXORD) GO TO 120
02583       IALTH = L
02584       IRET = 1
02585       GO TO 150
02586  110  IF (NQ .LE. MAXORD) GO TO 160
02587  120  NQ = MAXORD
02588       L = LMAX
02589       DO 125 I = 1,L
02590  125    EL(I) = ELCO(I,NQ)
02591       NQNYH = NQ*NYH
02592       RC = RC*EL(1)/EL0
02593       EL0 = EL(1)
02594       CONIT = 0.5D0/(NQ+2)
02595       DDN = DVNORM (N, SAVF, EWT)/TESCO(1,L)
02596       EXDN = 1.0D0/L
02597       RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
02598       RH = DMIN1(RHDN,1.0D0)
02599       IREDO = 3
02600       IF (H .EQ. HOLD) GO TO 170
02601       RH = DMIN1(RH,DABS(H/HOLD))
02602       H = HOLD
02603       GO TO 175
02604 !C----------------------------------------------------------------------
02605 !CDCFODE is called to get all the integration coefficients for the
02606 !Ccurrent METH.  Then the EL vector and related constants are reset
02607 !Cwhenever the order NQ is changed, or at the start of the problem.
02608 !C----------------------------------------------------------------------
02609  140  CALL DCFODE (METH, ELCO, TESCO)
02610  150  DO 155 I = 1,L
02611  155    EL(I) = ELCO(I,NQ)
02612       NQNYH = NQ*NYH
02613       RC = RC*EL(1)/EL0
02614       EL0 = EL(1)
02615       CONIT = 0.5D0/(NQ+2)
02616       GO TO (160, 170, 200), IRET
02617 !C----------------------------------------------------------------------
02618 !CIf H is being changed, the H ratio RH is checked against
02619 !CRMAX, HMIN, and HMXI, and the YH array rescaled.  IALTH is set to
02620 !CL = NQ + 1 to prevent a change of H for that many steps, unless
02621 !Cforced by a convergence or error test failure.
02622 !C----------------------------------------------------------------------
02623  160  IF (H .EQ. HOLD) GO TO 200
02624       RH = H/HOLD
02625       H = HOLD
02626       IREDO = 3
02627       GO TO 175
02628  170  RH = DMAX1(RH,HMIN/DABS(H))
02629  175  RH = DMIN1(RH,RMAX)
02630       RH = RH/DMAX1(1.0D0,DABS(H)*HMXI*RH)
02631       R = 1.0D0
02632       DO 180 J = 2,L
02633         R = R*RH
02634         DO 180 I = 1,N
02635  180      YH(I,J) = YH(I,J)*R
02636       H = H*RH
02637       RC = RC*RH
02638       IALTH = L
02639       IF (IREDO .EQ. 0) GO TO 690
02640 !C----------------------------------------------------------------------
02641 !CThis section computes the predicted values by effectively
02642 !Cmultiplying the YH array by the Pascal Triangle matrix.
02643 !CRC is the ratio of new to old values of the coefficient  H*EL(1).
02644 !CWhen RC differs from 1 by more than CCMAX, IPUP is set to MITER
02645 !Cto force PJAC to be called, if a Jacobian is involved.
02646 !CIn any case, PJAC is called at least every MSBP steps.
02647 !C----------------------------------------------------------------------
02648  200  IF (DABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER
02649       IF (NST .GE. NSLP+MSBP) IPUP = MITER
02650       TN = TN + H
02651       I1 = NQNYH + 1
02652       DO 215 JB = 1,NQ
02653         I1 = I1 - NYH
02654 !Cir$ ivdep
02655         DO 210 I = I1,NQNYH
02656  210      YH1(I) = YH1(I) + YH1(I+NYH)
02657  215    CONTINUE
02658 !C----------------------------------------------------------------------
02659 !CUp to MAXCOR corrector iterations are taken.  A convergence test is
02660 !Cmade on the R.M.S. norm of each correction, weighted by the error
02661 !Cweight vector EWT.  The sum of the corrections is accumulated in the
02662 !Cvector ACOR(i).  The YH array is not altered in the corrector loop.
02663 !C----------------------------------------------------------------------
02664  220  M = 0
02665       DO 230 I = 1,N
02666  230    Y(I) = YH(I,1)
02667       CALL F (NEQ, TN, Y, SAVF)
02668       NFE = NFE + 1
02669       IF (IPUP .LE. 0) GO TO 250
02670 !C----------------------------------------------------------------------
02671 !CIf indicated, the matrix P = I - h*el(1)*J is reevaluated and
02672 !Cpreprocessed before starting the corrector iteration.  IPUP is set
02673 !Cto 0 as an indicator that this has been done.
02674 !C----------------------------------------------------------------------
02675       CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC)
02676       IPUP = 0
02677       RC = 1.0D0
02678       NSLP = NST
02679       CRATE = 0.7D0
02680       IF (IERPJ .NE. 0) GO TO 430
02681  250  DO 260 I = 1,N
02682  260    ACOR(I) = 0.0D0
02683  270  IF (MITER .NE. 0) GO TO 350
02684 !C----------------------------------------------------------------------
02685 !CIn the case of functional iteration, update Y directly from
02686 !Cthe result of the last function evaluation.
02687 !C----------------------------------------------------------------------
02688       DO 290 I = 1,N
02689         SAVF(I) = H*SAVF(I) - YH(I,2)
02690  290    Y(I) = SAVF(I) - ACOR(I)
02691       DEL = DVNORM (N, Y, EWT)
02692       DO 300 I = 1,N
02693         Y(I) = YH(I,1) + EL(1)*SAVF(I)
02694  300    ACOR(I) = SAVF(I)
02695       GO TO 400
02696 !C----------------------------------------------------------------------
02697 !CIn the case of the chord method, compute the corrector error,
02698 !Cand solve the linear system with that as right-hand side and
02699 !CP as coefficient matrix.
02700 !C----------------------------------------------------------------------
02701  350  DO 360 I = 1,N
02702  360    Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
02703       CALL SLVS (WM, IWM, Y, SAVF)
02704       IF (IERSL .LT. 0) GO TO 430
02705       IF (IERSL .GT. 0) GO TO 410
02706       DEL = DVNORM (N, Y, EWT)
02707       DO 380 I = 1,N
02708         ACOR(I) = ACOR(I) + Y(I)
02709  380    Y(I) = YH(I,1) + EL(1)*ACOR(I)
02710 !C----------------------------------------------------------------------
02711 !CTest for convergence.  If M.gt.0, an estimate of the convergence
02712 !Crate constant is stored in CRATE, and this is used in the test.
02713 !C----------------------------------------------------------------------
02714  400  IF (M .NE. 0) CRATE = DMAX1(0.2D0*CRATE,DEL/DELP)
02715       DCON = DEL*DMIN1(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT)
02716       IF (DCON .LE. 1.0D0) GO TO 450
02717       M = M + 1
02718       IF (M .EQ. MAXCOR) GO TO 410
02719       IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
02720       DELP = DEL
02721       CALL F (NEQ, TN, Y, SAVF)
02722       NFE = NFE + 1
02723       GO TO 270
02724 !C----------------------------------------------------------------------
02725 !CThe corrector iteration failed to converge.
02726 !CIf MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
02727 !Cthe next try.  Otherwise the YH array is retracted to its values
02728 !Cbefore prediction, and H is reduced, if possible.  If H cannot be
02729 !Creduced or MXNCF failures have occurred, exit with KFLAG = -2.
02730 !C----------------------------------------------------------------------
02731  410  IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
02732       ICF = 1
02733       IPUP = MITER
02734       GO TO 220
02735  430  ICF = 2
02736       NCF = NCF + 1
02737       RMAX = 2.0D0
02738       TN = TOLD
02739       I1 = NQNYH + 1
02740       DO 445 JB = 1,NQ
02741         I1 = I1 - NYH
02742 !Cir$ ivdep
02743         DO 440 I = I1,NQNYH
02744  440      YH1(I) = YH1(I) - YH1(I+NYH)
02745  445    CONTINUE
02746       IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
02747       IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 670
02748       IF (NCF .EQ. MXNCF) GO TO 670
02749       RH = 0.25D0
02750       IPUP = MITER
02751       IREDO = 1
02752       GO TO 170
02753 !C----------------------------------------------------------------------
02754 !CThe corrector has converged.  JCUR is set to 0
02755 !Cto signal that the Jacobian involved may need updating later.
02756 !CThe local error test is made and control passes to statement 500
02757 !Cif it fails.
02758 !C----------------------------------------------------------------------
02759  450  JCUR = 0
02760       IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
02761       IF (M .GT. 0) DSM = DVNORM (N, ACOR, EWT)/TESCO(2,NQ)
02762       IF (DSM .GT. 1.0D0) GO TO 500
02763 !C----------------------------------------------------------------------
02764 !CAfter a successful step, update the YH array.
02765 !CConsider changing H if IALTH = 1.  Otherwise decrease IALTH by 1.
02766 !CIf IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
02767 !Cuse in a possible order increase on the next step.
02768 !CIf a change in H is considered, an increase or decrease in order
02769 !Cby one is considered also.  A change in H is made only if it is by a
02770 !Cfactor of at least 1.1.  If not, IALTH is set to 3 to prevent
02771 !Ctesting for that many steps.
02772 !C----------------------------------------------------------------------
02773       KFLAG = 0
02774       IREDO = 0
02775       NST = NST + 1
02776       HU = H
02777       NQU = NQ
02778       DO 470 J = 1,L
02779         DO 470 I = 1,N
02780  470      YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
02781       IALTH = IALTH - 1
02782       IF (IALTH .EQ. 0) GO TO 520
02783       IF (IALTH .GT. 1) GO TO 700
02784       IF (L .EQ. LMAX) GO TO 700
02785       DO 490 I = 1,N
02786  490    YH(I,LMAX) = ACOR(I)
02787       GO TO 700
02788 !C----------------------------------------------------------------------
02789 !CThe error test failed.  KFLAG keeps track of multiple failures.
02790 !CRestore TN and the YH array to their previous values, and prepare
02791 !Cto try the step again.  Compute the optimum step size for this or
02792 !Cone lower order.  After 2 or more failures, H is forced to decrease
02793 !Cby a factor of 0.2 or less.
02794 !C----------------------------------------------------------------------
02795  500  KFLAG = KFLAG - 1
02796       TN = TOLD
02797       I1 = NQNYH + 1
02798       DO 515 JB = 1,NQ
02799         I1 = I1 - NYH
02800 !Cir$ ivdep
02801         DO 510 I = I1,NQNYH
02802  510      YH1(I) = YH1(I) - YH1(I+NYH)
02803  515    CONTINUE
02804       RMAX = 2.0D0
02805       IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 660
02806       IF (KFLAG .LE. -3) GO TO 640
02807       IREDO = 2
02808       RHUP = 0.0D0
02809       GO TO 540
02810 !C----------------------------------------------------------------------
02811 !CRegardless of the success or failure of the step, factors
02812 !CRHDN, RHSM, and RHUP are computed, by which H could be multiplied
02813 !Cat order NQ - 1, order NQ, or order NQ + 1, respectively.
02814 !CIn the case of failure, RHUP = 0.0 to avoid an order increase.
02815 !CThe largest of these is determined and the new order chosen
02816 !Caccordingly.  If the order is to be increased, we compute one
02817 !Cadditional scaled derivative.
02818 !C----------------------------------------------------------------------
02819  520  RHUP = 0.0D0
02820       IF (L .EQ. LMAX) GO TO 540
02821       DO 530 I = 1,N
02822  530    SAVF(I) = ACOR(I) - YH(I,LMAX)
02823       DUP = DVNORM (N, SAVF, EWT)/TESCO(3,NQ)
02824       EXUP = 1.0D0/(L+1)
02825       RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0)
02826  540  EXSM = 1.0D0/L
02827       RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
02828       RHDN = 0.0D0
02829       IF (NQ .EQ. 1) GO TO 560
02830       DDN = DVNORM (N, YH(1,L), EWT)/TESCO(1,NQ)
02831       EXDN = 1.0D0/NQ
02832       RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
02833  560  IF (RHSM .GE. RHUP) GO TO 570
02834       IF (RHUP .GT. RHDN) GO TO 590
02835       GO TO 580
02836  570  IF (RHSM .LT. RHDN) GO TO 580
02837       NEWQ = NQ
02838       RH = RHSM
02839       GO TO 620
02840  580  NEWQ = NQ - 1
02841       RH = RHDN
02842       IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0
02843       GO TO 620
02844  590  NEWQ = L
02845       RH = RHUP
02846       IF (RH .LT. 1.1D0) GO TO 610
02847       R = EL(L)/L
02848       DO 600 I = 1,N
02849  600    YH(I,NEWQ+1) = ACOR(I)*R
02850       GO TO 630
02851  610  IALTH = 3
02852       GO TO 700
02853  620  IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610
02854       IF (KFLAG .LE. -2) RH = DMIN1(RH,0.2D0)
02855 !C----------------------------------------------------------------------
02856 !CIf there is a change of order, reset NQ, l, and the coefficients.
02857 !CIn any case H is reset according to RH and the YH array is rescaled.
02858 !CThen exit from 690 if the step was OK, or redo the step otherwise.
02859 !C----------------------------------------------------------------------
02860       IF (NEWQ .EQ. NQ) GO TO 170
02861  630  NQ = NEWQ
02862       L = NQ + 1
02863       IRET = 2
02864       GO TO 150
02865 !C----------------------------------------------------------------------
02866 !CControl reaches this section if 3 or more failures have occured.
02867 !CIf 10 failures have occurred, exit with KFLAG = -1.
02868 !CIt is assumed that the derivatives that have accumulated in the
02869 !CYH array have errors of the wrong order.  Hence the first
02870 !Cderivative is recomputed, and the order is set to 1.  Then
02871 !CH is reduced by a factor of 10, and the step is retried,
02872 !Cuntil it succeeds or H reaches HMIN.
02873 !C----------------------------------------------------------------------
02874  640  IF (KFLAG .EQ. -10) GO TO 660
02875       RH = 0.1D0
02876       RH = DMAX1(HMIN/DABS(H),RH)
02877       H = H*RH
02878       DO 645 I = 1,N
02879  645    Y(I) = YH(I,1)
02880       CALL F (NEQ, TN, Y, SAVF)
02881       NFE = NFE + 1
02882       DO 650 I = 1,N
02883  650    YH(I,2) = H*SAVF(I)
02884       IPUP = MITER
02885       IALTH = 5
02886       IF (NQ .EQ. 1) GO TO 200
02887       NQ = 1
02888       L = 2
02889       IRET = 3
02890       GO TO 150
02891 !C----------------------------------------------------------------------
02892 !CAll returns are made through this section.  H is saved in HOLD
02893 !Cto allow the caller to change H on the next step.
02894 !C----------------------------------------------------------------------
02895  660  KFLAG = -1
02896       GO TO 720
02897  670  KFLAG = -2
02898       GO TO 720
02899  680  KFLAG = -3
02900       GO TO 720
02901  690  RMAX = 10.0D0
02902  700  R = 1.0D0/TESCO(2,NQU)
02903       DO 710 I = 1,N
02904  710    ACOR(I) = ACOR(I)*R
02905  720  HOLD = H
02906       JSTART = 1
02907       RETURN
02908 !C---------------------- END OF SUBROUTINE DSTODE ----------------------
02909       END
02910 !*DECK DEWSET
02911       SUBROUTINE DEWSET (N, ITOL, RTOL, ATOL, YCUR, EWT)
02912 !C**BEGIN PROLOGUE  DEWSET
02913 !C**SUBSIDIARY
02914 !C**PURPOSE  Set error weight vector.
02915 !C**TYPE      DOUBLE PRECISION (SEWSET-S, DEWSET-D)
02916 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
02917 !C**DESCRIPTION
02918 !C
02919 !C This subroutine sets the error weight vector EWT according to
02920 !C     EWT(i) = RTOL(i)*ABS(YCUR(i)) + ATOL(i),  i = 1,...,N,
02921 !C with the subscript on RTOL and/or ATOL possibly replaced by 1 above,
02922 !C depending on the value of ITOL.
02923 !C
02924 !C**SEE ALSO  DLSODE
02925 !C**ROUTINES CALLED  (NONE)
02926 !C**REVISION HISTORY  (YYMMDD)
02927 !C  791129  DATE WRITTEN
02928 !C  890501  Modified prologue to SLATEC/LDOC format.  (FNF)
02929 !C  890503  Minor cosmetic changes.  (FNF)
02930 !C  930809  Renamed to allow single/double precision versions. (ACH)
02931 !C**END PROLOGUE  DEWSET
02932 !C*End
02933       INTEGER N, ITOL
02934       INTEGER I
02935       DOUBLE PRECISION RTOL, ATOL, YCUR, EWT
02936       DIMENSION RTOL(*), ATOL(*), YCUR(N), EWT(N)
02937 !C
02938 !C**FIRST EXECUTABLE STATEMENT  DEWSET
02939       GO TO (10, 20, 30, 40), ITOL
02940  10   CONTINUE
02941       DO 15 I = 1,N
02942  15     EWT(I) = RTOL(1)*DABS(YCUR(I)) + ATOL(1)
02943       RETURN
02944  20   CONTINUE
02945       DO 25 I = 1,N
02946  25     EWT(I) = RTOL(1)*DABS(YCUR(I)) + ATOL(I)
02947       RETURN
02948  30   CONTINUE
02949       DO 35 I = 1,N
02950  35     EWT(I) = RTOL(I)*DABS(YCUR(I)) + ATOL(1)
02951       RETURN
02952  40   CONTINUE
02953       DO 45 I = 1,N
02954  45     EWT(I) = RTOL(I)*DABS(YCUR(I)) + ATOL(I)
02955       RETURN
02956 !C---------------------- END OF SUBROUTINE DEWSET ----------------------
02957       END
02958 !*DECK DVNORM
02959       DOUBLE PRECISION FUNCTION DVNORM (N, V, W)
02960 !C**BEGIN PROLOGUE  DVNORM
02961 !C**SUBSIDIARY
02962 !C**PURPOSE  Weighted root-mean-square vector norm.
02963 !C**TYPE      DOUBLE PRECISION (SVNORM-S, DVNORM-D)
02964 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
02965 !C**DESCRIPTION
02966 !C
02967 !C This function routine computes the weighted root-mean-square norm
02968 !C of the vector of length N contained in the array V, with weights
02969 !C contained in the array W of length N:
02970 !C   DVNORM = SQRT( (1/N) * SUM( V(i)*W(i) )**2 )
02971 !C
02972 !C**SEE ALSO  DLSODE
02973 !C**ROUTINES CALLED  (NONE)
02974 !C**REVISION HISTORY  (YYMMDD)
02975 !C  791129  DATE WRITTEN
02976 !C  890501  Modified prologue to SLATEC/LDOC format.  (FNF)
02977 !C  890503  Minor cosmetic changes.  (FNF)
02978 !C  930809  Renamed to allow single/double precision versions. (ACH)
02979 !C**END PROLOGUE  DVNORM
02980 !C*End
02981       INTEGER N,   I
02982       DOUBLE PRECISION V, W,   SUM
02983       DIMENSION V(N), W(N)
02984 !C
02985 !C**FIRST EXECUTABLE STATEMENT  DVNORM
02986       SUM = 0.0D0
02987       DO 10 I = 1,N
02988  10     SUM = SUM + (V(I)*W(I))**2
02989       DVNORM = DSQRT(SUM/dble(N))
02990       RETURN
02991 !C---------------------- END OF FUNCTION DVNORM ------------------------
02992       END
02993 !*DECK DGEFA
02994       SUBROUTINE DGEFA (A, LDA, N, IPVT, INFO)
02995 !C**BEGIN PROLOGUE  DGEFA
02996 !C**PURPOSE  Factor a matrix using Gaussian elimination.
02997 !C**CATEGORY  D2A1
02998 !C**TYPE      DOUBLE PRECISION (SGEFA-S, DGEFA-D, CGEFA-C)
02999 !C**KEYWORDS  GENERAL MATRIX, LINEAR ALGEBRA, LINPACK,
03000 !C            MATRIX FACTORIZATION
03001 !C**AUTHOR  Moler, C. B., (U. of New Mexico)
03002 !C**DESCRIPTION
03003 !C
03004 !C    DGEFA factors a double precision matrix by Gaussian elimination.
03005 !C
03006 !C    DGEFA is usually called by DGECO, but it can be called
03007 !C    directly with a saving in time if  RCOND  is not needed.
03008 !C    (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) .
03009 !C
03010 !C    On Entry
03011 !C
03012 !C       A       DOUBLE PRECISION(LDA, N)
03013 !C               the matrix to be factored.
03014 !C
03015 !C       LDA     INTEGER
03016 !C               the leading dimension of the array  A .
03017 !C
03018 !C       N       INTEGER
03019 !C               the order of the matrix  A .
03020 !C
03021 !C    On Return
03022 !C
03023 !C       A       an upper triangular matrix and the multipliers
03024 !C               which were used to obtain it.
03025 !C               The factorization can be written  A = L*U  where
03026 !C               L  is a product of permutation and unit lower
03027 !C               triangular matrices and  U  is upper triangular.
03028 !C
03029 !C       IPVT    INTEGER(N)
03030 !C               an integer vector of pivot indices.
03031 !C
03032 !C       INFO    INTEGER
03033 !C               = 0  normal value.
03034 !C               = K  if  U(K,K) .EQ. 0.0 .  This is not an error
03035 !C                    condition for this subroutine, but it does
03036 !C                    indicate that DGESL or DGEDI will divide by zero
03037 !C                    if called.  Use  RCOND  in DGECO for a reliable
03038 !C                    indication of singularity.
03039 !C
03040 !C**REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
03041 !C                Stewart, LINPACK Users' Guide, SIAM, 1979.
03042 !C**ROUTINES CALLED  DAXPY, DSCAL, IDAMAX
03043 !C**REVISION HISTORY  (YYMMDD)
03044 !C  780814  DATE WRITTEN
03045 !C  890831  Modified array declarations.  (WRB)
03046 !C  890831  REVISION DATE from Version 3.2
03047 !C  891214  Prologue converted to Version 4.0 format.  (BAB)
03048 !C  900326  Removed duplicate information from DESCRIPTION section.
03049 !C          (WRB)
03050 !C  920501  Reformatted the REFERENCES section.  (WRB)
03051 !C**END PROLOGUE  DGEFA
03052       INTEGER LDA,N,IPVT(*),INFO
03053       DOUBLE PRECISION A(LDA,*)
03054 !C
03055       DOUBLE PRECISION T
03056       INTEGER IDAMAX,J,K,KP1,L,NM1
03057 !C
03058 !C    GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
03059 !C
03060 !C**FIRST EXECUTABLE STATEMENT  DGEFA
03061       INFO = 0
03062       NM1 = N - 1
03063       IF (NM1 .LT. 1) GO TO 70
03064       DO 60 K = 1, NM1
03065          KP1 = K + 1
03066 !C
03067 !C       FIND L = PIVOT INDEX
03068 !C
03069          L = IDAMAX(N-K+1,A(K,K),1) + K - 1
03070          IPVT(K) = L
03071 !C
03072 !C       ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
03073 !C
03074          IF (A(L,K) .EQ. 0.0D0) GO TO 40
03075 !C
03076 !C          INTERCHANGE IF NECESSARY
03077 !C
03078             IF (L .EQ. K) GO TO 10
03079                T = A(L,K)
03080                A(L,K) = A(K,K)
03081                A(K,K) = T
03082    10       CONTINUE
03083 !C
03084 !C          COMPUTE MULTIPLIERS
03085 !C
03086             T = -1.0D0/A(K,K)
03087             CALL DSCAL(N-K,T,A(K+1,K),1)
03088 !C
03089 !C          ROW ELIMINATION WITH COLUMN INDEXING
03090 !C
03091             DO 30 J = KP1, N
03092                T = A(L,J)
03093                IF (L .EQ. K) GO TO 20
03094                   A(L,J) = A(K,J)
03095                   A(K,J) = T
03096    20          CONTINUE
03097                CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
03098    30       CONTINUE
03099          GO TO 50
03100    40    CONTINUE
03101             INFO = K
03102    50    CONTINUE
03103    60 CONTINUE
03104    70 CONTINUE
03105       IPVT(N) = N
03106       IF (A(N,N) .EQ. 0.0D0) INFO = N
03107       RETURN
03108       END
03109 !*DECK DGESL
03110       SUBROUTINE DGESL (A, LDA, N, IPVT, B, JOB)
03111 !C**BEGIN PROLOGUE  DGESL
03112 !C**PURPOSE  Solve the real system A*X=B or TRANS(A)*X=B using the
03113 !C           factors computed by DGECO or DGEFA.
03114 !C**CATEGORY  D2A1
03115 !C**TYPE      DOUBLE PRECISION (SGESL-S, DGESL-D, CGESL-C)
03116 !C**KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
03117 !C**AUTHOR  Moler, C. B., (U. of New Mexico)
03118 !C**DESCRIPTION
03119 !C
03120 !C    DGESL solves the double precision system
03121 !C    A * X = B  or  TRANS(A) * X = B
03122 !C    using the factors computed by DGECO or DGEFA.
03123 !C
03124 !C    On Entry
03125 !C
03126 !C       A       DOUBLE PRECISION(LDA, N)
03127 !C               the output from DGECO or DGEFA.
03128 !C
03129 !C       LDA     INTEGER
03130 !C               the leading dimension of the array  A .
03131 !C
03132 !C       N       INTEGER
03133 !C               the order of the matrix  A .
03134 !C
03135 !C       IPVT    INTEGER(N)
03136 !C               the pivot vector from DGECO or DGEFA.
03137 !C
03138 !C       B       DOUBLE PRECISION(N)
03139 !C               the right hand side vector.
03140 !C
03141 !C       JOB     INTEGER
03142 !C               = 0         to solve  A*X = B ,
03143 !C               = nonzero   to solve  TRANS(A)*X = B  where
03144 !C                           TRANS(A)  is the transpose.
03145 !C
03146 !C    On Return
03147 !C
03148 !C       B       the solution vector  X .
03149 !C
03150 !C    Error Condition
03151 !C
03152 !C       A division by zero will occur if the input factor contains a
03153 !C       zero on the diagonal.  Technically this indicates singularity
03154 !C       but it is often caused by improper arguments or improper
03155 !C       setting of LDA .  It will not occur if the subroutines are
03156 !C       called correctly and if DGECO has set RCOND .GT. 0.0
03157 !C       or DGEFA has set INFO .EQ. 0 .
03158 !C
03159 !C    To compute  INVERSE(A) * C  where  C  is a matrix
03160 !C    with  P  columns
03161 !C          CALL DGECO(A,LDA,N,IPVT,RCOND,Z)
03162 !C          IF (RCOND is too small) GO TO ...
03163 !C          DO 10 J = 1, P
03164 !C             CALL DGESL(A,LDA,N,IPVT,C(1,J),0)
03165 !C       10 CONTINUE
03166 !C
03167 !C**REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
03168 !C                Stewart, LINPACK Users' Guide, SIAM, 1979.
03169 !C**ROUTINES CALLED  DAXPY, DDOT
03170 !C**REVISION HISTORY  (YYMMDD)
03171 !C  780814  DATE WRITTEN
03172 !C  890831  Modified array declarations.  (WRB)
03173 !C  890831  REVISION DATE from Version 3.2
03174 !C  891214  Prologue converted to Version 4.0 format.  (BAB)
03175 !C  900326  Removed duplicate information from DESCRIPTION section.
03176 !C          (WRB)
03177 !C  920501  Reformatted the REFERENCES section.  (WRB)
03178 !C**END PROLOGUE  DGESL
03179       INTEGER LDA,N,IPVT(*),JOB
03180       DOUBLE PRECISION A(LDA,*),B(*)
03181 !C
03182       DOUBLE PRECISION DDOT,T
03183       INTEGER K,KB,L,NM1
03184 !C**FIRST EXECUTABLE STATEMENT  DGESL
03185       NM1 = N - 1
03186       IF (JOB .NE. 0) GO TO 50
03187 !C
03188 !C       JOB = 0 , SOLVE  A * X = B
03189 !C       FIRST SOLVE  L*Y = B
03190 !C
03191          IF (NM1 .LT. 1) GO TO 30
03192          DO 20 K = 1, NM1
03193             L = IPVT(K)
03194             T = B(L)
03195             IF (L .EQ. K) GO TO 10
03196                B(L) = B(K)
03197                B(K) = T
03198    10       CONTINUE
03199             CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
03200    20    CONTINUE
03201    30    CONTINUE
03202 !C
03203 !C       NOW SOLVE  U*X = Y
03204 !C
03205          DO 40 KB = 1, N
03206             K = N + 1 - KB
03207             B(K) = B(K)/A(K,K)
03208             T = -B(K)
03209             CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
03210    40    CONTINUE
03211       GO TO 100
03212    50 CONTINUE
03213 !C
03214 !C       JOB = NONZERO, SOLVE  TRANS(A) * X = B
03215 !C       FIRST SOLVE  TRANS(U)*Y = B
03216 !C
03217          DO 60 K = 1, N
03218             T = DDOT(K-1,A(1,K),1,B(1),1)
03219             B(K) = (B(K) - T)/A(K,K)
03220    60    CONTINUE
03221 !C
03222 !C       NOW SOLVE TRANS(L)*X = Y
03223 !C
03224          IF (NM1 .LT. 1) GO TO 90
03225          DO 80 KB = 1, NM1
03226             K = N - KB
03227             B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1)
03228             L = IPVT(K)
03229             IF (L .EQ. K) GO TO 70
03230                T = B(L)
03231                B(L) = B(K)
03232                B(K) = T
03233    70       CONTINUE
03234    80    CONTINUE
03235    90    CONTINUE
03236   100 CONTINUE
03237       RETURN
03238       END
03239 !*DECK DGBFA
03240       SUBROUTINE DGBFA (ABD, LDA, N, ML, MU, IPVT, INFO)
03241 !C**BEGIN PROLOGUE  DGBFA
03242 !C**PURPOSE  Factor a band matrix using Gaussian elimination.
03243 !C**CATEGORY  D2A2
03244 !C**TYPE      DOUBLE PRECISION (SGBFA-S, DGBFA-D, CGBFA-C)
03245 !C**KEYWORDS  BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION
03246 !C**AUTHOR  Moler, C. B., (U. of New Mexico)
03247 !C**DESCRIPTION
03248 !C
03249 !C    DGBFA factors a double precision band matrix by elimination.
03250 !C
03251 !C    DGBFA is usually called by DGBCO, but it can be called
03252 !C    directly with a saving in time if  RCOND  is not needed.
03253 !C
03254 !C    On Entry
03255 !C
03256 !C       ABD     DOUBLE PRECISION(LDA, N)
03257 !C               contains the matrix in band storage.  The columns
03258 !C               of the matrix are stored in the columns of  ABD  and
03259 !C               the diagonals of the matrix are stored in rows
03260 !C               ML+1 through 2*ML+MU+1 of  ABD .
03261 !C               See the comments below for details.
03262 !C
03263 !C       LDA     INTEGER
03264 !C               the leading dimension of the array  ABD .
03265 !C               LDA must be .GE. 2*ML + MU + 1 .
03266 !C
03267 !C       N       INTEGER
03268 !C               the order of the original matrix.
03269 !C
03270 !C       ML      INTEGER
03271 !C               number of diagonals below the main diagonal.
03272 !C               0 .LE. ML .LT.  N .
03273 !C
03274 !C       MU      INTEGER
03275 !C               number of diagonals above the main diagonal.
03276 !C               0 .LE. MU .LT.  N .
03277 !C               More efficient if  ML .LE. MU .
03278 !C    On Return
03279 !C
03280 !C       ABD     an upper triangular matrix in band storage and
03281 !C               the multipliers which were used to obtain it.
03282 !C               The factorization can be written  A = L*U  where
03283 !C               L  is a product of permutation and unit lower
03284 !C               triangular matrices and  U  is upper triangular.
03285 !C
03286 !C       IPVT    INTEGER(N)
03287 !C               an integer vector of pivot indices.
03288 !C
03289 !C       INFO    INTEGER
03290 !C               = 0  normal value.
03291 !C               = K  if  U(K,K) .EQ. 0.0 .  This is not an error
03292 !C                    condition for this subroutine, but it does
03293 !C                    indicate that DGBSL will divide by zero if
03294 !C                    called.  Use  RCOND  in DGBCO for a reliable
03295 !C                    indication of singularity.
03296 !C
03297 !C    Band Storage
03298 !C
03299 !C          If  A  is a band matrix, the following program segment
03300 !C          will set up the input.
03301 !C
03302 !C                  ML = (band width below the diagonal)
03303 !C                  MU = (band width above the diagonal)
03304 !C                  M = ML + MU + 1
03305 !C                  DO 20 J = 1, N
03306 !C                     I1 = MAX(1, J-MU)
03307 !C                     I2 = MIN(N, J+ML)
03308 !C                     DO 10 I = I1, I2
03309 !C                        K = I - J + M
03310 !C                        ABD(K,J) = A(I,J)
03311 !C               10    CONTINUE
03312 !C               20 CONTINUE
03313 !C
03314 !C          This uses rows  ML+1  through  2*ML+MU+1  of  ABD .
03315 !C          In addition, the first  ML  rows in  ABD  are used for
03316 !C          elements generated during the triangularization.
03317 !C          The total number of rows needed in  ABD  is  2*ML+MU+1 .
03318 !C          The  ML+MU by ML+MU  upper left triangle and the
03319 !C          ML by ML  lower right triangle are not referenced.
03320 !C
03321 !C**REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
03322 !C                Stewart, LINPACK Users' Guide, SIAM, 1979.
03323 !C**ROUTINES CALLED  DAXPY, DSCAL, IDAMAX
03324 !C**REVISION HISTORY  (YYMMDD)
03325 !C  780814  DATE WRITTEN
03326 !C  890531  Changed all specific intrinsics to generic.  (WRB)
03327 !C  890831  Modified array declarations.  (WRB)
03328 !C  890831  REVISION DATE from Version 3.2
03329 !C  891214  Prologue converted to Version 4.0 format.  (BAB)
03330 !C  900326  Removed duplicate information from DESCRIPTION section.
03331 !C          (WRB)
03332 !C  920501  Reformatted the REFERENCES section.  (WRB)
03333 !C**END PROLOGUE  DGBFA
03334       INTEGER LDA,N,ML,MU,IPVT(*),INFO
03335       DOUBLE PRECISION ABD(LDA,*)
03336 !C
03337       DOUBLE PRECISION T
03338       INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1
03339 !C
03340 !C**FIRST EXECUTABLE STATEMENT  DGBFA
03341       M = ML + MU + 1
03342       INFO = 0
03343 !C
03344 !C    ZERO INITIAL FILL-IN COLUMNS
03345 !C
03346       J0 = MU + 2
03347       J1 = MIN(N,M) - 1
03348       IF (J1 .LT. J0) GO TO 30
03349       DO 20 JZ = J0, J1
03350          I0 = M + 1 - JZ
03351          DO 10 I = I0, ML
03352             ABD(I,JZ) = 0.0D0
03353    10    CONTINUE
03354    20 CONTINUE
03355    30 CONTINUE
03356       JZ = J1
03357       JU = 0
03358 !C
03359 !C    GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
03360 !C
03361       NM1 = N - 1
03362       IF (NM1 .LT. 1) GO TO 130
03363       DO 120 K = 1, NM1
03364          KP1 = K + 1
03365 !C
03366 !C       ZERO NEXT FILL-IN COLUMN
03367 !C
03368          JZ = JZ + 1
03369          IF (JZ .GT. N) GO TO 50
03370          IF (ML .LT. 1) GO TO 50
03371             DO 40 I = 1, ML
03372                ABD(I,JZ) = 0.0D0
03373    40       CONTINUE
03374    50    CONTINUE
03375 !C
03376 !C       FIND L = PIVOT INDEX
03377 !C
03378          LM = MIN(ML,N-K)
03379          L = IDAMAX(LM+1,ABD(M,K),1) + M - 1
03380          IPVT(K) = L + K - M
03381 !C
03382 !C       ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
03383 !C
03384          IF (ABD(L,K) .EQ. 0.0D0) GO TO 100
03385 !C
03386 !C          INTERCHANGE IF NECESSARY
03387 !C
03388             IF (L .EQ. M) GO TO 60
03389                T = ABD(L,K)
03390                ABD(L,K) = ABD(M,K)
03391                ABD(M,K) = T
03392    60       CONTINUE
03393 !C
03394 !C          COMPUTE MULTIPLIERS
03395 !C
03396             T = -1.0D0/ABD(M,K)
03397             CALL DSCAL(LM,T,ABD(M+1,K),1)
03398 !C
03399 !C          ROW ELIMINATION WITH COLUMN INDEXING
03400 !C
03401             JU = MIN(MAX(JU,MU+IPVT(K)),N)
03402             MM = M
03403             IF (JU .LT. KP1) GO TO 90
03404             DO 80 J = KP1, JU
03405                L = L - 1
03406                MM = MM - 1
03407                T = ABD(L,J)
03408                IF (L .EQ. MM) GO TO 70
03409                   ABD(L,J) = ABD(MM,J)
03410                   ABD(MM,J) = T
03411    70          CONTINUE
03412                CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
03413    80       CONTINUE
03414    90       CONTINUE
03415          GO TO 110
03416   100    CONTINUE
03417             INFO = K
03418   110    CONTINUE
03419   120 CONTINUE
03420   130 CONTINUE
03421       IPVT(N) = N
03422       IF (ABD(M,N) .EQ. 0.0D0) INFO = N
03423       RETURN
03424       END
03425 !*DECK DGBSL
03426       SUBROUTINE DGBSL (ABD, LDA, N, ML, MU, IPVT, B, JOB)
03427 !C**BEGIN PROLOGUE  DGBSL
03428 !C**PURPOSE  Solve the real band system A*X=B or TRANS(A)*X=B using
03429 !C           the factors computed by DGBCO or DGBFA.
03430 !C**CATEGORY  D2A2
03431 !C**TYPE      DOUBLE PRECISION (SGBSL-S, DGBSL-D, CGBSL-C)
03432 !C**KEYWORDS  BANDED, LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE
03433 !C**AUTHOR  Moler, C. B., (U. of New Mexico)
03434 !C**DESCRIPTION
03435 !C
03436 !C    DGBSL solves the double precision band system
03437 !C    A * X = B  or  TRANS(A) * X = B
03438 !C    using the factors computed by DGBCO or DGBFA.
03439 !C
03440 !C    On Entry
03441 !C
03442 !C       ABD     DOUBLE PRECISION(LDA, N)
03443 !C               the output from DGBCO or DGBFA.
03444 !C
03445 !C       LDA     INTEGER
03446 !C               the leading dimension of the array  ABD .
03447 !C
03448 !C       N       INTEGER
03449 !C               the order of the original matrix.
03450 !C
03451 !C       ML      INTEGER
03452 !C               number of diagonals below the main diagonal.
03453 !C
03454 !C       MU      INTEGER
03455 !C               number of diagonals above the main diagonal.
03456 !C
03457 !C       IPVT    INTEGER(N)
03458 !C               the pivot vector from DGBCO or DGBFA.
03459 !C
03460 !C       B       DOUBLE PRECISION(N)
03461 !C               the right hand side vector.
03462 !C
03463 !C       JOB     INTEGER
03464 !C               = 0         to solve  A*X = B ,
03465 !C               = nonzero   to solve  TRANS(A)*X = B , where
03466 !C                           TRANS(A)  is the transpose.
03467 !C
03468 !C    On Return
03469 !C
03470 !C       B       the solution vector  X .
03471 !C
03472 !C    Error Condition
03473 !C
03474 !C       A division by zero will occur if the input factor contains a
03475 !C       zero on the diagonal.  Technically this indicates singularity
03476 !C       but it is often caused by improper arguments or improper
03477 !C       setting of LDA .  It will not occur if the subroutines are
03478 !C       called correctly and if DGBCO has set RCOND .GT. 0.0
03479 !C       or DGBFA has set INFO .EQ. 0 .
03480 !C
03481 !C    To compute  INVERSE(A) * C  where  C  is a matrix
03482 !C    with  P  columns
03483 !C          CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
03484 !C          IF (RCOND is too small) GO TO ...
03485 !C          DO 10 J = 1, P
03486 !C             CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
03487 !C       10 CONTINUE
03488 !C
03489 !C**REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
03490 !C                Stewart, LINPACK Users' Guide, SIAM, 1979.
03491 !C**ROUTINES CALLED  DAXPY, DDOT
03492 !C**REVISION HISTORY  (YYMMDD)
03493 !C  780814  DATE WRITTEN
03494 !C  890531  Changed all specific intrinsics to generic.  (WRB)
03495 !C  890831  Modified array declarations.  (WRB)
03496 !C  890831  REVISION DATE from Version 3.2
03497 !C  891214  Prologue converted to Version 4.0 format.  (BAB)
03498 !C  900326  Removed duplicate information from DESCRIPTION section.
03499 !C          (WRB)
03500 !C  920501  Reformatted the REFERENCES section.  (WRB)
03501 !C**END PROLOGUE  DGBSL
03502       INTEGER LDA,N,ML,MU,IPVT(*),JOB
03503       DOUBLE PRECISION ABD(LDA,*),B(*)
03504 !C
03505       DOUBLE PRECISION DDOT,T
03506       INTEGER K,KB,L,LA,LB,LM,M,NM1
03507 !C**FIRST EXECUTABLE STATEMENT  DGBSL
03508       M = MU + ML + 1
03509       NM1 = N - 1
03510       IF (JOB .NE. 0) GO TO 50
03511 !C
03512 !C       JOB = 0 , SOLVE  A * X = B
03513 !C       FIRST SOLVE L*Y = B
03514 !C
03515          IF (ML .EQ. 0) GO TO 30
03516          IF (NM1 .LT. 1) GO TO 30
03517             DO 20 K = 1, NM1
03518                LM = MIN(ML,N-K)
03519                L = IPVT(K)
03520                T = B(L)
03521                IF (L .EQ. K) GO TO 10
03522                   B(L) = B(K)
03523                   B(K) = T
03524    10          CONTINUE
03525                CALL DAXPY(LM,T,ABD(M+1,K),1,B(K+1),1)
03526    20       CONTINUE
03527    30    CONTINUE
03528 !C
03529 !C       NOW SOLVE  U*X = Y
03530 !C
03531          DO 40 KB = 1, N
03532             K = N + 1 - KB
03533             B(K) = B(K)/ABD(M,K)
03534             LM = MIN(K,M) - 1
03535             LA = M - LM
03536             LB = K - LM
03537             T = -B(K)
03538             CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1)
03539    40    CONTINUE
03540       GO TO 100
03541    50 CONTINUE
03542 !C
03543 !C       JOB = NONZERO, SOLVE  TRANS(A) * X = B
03544 !C       FIRST SOLVE  TRANS(U)*Y = B
03545 !C
03546          DO 60 K = 1, N
03547             LM = MIN(K,M) - 1
03548             LA = M - LM
03549             LB = K - LM
03550             T = DDOT(LM,ABD(LA,K),1,B(LB),1)
03551             B(K) = (B(K) - T)/ABD(M,K)
03552    60    CONTINUE
03553 !C
03554 !C       NOW SOLVE TRANS(L)*X = Y
03555 !C
03556          IF (ML .EQ. 0) GO TO 90
03557          IF (NM1 .LT. 1) GO TO 90
03558             DO 80 KB = 1, NM1
03559                K = N - KB
03560                LM = MIN(ML,N-K)
03561                B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1)
03562                L = IPVT(K)
03563                IF (L .EQ. K) GO TO 70
03564                   T = B(L)
03565                   B(L) = B(K)
03566                   B(K) = T
03567    70          CONTINUE
03568    80       CONTINUE
03569    90    CONTINUE
03570   100 CONTINUE
03571       RETURN
03572       END
03573 !*DECK DAXPY
03574       SUBROUTINE DAXPY (N, DA, DX, INCX, DY, INCY)
03575 !C**BEGIN PROLOGUE  DAXPY
03576 !C**PURPOSE  Compute a constant times a vector plus a vector.
03577 !C**CATEGORY  D1A7
03578 !C**TYPE      DOUBLE PRECISION (SAXPY-S, DAXPY-D, CAXPY-C)
03579 !C**KEYWORDS  BLAS, LINEAR ALGEBRA, TRIAD, VECTOR
03580 !C**AUTHOR  Lawson, C. L., (JPL)
03581 !C          Hanson, R. J., (SNLA)
03582 !C          Kincaid, D. R., (U. of Texas)
03583 !C          Krogh, F. T., (JPL)
03584 !C**DESCRIPTION
03585 !C
03586 !C               B L A S  Subprogram
03587 !C   Description of Parameters
03588 !C
03589 !C    --Input--
03590 !C       N  number of elements in input vector(s)
03591 !C      DA  double precision scalar multiplier
03592 !C      DX  double precision vector with N elements
03593 !C    INCX  storage spacing between elements of DX
03594 !C      DY  double precision vector with N elements
03595 !C    INCY  storage spacing between elements of DY
03596 !C
03597 !C    --Output--
03598 !C      DY  double precision result (unchanged if N .LE. 0)
03599 !C
03600 !C    Overwrite double precision DY with double precision DA*DX + DY.
03601 !C    For I = 0 to N-1, replace  DY(LY+I*INCY) with DA*DX(LX+I*INCX) +
03602 !C      DY(LY+I*INCY),
03603 !C    where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
03604 !C    defined in a similar way using INCY.
03605 !C
03606 !C**REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
03607 !C                Krogh, Basic linear algebra subprograms for Fortran
03608 !C                usage, Algorithm No. 539, Transactions on Mathematical
03609 !C                Software 5, 3 (September 1979), pp. 308-323.
03610 !C**ROUTINES CALLED  (NONE)
03611 !C**REVISION HISTORY  (YYMMDD)
03612 !C  791001  DATE WRITTEN
03613 !C  890831  Modified array declarations.  (WRB)
03614 !C  890831  REVISION DATE from Version 3.2
03615 !C  891214  Prologue converted to Version 4.0 format.  (BAB)
03616 !C  920310  Corrected definition of LX in DESCRIPTION.  (WRB)
03617 !C  920501  Reformatted the REFERENCES section.  (WRB)
03618 !C**END PROLOGUE  DAXPY
03619       DOUBLE PRECISION DX(*), DY(*), DA
03620 !C**FIRST EXECUTABLE STATEMENT  DAXPY
03621       IF (N.LE.0 .OR. DA.EQ.0.0D0) RETURN
03622       IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60
03623 !C
03624 !C    Code for unequal or nonpositive increments.
03625 !C
03626     5 IX = 1
03627       IY = 1
03628       IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
03629       IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
03630       DO 10 I = 1,N
03631         DY(IY) = DY(IY) + DA*DX(IX)
03632         IX = IX + INCX
03633         IY = IY + INCY
03634    10 CONTINUE
03635       RETURN
03636 !C
03637 !C    Code for both increments equal to 1.
03638 !C
03639 !C    Clean-up loop so remaining vector length is a multiple of 4.
03640 !C
03641    20 M = MOD(N,4)
03642       IF (M .EQ. 0) GO TO 40
03643       DO 30 I = 1,M
03644         DY(I) = DY(I) + DA*DX(I)
03645    30 CONTINUE
03646       IF (N .LT. 4) RETURN
03647    40 MP1 = M + 1
03648       DO 50 I = MP1,N,4
03649         DY(I) = DY(I) + DA*DX(I)
03650         DY(I+1) = DY(I+1) + DA*DX(I+1)
03651         DY(I+2) = DY(I+2) + DA*DX(I+2)
03652         DY(I+3) = DY(I+3) + DA*DX(I+3)
03653    50 CONTINUE
03654       RETURN
03655 !C
03656 !C    Code for equal, positive, non-unit increments.
03657 !C
03658    60 NS = N*INCX
03659       DO 70 I = 1,NS,INCX
03660         DY(I) = DA*DX(I) + DY(I)
03661    70 CONTINUE
03662       RETURN
03663       END
03664 !*DECK DDOT
03665       DOUBLE PRECISION FUNCTION DDOT (N, DX, INCX, DY, INCY)
03666 !C**BEGIN PROLOGUE  DDOT
03667 !C**PURPOSE  Compute the inner product of two vectors.
03668 !C**CATEGORY  D1A4
03669 !C**TYPE      DOUBLE PRECISION (SDOT-S, DDOT-D, CDOTU-C)
03670 !C**KEYWORDS  BLAS, INNER PRODUCT, LINEAR ALGEBRA, VECTOR
03671 !C**AUTHOR  Lawson, C. L., (JPL)
03672 !C          Hanson, R. J., (SNLA)
03673 !C          Kincaid, D. R., (U. of Texas)
03674 !C          Krogh, F. T., (JPL)
03675 !C**DESCRIPTION
03676 !C
03677 !C               B L A S  Subprogram
03678 !C   Description of Parameters
03679 !C
03680 !C    --Input--
03681 !C       N  number of elements in input vector(s)
03682 !C      DX  double precision vector with N elements
03683 !C    INCX  storage spacing between elements of DX
03684 !C      DY  double precision vector with N elements
03685 !C    INCY  storage spacing between elements of DY
03686 !C
03687 !C    --Output--
03688 !C    DDOT  double precision dot product (zero if N .LE. 0)
03689 !C
03690 !C    Returns the dot product of double precision DX and DY.
03691 !C    DDOT = sum for I = 0 to N-1 of  DX(LX+I*INCX) * DY(LY+I*INCY),
03692 !C    where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
03693 !C    defined in a similar way using INCY.
03694 !C
03695 !C**REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
03696 !C                Krogh, Basic linear algebra subprograms for Fortran
03697 !C                usage, Algorithm No. 539, Transactions on Mathematical
03698 !C                Software 5, 3 (September 1979), pp. 308-323.
03699 !C**ROUTINES CALLED  (NONE)
03700 !C**REVISION HISTORY  (YYMMDD)
03701 !C  791001  DATE WRITTEN
03702 !C  890831  Modified array declarations.  (WRB)
03703 !C  890831  REVISION DATE from Version 3.2
03704 !C  891214  Prologue converted to Version 4.0 format.  (BAB)
03705 !C  920310  Corrected definition of LX in DESCRIPTION.  (WRB)
03706 !C  920501  Reformatted the REFERENCES section.  (WRB)
03707 !C**END PROLOGUE  DDOT
03708       DOUBLE PRECISION DX(*), DY(*)
03709 !C**FIRST EXECUTABLE STATEMENT  DDOT
03710       DDOT = 0.0D0
03711       IF (N .LE. 0) RETURN
03712       IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60
03713 !C
03714 !C    Code for unequal or nonpositive increments.
03715 !C
03716     5 IX = 1
03717       IY = 1
03718       IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
03719       IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
03720       DO 10 I = 1,N
03721         DDOT = DDOT + DX(IX)*DY(IY)
03722         IX = IX + INCX
03723         IY = IY + INCY
03724    10 CONTINUE
03725       RETURN
03726 !C
03727 !C    Code for both increments equal to 1.
03728 !C
03729 !C    Clean-up loop so remaining vector length is a multiple of 5.
03730 !C
03731    20 M = MOD(N,5)
03732       IF (M .EQ. 0) GO TO 40
03733       DO 30 I = 1,M
03734          DDOT = DDOT + DX(I)*DY(I)
03735    30 CONTINUE
03736       IF (N .LT. 5) RETURN
03737    40 MP1 = M + 1
03738       DO 50 I = MP1,N,5
03739       DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) + DX(I+2)*DY(I+2) + &
03740                    DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4)
03741    50 CONTINUE
03742       RETURN
03743 !C
03744 !C    Code for equal, positive, non-unit increments.
03745 !C
03746    60 NS = N*INCX
03747       DO 70 I = 1,NS,INCX
03748         DDOT = DDOT + DX(I)*DY(I)
03749    70 CONTINUE
03750       RETURN
03751       END
03752 !*DECK DSCAL
03753       SUBROUTINE DSCAL (N, DA, DX, INCX)
03754 !C**BEGIN PROLOGUE  DSCAL
03755 !C**PURPOSE  Multiply a vector by a constant.
03756 !C**CATEGORY  D1A6
03757 !C**TYPE      DOUBLE PRECISION (SSCAL-S, DSCAL-D, CSCAL-C)
03758 !C**KEYWORDS  BLAS, LINEAR ALGEBRA, SCALE, VECTOR
03759 !C**AUTHOR  Lawson, C. L., (JPL)
03760 !C          Hanson, R. J., (SNLA)
03761 !C          Kincaid, D. R., (U. of Texas)
03762 !C          Krogh, F. T., (JPL)
03763 !C**DESCRIPTION
03764 !C
03765 !C               B L A S  Subprogram
03766 !C   Description of Parameters
03767 !C
03768 !C    --Input--
03769 !C       N  number of elements in input vector(s)
03770 !C      DA  double precision scale factor
03771 !C      DX  double precision vector with N elements
03772 !C    INCX  storage spacing between elements of DX
03773 !C
03774 !C    --Output--
03775 !C      DX  double precision result (unchanged if N.LE.0)
03776 !C
03777 !C    Replace double precision DX by double precision DA*DX.
03778 !C    For I = 0 to N-1, replace DX(IX+I*INCX) with  DA * DX(IX+I*INCX),
03779 !C    where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX.
03780 !C
03781 !C**REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
03782 !C                Krogh, Basic linear algebra subprograms for Fortran
03783 !C                usage, Algorithm No. 539, Transactions on Mathematical
03784 !C                Software 5, 3 (September 1979), pp. 308-323.
03785 !C**ROUTINES CALLED  (NONE)
03786 !C**REVISION HISTORY  (YYMMDD)
03787 !C  791001  DATE WRITTEN
03788 !C  890831  Modified array declarations.  (WRB)
03789 !C  890831  REVISION DATE from Version 3.2
03790 !C  891214  Prologue converted to Version 4.0 format.  (BAB)
03791 !C  900821  Modified to correct problem with a negative increment.
03792 !C          (WRB)
03793 !C  920501  Reformatted the REFERENCES section.  (WRB)
03794 !C**END PROLOGUE  DSCAL
03795       DOUBLE PRECISION DA, DX(*)
03796       INTEGER I, INCX, IX, M, MP1, N
03797 !C**FIRST EXECUTABLE STATEMENT  DSCAL
03798       IF (N .LE. 0) RETURN
03799       IF (INCX .EQ. 1) GOTO 20
03800 !C
03801 !C    Code for increment not equal to 1.
03802 !C
03803       IX = 1
03804       IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
03805       DO 10 I = 1,N
03806         DX(IX) = DA*DX(IX)
03807         IX = IX + INCX
03808    10 CONTINUE
03809       RETURN
03810 !C
03811 !C    Code for increment equal to 1.
03812 !C
03813 !C    Clean-up loop so remaining vector length is a multiple of 5.
03814 !C
03815    20 M = MOD(N,5)
03816       IF (M .EQ. 0) GOTO 40
03817       DO 30 I = 1,M
03818         DX(I) = DA*DX(I)
03819    30 CONTINUE
03820       IF (N .LT. 5) RETURN
03821    40 MP1 = M + 1
03822       DO 50 I = MP1,N,5
03823         DX(I) = DA*DX(I)
03824         DX(I+1) = DA*DX(I+1)
03825         DX(I+2) = DA*DX(I+2)
03826         DX(I+3) = DA*DX(I+3)
03827         DX(I+4) = DA*DX(I+4)
03828    50 CONTINUE
03829       RETURN
03830       END
03831 !*DECK IDAMAX
03832       INTEGER FUNCTION IDAMAX (N, DX, INCX)
03833 !C**BEGIN PROLOGUE  IDAMAX
03834 !C**PURPOSE  Find the smallest index of that component of a vector
03835 !C           having the maximum magnitude.
03836 !C**CATEGORY  D1A2
03837 !C**TYPE      DOUBLE PRECISION (ISAMAX-S, IDAMAX-D, ICAMAX-C)
03838 !C**KEYWORDS  BLAS, LINEAR ALGEBRA, MAXIMUM COMPONENT, VECTOR
03839 !C**AUTHOR  Lawson, C. L., (JPL)
03840 !C          Hanson, R. J., (SNLA)
03841 !C          Kincaid, D. R., (U. of Texas)
03842 !C          Krogh, F. T., (JPL)
03843 !C**DESCRIPTION
03844 !C
03845 !C               B L A S  Subprogram
03846 !C   Description of Parameters
03847 !C
03848 !C    --Input--
03849 !C       N  number of elements in input vector(s)
03850 !C      DX  double precision vector with N elements
03851 !C    INCX  storage spacing between elements of DX
03852 !C
03853 !C    --Output--
03854 !C  IDAMAX  smallest index (zero if N .LE. 0)
03855 !C
03856 !C    Find smallest index of maximum magnitude of double precision DX.
03857 !C    IDAMAX = first I, I = 1 to N, to maximize ABS(DX(IX+(I-1)*INCX)),
03858 !C    where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX.
03859 !C
03860 !C**REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
03861 !C                Krogh, Basic linear algebra subprograms for Fortran
03862 !C                usage, Algorithm No. 539, Transactions on Mathematical
03863 !C                Software 5, 3 (September 1979), pp. 308-323.
03864 !C**ROUTINES CALLED  (NONE)
03865 !C**REVISION HISTORY  (YYMMDD)
03866 !C  791001  DATE WRITTEN
03867 !C  890531  Changed all specific intrinsics to generic.  (WRB)
03868 !C  890531  REVISION DATE from Version 3.2
03869 !C  891214  Prologue converted to Version 4.0 format.  (BAB)
03870 !C  900821  Modified to correct problem with a negative increment.
03871 !C          (WRB)
03872 !C  920501  Reformatted the REFERENCES section.  (WRB)
03873 !C**END PROLOGUE  IDAMAX
03874       DOUBLE PRECISION DX(*), DMAX, XMAG
03875       INTEGER I, INCX, IX, N
03876 !C**FIRST EXECUTABLE STATEMENT  IDAMAX
03877       IDAMAX = 0
03878       IF (N .LE. 0) RETURN
03879       IDAMAX = 1
03880       IF (N .EQ. 1) RETURN
03881 !C
03882       IF (INCX .EQ. 1) GOTO 20
03883 !C
03884 !C    Code for increments not equal to 1.
03885 !C
03886       IX = 1
03887       IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
03888       DMAX = DABS(DX(IX))
03889       IX = IX + INCX
03890       DO 10 I = 2,N
03891         XMAG = DABS(DX(IX))
03892         IF (XMAG .GT. DMAX) THEN
03893           IDAMAX = I
03894           DMAX = XMAG
03895         ENDIF
03896         IX = IX + INCX
03897    10 CONTINUE
03898       RETURN
03899 !C
03900 !C    Code for increments equal to 1.
03901 !C
03902    20 DMAX = DABS(DX(1))
03903       DO 30 I = 2,N
03904         XMAG = DABS(DX(I))
03905         IF (XMAG .GT. DMAX) THEN
03906           IDAMAX = I
03907           DMAX = XMAG
03908         ENDIF
03909    30 CONTINUE
03910       RETURN
03911       END
03912 !*DECK XERRWD
03913       SUBROUTINE XERRWD (MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
03914 !C**BEGIN PROLOGUE  XERRWD
03915 !C**SUBSIDIARY
03916 !C**PURPOSE  Write error message with values.
03917 !C**CATEGORY  R3C
03918 !C**TYPE      DOUBLE PRECISION (XERRWV-S, XERRWD-D)
03919 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
03920 !C**DESCRIPTION
03921 !C
03922 !C Subroutines XERRWD, XSETF, XSETUN, and the function routine IXSAV,
03923 !C as given here, constitute a simplified version of the SLATEC error
03924 !C handling package.
03925 !C
03926 !C All arguments are input arguments.
03927 !C
03928 !C MSG    = The message (character array).
03929 !C NMES   = The length of MSG (number of characters).
03930 !C NERR   = The error number (not used).
03931 !C LEVEL  = The error level..
03932 !C          0 or 1 means recoverable (control returns to caller).
03933 !C          2 means fatal (run is aborted--see note below).
03934 !C NI     = Number of integers (0, 1, or 2) to be printed with message.
03935 !C I1,I2  = Integers to be printed, depending on NI.
03936 !C NR     = Number of reals (0, 1, or 2) to be printed with message.
03937 !C R1,R2  = Reals to be printed, depending on NR.
03938 !C
03939 !C Note..  this routine is machine-dependent and specialized for use
03940 !C in limited context, in the following ways..
03941 !C 1. The argument MSG is assumed to be of type CHARACTER, and
03942 !C    the message is printed with a format of (1X,A).
03943 !C 2. The message is assumed to take only one line.
03944 !C    Multi-line messages are generated by repeated calls.
03945 !C 3. If LEVEL = 2, control passes to the statement   STOP
03946 !C    to abort the run.  This statement may be machine-dependent.
03947 !C 4. R1 and R2 are assumed to be in double precision and are printed
03948 !C    in D21.13 format.
03949 !C
03950 !C**ROUTINES CALLED  IXSAV
03951 !C**REVISION HISTORY  (YYMMDD)
03952 !C  920831  DATE WRITTEN
03953 !C  921118  Replaced MFLGSV/LUNSAV by IXSAV. (ACH)
03954 !C  930329  Modified prologue to SLATEC format. (FNF)
03955 !C  930407  Changed MSG from CHARACTER*1 array to variable. (FNF)
03956 !C  930922  Minor cosmetic change. (FNF)
03957 !C**END PROLOGUE  XERRWD
03958 !C
03959 !CInternal Notes:
03960 !C
03961 !CFor a different default logical unit number, IXSAV (or a subsidiary
03962 !Croutine that it calls) will need to be modified.
03963 !CFor a different run-abort command, change the statement following
03964 !Cstatement 100 at the end.
03965 !C----------------------------------------------------------------------
03966 !CSubroutines called by XERRWD.. None
03967 !CFunction routine called by XERRWD.. IXSAV
03968 !C----------------------------------------------------------------------
03969 !C*End
03970 !C
03971 !C Declare arguments.
03972 !C
03973       DOUBLE PRECISION R1, R2
03974       INTEGER NMES, NERR, LEVEL, NI, I1, I2, NR
03975       CHARACTER*(*) MSG
03976 !C
03977 !C Declare local variables.
03978 !C
03979       INTEGER LUNIT, IXSAV, MESFLG
03980 !C
03981 !C Get logical unit number and message print flag.
03982 !C
03983 !C**FIRST EXECUTABLE STATEMENT  XERRWD
03984       LUNIT = IXSAV (1, 0, .FALSE.)
03985       MESFLG = IXSAV (2, 0, .FALSE.)
03986       IF (MESFLG .EQ. 0) GO TO 100
03987 !C
03988 !C Write the message.
03989 !C
03990 !      WRITE (LUNIT,10)  MSG
03991 ! 10   FORMAT(1X,A)
03992 !      IF (NI .EQ. 1) WRITE (LUNIT, 20) I1
03993 ! 20   FORMAT(6X,'In above message,  I1 =',I10)
03994 !      IF (NI .EQ. 2) WRITE (LUNIT, 30) I1,I2
03995 ! 30   FORMAT(6X,'In above message,  I1 =',I10,3X,'I2 =',I10)
03996 !      IF (NR .EQ. 1) WRITE (LUNIT, 40) R1
03997 ! 40   FORMAT(6X,'In above message,  R1 =',D21.13)
03998 !      IF (NR .EQ. 2) WRITE (LUNIT, 50) R1,R2
03999 ! 50   FORMAT(6X,'In above,  R1 =',D21.13,3X,'R2 =',D21.13)
04000 !C
04001 !C Abort the run if LEVEL = 2.
04002 !C
04003  100  IF (LEVEL .NE. 2) RETURN
04004 !Cmodif 07/10/06 par guedj
04005 !C     STOP
04006 !C---------------------- End of Subroutine XERRWD ----------------------
04007       END
04008 !*DECK XSETF
04009       SUBROUTINE XSETF (MFLAG)
04010 !C**BEGIN PROLOGUE  XSETF
04011 !C**PURPOSE  Reset the error print control flag.
04012 !C**CATEGORY  R3A
04013 !C**TYPE      ALL (XSETF-A)
04014 !C**KEYWORDS  ERROR CONTROL
04015 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
04016 !C**DESCRIPTION
04017 !C
04018 !C  XSETF sets the error print control flag to MFLAG:
04019 !C     MFLAG=1 means print all messages (the default).
04020 !C     MFLAG=0 means no printing.
04021 !C
04022 !C**SEE ALSO  XERRWD, XERRWV
04023 !C**REFERENCES  (NONE)
04024 !C**ROUTINES CALLED  IXSAV
04025 !C**REVISION HISTORY  (YYMMDD)
04026 !C  921118  DATE WRITTEN
04027 !C  930329  Added SLATEC format prologue. (FNF)
04028 !C  930407  Corrected SEE ALSO section. (FNF)
04029 !C  930922  Made user-callable, and other cosmetic changes. (FNF)
04030 !C**END PROLOGUE  XSETF
04031 !C
04032 !CSubroutines called by XSETF.. None
04033 !CFunction routine called by XSETF.. IXSAV
04034 !C----------------------------------------------------------------------
04035 !C*End
04036       INTEGER MFLAG, JUNK, IXSAV
04037 !C
04038 !C**FIRST EXECUTABLE STATEMENT  XSETF
04039       IF (MFLAG .EQ. 0 .OR. MFLAG .EQ. 1) JUNK = IXSAV (2,MFLAG,.TRUE.)
04040       RETURN
04041 !C---------------------- End of Subroutine XSETF -----------------------
04042       END
04043 !*DECK XSETUN
04044       SUBROUTINE XSETUN (LUN)
04045 !C**BEGIN PROLOGUE  XSETUN
04046 !C**PURPOSE  Reset the logical unit number for error messages.
04047 !C**CATEGORY  R3B
04048 !C**TYPE      ALL (XSETUN-A)
04049 !C**KEYWORDS  ERROR CONTROL
04050 !C**DESCRIPTION
04051 !C
04052 !C  XSETUN sets the logical unit number for error messages to LUN.
04053 !C
04054 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
04055 !C**SEE ALSO  XERRWD, XERRWV
04056 !C**REFERENCES  (NONE)
04057 !C**ROUTINES CALLED  IXSAV
04058 !C**REVISION HISTORY  (YYMMDD)
04059 !C  921118  DATE WRITTEN
04060 !C  930329  Added SLATEC format prologue. (FNF)
04061 !C  930407  Corrected SEE ALSO section. (FNF)
04062 !C  930922  Made user-callable, and other cosmetic changes. (FNF)
04063 !C**END PROLOGUE  XSETUN
04064 !C
04065 !CSubroutines called by XSETUN.. None
04066 !CFunction routine called by XSETUN.. IXSAV
04067 !C----------------------------------------------------------------------
04068 !C*End
04069       INTEGER LUN, JUNK, IXSAV
04070 !C
04071 !C**FIRST EXECUTABLE STATEMENT  XSETUN
04072       IF (LUN .GT. 0) JUNK = IXSAV (1,LUN,.TRUE.)
04073       RETURN
04074 !C---------------------- End of Subroutine XSETUN ----------------------
04075       END
04076 !*DECK IXSAV
04077       INTEGER FUNCTION IXSAV (IPAR, IVALUE, ISET)
04078 !C**BEGIN PROLOGUE  IXSAV
04079 !C**SUBSIDIARY
04080 !C**PURPOSE  Save and recall error message control parameters.
04081 !C**CATEGORY  R3C
04082 !C**TYPE      ALL (IXSAV-A)
04083 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
04084 !C**DESCRIPTION
04085 !C
04086 !C IXSAV saves and recalls one of two error message parameters:
04087 !C   LUNIT, the logical unit number to which messages are printed, and
04088 !C   MESFLG, the message print flag.
04089 !C This is a modification of the SLATEC library routine J4SAVE.
04090 !C
04091 !C Saved local variables..
04092 !C  LUNIT  = Logical unit number for messages.  The default is obtained
04093 !C           by a call to IUMACH (may be machine-dependent).
04094 !C  MESFLG = Print control flag..
04095 !C           1 means print all messages (the default).
04096 !C           0 means no printing.
04097 !C
04098 !C On input..
04099 !C   IPAR   = Parameter indicator (1 for LUNIT, 2 for MESFLG).
04100 !C   IVALUE = The value to be set for the parameter, if ISET = .TRUE.
04101 !C   ISET   = Logical flag to indicate whether to read or write.
04102 !C            If ISET = .TRUE., the parameter will be given
04103 !C            the value IVALUE.  If ISET = .FALSE., the parameter
04104 !C            will be unchanged, and IVALUE is a dummy argument.
04105 !C
04106 !C On return..
04107 !C   IXSAV = The (old) value of the parameter.
04108 !C
04109 !C**SEE ALSO  XERRWD, XERRWV
04110 !C**ROUTINES CALLED  IUMACH
04111 !C**REVISION HISTORY  (YYMMDD)
04112 !C  921118  DATE WRITTEN
04113 !C  930329  Modified prologue to SLATEC format. (FNF)
04114 !C  930915  Added IUMACH call to get default output unit.  (ACH)
04115 !C  930922  Minor cosmetic changes. (FNF)
04116 !C  010425  Type declaration for IUMACH added. (ACH)
04117 !C**END PROLOGUE  IXSAV
04118 !C
04119 !CSubroutines called by IXSAV.. None
04120 !CFunction routine called by IXSAV.. IUMACH
04121 !C----------------------------------------------------------------------
04122 !C*End
04123       LOGICAL ISET
04124       INTEGER IPAR, IVALUE
04125 !C----------------------------------------------------------------------
04126       INTEGER IUMACH, LUNIT, MESFLG
04127 !C----------------------------------------------------------------------
04128 !CThe following Fortran-77 declaration is to cause the values of the
04129 !Clisted (local) variables to be saved between calls to this routine.
04130 !C----------------------------------------------------------------------
04131       SAVE LUNIT, MESFLG
04132       DATA LUNIT/-1/, MESFLG/1/
04133 !C
04134 !C**FIRST EXECUTABLE STATEMENT  IXSAV
04135       IF (IPAR .EQ. 1) THEN
04136         IF (LUNIT .EQ. -1) LUNIT = IUMACH()
04137         IXSAV = LUNIT
04138         IF (ISET) LUNIT = IVALUE
04139         ENDIF
04140 !C
04141       IF (IPAR .EQ. 2) THEN
04142         IXSAV = MESFLG
04143         IF (ISET) MESFLG = IVALUE
04144         ENDIF
04145 !C
04146       RETURN
04147 !C---------------------- End of Function IXSAV -------------------------
04148       END
04149 !*DECK IUMACH
04150       INTEGER FUNCTION IUMACH()
04151 !C**BEGIN PROLOGUE  IUMACH
04152 !C**PURPOSE  Provide standard output unit number.
04153 !C**CATEGORY  R1
04154 !C**TYPE      INTEGER (IUMACH-I)
04155 !C**KEYWORDS  MACHINE CONSTANTS
04156 !C**AUTHOR  Hindmarsh, Alan C., (LLNL)
04157 !C**DESCRIPTION
04158 !C*Usage:
04159 !C       INTEGER  LOUT, IUMACH
04160 !C       LOUT = IUMACH()
04161 !C
04162 !C*Function Return Values:
04163 !C    LOUT : the standard logical unit for Fortran output.
04164 !C
04165 !C**REFERENCES  (NONE)
04166 !C**ROUTINES CALLED  (NONE)
04167 !C**REVISION HISTORY  (YYMMDD)
04168 !C  930915  DATE WRITTEN
04169 !C  930922  Made user-callable, and other cosmetic changes. (FNF)
04170 !C**END PROLOGUE  IUMACH
04171 !C
04172 !CInternal Notes:
04173 !C The built-in value of 6 is standard on a wide range of Fortran
04174 !C systems.  This may be machine-dependent.
04175 !C*End
04176 !C**FIRST EXECUTABLE STATEMENT  IUMACH
04177       IUMACH = 6
04178 !C
04179       RETURN
04180 !C---------------------- End of Function IUMACH ------------------------
04181       END