ddaspk.f

Go to the documentation of this file.
00001 C Work performed under the auspices of the U.S. Department of Energy
00002 C by Lawrence Livermore National Laboratory under contract number 
00003 C W-7405-Eng-48.
00004 C
00005       SUBROUTINE DDASPK (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
00006      *   IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL)
00007 C
00008 C***BEGIN PROLOGUE  DDASPK
00009 C***DATE WRITTEN   890101   (YYMMDD)
00010 C***REVISION DATE  910624   
00011 C***REVISION DATE  920929   (CJ in RES call, RES counter fix.)
00012 C***REVISION DATE  921215   (Warnings on poor iteration performance)
00013 C***REVISION DATE  921216   (NRMAX as optional input)
00014 C***REVISION DATE  930315   (Name change: DDINI to DDINIT)
00015 C***REVISION DATE  940822   (Replaced initial condition calculation)
00016 C***REVISION DATE  941101   (Added linesearch in I.C. calculations)
00017 C***REVISION DATE  941220   (Misc. corrections throughout)
00018 C***REVISION DATE  950125   (Added DINVWT routine)
00019 C***REVISION DATE  950714   (Misc. corrections throughout)
00020 C***REVISION DATE  950802   (Default NRMAX = 5, based on tests.)
00021 C***REVISION DATE  950808   (Optional error test added.)
00022 C***REVISION DATE  950814   (Added I.C. constraints and INFO(14))
00023 C***REVISION DATE  950828   (Various minor corrections.)
00024 C***REVISION DATE  951006   (Corrected WT scaling in DFNRMK.)
00025 C***REVISION DATE  960129   (Corrected RL bug in DLINSD, DLINSK.)
00026 C***REVISION DATE  960301   (Added NONNEG to SAVE statement.)
00027 C***CATEGORY NO.  I1A2
00028 C***KEYWORDS  DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS,
00029 C             IMPLICIT DIFFERENTIAL SYSTEMS, KRYLOV ITERATION
00030 C***AUTHORS   Linda R. Petzold, Peter N. Brown, Alan C. Hindmarsh, and
00031 C                  Clement W. Ulrich
00032 C             Center for Computational Sciences & Engineering, L-316
00033 C             Lawrence Livermore National Laboratory
00034 C             P.O. Box 808,
00035 C             Livermore, CA 94551
00036 C***PURPOSE  This code solves a system of differential/algebraic 
00037 C            equations of the form 
00038 C               G(t,y,y') = 0 , 
00039 C            using a combination of Backward Differentiation Formula 
00040 C            (BDF) methods and a choice of two linear system solution 
00041 C            methods: direct (dense or band) or Krylov (iterative).
00042 C            This version is in double precision.
00043 C-----------------------------------------------------------------------
00044 C***DESCRIPTION
00045 C
00046 C *Usage:
00047 C
00048 C      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
00049 C      INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR(*)
00050 C      DOUBLE PRECISION T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*),
00051 C         RWORK(LRW), RPAR(*)
00052 C      EXTERNAL  RES, JAC, PSOL
00053 C
00054 C      CALL DDASPK (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
00055 C     *   IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL)
00056 C
00057 C  Quantities which may be altered by the code are:
00058 C     T, Y(*), YPRIME(*), INFO(*), RTOL, ATOL, IDID, RWORK(*), IWORK(*)
00059 C
00060 C
00061 C *Arguments:
00062 C
00063 C  RES:EXT          This is the name of a subroutine which you
00064 C                   provide to define the residual function G(t,y,y')
00065 C                   of the differential/algebraic system.
00066 C
00067 C  NEQ:IN           This is the number of equations in the system.
00068 C
00069 C  T:INOUT          This is the current value of the independent 
00070 C                   variable.
00071 C
00072 C  Y(*):INOUT       This array contains the solution components at T.
00073 C
00074 C  YPRIME(*):INOUT  This array contains the derivatives of the solution
00075 C                   components at T.
00076 C
00077 C  TOUT:IN          This is a point at which a solution is desired.
00078 C
00079 C  INFO(N):IN       This is an integer array used to communicate details
00080 C                   of how the solution is to be carried out, such as
00081 C                   tolerance type, matrix structure, step size and
00082 C                   order limits, and choice of nonlinear system method.
00083 C                   N must be at least 20.
00084 C
00085 C  RTOL,ATOL:INOUT  These quantities represent absolute and relative
00086 C                   error tolerances (on local error) which you provide
00087 C                   to indicate how accurately you wish the solution to
00088 C                   be computed.  You may choose them to be both scalars
00089 C                   or else both arrays of length NEQ.
00090 C
00091 C  IDID:OUT         This integer scalar is an indicator reporting what
00092 C                   the code did.  You must monitor this variable to
00093 C                   decide what action to take next.
00094 C
00095 C  RWORK:WORK       A real work array of length LRW which provides the
00096 C                   code with needed storage space.
00097 C
00098 C  LRW:IN           The length of RWORK.
00099 C
00100 C  IWORK:WORK       An integer work array of length LIW which provides
00101 C                   the code with needed storage space.
00102 C
00103 C  LIW:IN           The length of IWORK.
00104 C
00105 C  RPAR,IPAR:IN     These are real and integer parameter arrays which
00106 C                   you can use for communication between your calling
00107 C                   program and the RES, JAC, and PSOL subroutines.
00108 C
00109 C  JAC:EXT          This is the name of a subroutine which you may
00110 C                   provide (optionally) for calculating Jacobian 
00111 C                   (partial derivative) data involved in solving linear
00112 C                   systems within DDASPK.
00113 C
00114 C  PSOL:EXT         This is the name of a subroutine which you must
00115 C                   provide for solving linear systems if you selected
00116 C                   a Krylov method.  The purpose of PSOL is to solve
00117 C                   linear systems involving a left preconditioner P.
00118 C
00119 C *Overview
00120 C
00121 C  The DDASPK solver uses the backward differentiation formulas of
00122 C  orders one through five to solve a system of the form G(t,y,y') = 0
00123 C  for y = Y and y' = YPRIME.  Values for Y and YPRIME at the initial 
00124 C  time must be given as input.  These values should be consistent, 
00125 C  that is, if T, Y, YPRIME are the given initial values, they should 
00126 C  satisfy G(T,Y,YPRIME) = 0.  However, if consistent values are not
00127 C  known, in many cases you can have DDASPK solve for them -- see INFO(11).
00128 C  (This and other options are described in more detail below.)
00129 C
00130 C  Normally, DDASPK solves the system from T to TOUT.  It is easy to
00131 C  continue the solution to get results at additional TOUT.  This is
00132 C  the interval mode of operation.  Intermediate results can also be
00133 C  obtained easily by specifying INFO(3).
00134 C
00135 C  On each step taken by DDASPK, a sequence of nonlinear algebraic  
00136 C  systems arises.  These are solved by one of two types of
00137 C  methods:
00138 C    * a Newton iteration with a direct method for the linear
00139 C      systems involved (INFO(12) = 0), or
00140 C    * a Newton iteration with a preconditioned Krylov iterative 
00141 C      method for the linear systems involved (INFO(12) = 1).
00142 C
00143 C  The direct method choices are dense and band matrix solvers, 
00144 C  with either a user-supplied or an internal difference quotient 
00145 C  Jacobian matrix, as specified by INFO(5) and INFO(6).
00146 C  In the band case, INFO(6) = 1, you must supply half-bandwidths
00147 C  in IWORK(1) and IWORK(2).
00148 C
00149 C  The Krylov method is the Generalized Minimum Residual (GMRES) 
00150 C  method, in either complete or incomplete form, and with 
00151 C  scaling and preconditioning.  The method is implemented
00152 C  in an algorithm called SPIGMR.  Certain options in the Krylov 
00153 C  method case are specified by INFO(13) and INFO(15).
00154 C
00155 C  If the Krylov method is chosen, you may supply a pair of routines,
00156 C  JAC and PSOL, to apply preconditioning to the linear system.
00157 C  If the system is A*x = b, the matrix is A = dG/dY + CJ*dG/dYPRIME
00158 C  (of order NEQ).  This system can then be preconditioned in the form
00159 C  (P-inverse)*A*x = (P-inverse)*b, with left preconditioner P.
00160 C  (DDASPK does not allow right preconditioning.)
00161 C  Then the Krylov method is applied to this altered, but equivalent,
00162 C  linear system, hopefully with much better performance than without
00163 C  preconditioning.  (In addition, a diagonal scaling matrix based on
00164 C  the tolerances is also introduced into the altered system.)
00165 C
00166 C  The JAC routine evaluates any data needed for solving systems
00167 C  with coefficient matrix P, and PSOL carries out that solution.
00168 C  In any case, in order to improve convergence, you should try to
00169 C  make P approximate the matrix A as much as possible, while keeping
00170 C  the system P*x = b reasonably easy and inexpensive to solve for x,
00171 C  given a vector b.
00172 C
00173 C
00174 C *Description
00175 C
00176 C------INPUT - WHAT TO DO ON THE FIRST CALL TO DDASPK-------------------
00177 C
00178 C
00179 C  The first call of the code is defined to be the start of each new
00180 C  problem.  Read through the descriptions of all the following items,
00181 C  provide sufficient storage space for designated arrays, set
00182 C  appropriate variables for the initialization of the problem, and
00183 C  give information about how you want the problem to be solved.
00184 C
00185 C
00186 C  RES -- Provide a subroutine of the form
00187 C
00188 C             SUBROUTINE RES (T, Y, YPRIME, CJ, DELTA, IRES, RPAR, IPAR)
00189 C
00190 C         to define the system of differential/algebraic
00191 C         equations which is to be solved. For the given values
00192 C         of T, Y and YPRIME, the subroutine should return
00193 C         the residual of the differential/algebraic system
00194 C             DELTA = G(T,Y,YPRIME)
00195 C         DELTA is a vector of length NEQ which is output from RES.
00196 C
00197 C         Subroutine RES must not alter T, Y, YPRIME, or CJ.
00198 C         You must declare the name RES in an EXTERNAL
00199 C         statement in your program that calls DDASPK.
00200 C         You must dimension Y, YPRIME, and DELTA in RES.
00201 C
00202 C         The input argument CJ can be ignored, or used to rescale
00203 C         constraint equations in the system (see Ref. 2, p. 145).
00204 C         Note: In this respect, DDASPK is not downward-compatible
00205 C         with DDASSL, which does not have the RES argument CJ.
00206 C
00207 C         IRES is an integer flag which is always equal to zero
00208 C         on input.  Subroutine RES should alter IRES only if it
00209 C         encounters an illegal value of Y or a stop condition.
00210 C         Set IRES = -1 if an input value is illegal, and DDASPK
00211 C         will try to solve the problem without getting IRES = -1.
00212 C         If IRES = -2, DDASPK will return control to the calling
00213 C         program with IDID = -11.
00214 C
00215 C         RPAR and IPAR are real and integer parameter arrays which
00216 C         you can use for communication between your calling program
00217 C         and subroutine RES. They are not altered by DDASPK. If you
00218 C         do not need RPAR or IPAR, ignore these parameters by treat-
00219 C         ing them as dummy arguments. If you do choose to use them,
00220 C         dimension them in your calling program and in RES as arrays
00221 C         of appropriate length.
00222 C
00223 C  NEQ -- Set it to the number of equations in the system (NEQ .GE. 1).
00224 C
00225 C  T -- Set it to the initial point of the integration. (T must be
00226 C       a variable.)
00227 C
00228 C  Y(*) -- Set this array to the initial values of the NEQ solution
00229 C          components at the initial point.  You must dimension Y of
00230 C          length at least NEQ in your calling program.
00231 C
00232 C  YPRIME(*) -- Set this array to the initial values of the NEQ first
00233 C               derivatives of the solution components at the initial
00234 C               point.  You must dimension YPRIME at least NEQ in your
00235 C               calling program. 
00236 C
00237 C  TOUT - Set it to the first point at which a solution is desired.
00238 C         You cannot take TOUT = T.  Integration either forward in T
00239 C         (TOUT .GT. T) or backward in T (TOUT .LT. T) is permitted.
00240 C
00241 C         The code advances the solution from T to TOUT using step
00242 C         sizes which are automatically selected so as to achieve the
00243 C         desired accuracy.  If you wish, the code will return with the
00244 C         solution and its derivative at intermediate steps (the
00245 C         intermediate-output mode) so that you can monitor them,
00246 C         but you still must provide TOUT in accord with the basic
00247 C         aim of the code.
00248 C
00249 C         The first step taken by the code is a critical one because
00250 C         it must reflect how fast the solution changes near the
00251 C         initial point.  The code automatically selects an initial
00252 C         step size which is practically always suitable for the
00253 C         problem.  By using the fact that the code will not step past
00254 C         TOUT in the first step, you could, if necessary, restrict the
00255 C         length of the initial step.
00256 C
00257 C         For some problems it may not be permissible to integrate
00258 C         past a point TSTOP, because a discontinuity occurs there
00259 C         or the solution or its derivative is not defined beyond
00260 C         TSTOP.  When you have declared a TSTOP point (see INFO(4)
00261 C         and RWORK(1)), you have told the code not to integrate past
00262 C         TSTOP.  In this case any tout beyond TSTOP is invalid input.
00263 C
00264 C  INFO(*) - Use the INFO array to give the code more details about
00265 C            how you want your problem solved.  This array should be
00266 C            dimensioned of length 20, though DDASPK uses only the 
00267 C            first 15 entries.  You must respond to all of the following
00268 C            items, which are arranged as questions.  The simplest use
00269 C            of DDASPK corresponds to setting all entries of INFO to 0.
00270 C
00271 C       INFO(1) - This parameter enables the code to initialize itself.
00272 C              You must set it to indicate the start of every new 
00273 C              problem.
00274 C
00275 C          **** Is this the first call for this problem ...
00276 C                yes - set INFO(1) = 0
00277 C                 no - not applicable here.
00278 C                      See below for continuation calls.  ****
00279 C
00280 C       INFO(2) - How much accuracy you want of your solution
00281 C              is specified by the error tolerances RTOL and ATOL.
00282 C              The simplest use is to take them both to be scalars.
00283 C              To obtain more flexibility, they can both be arrays.
00284 C              The code must be told your choice.
00285 C
00286 C          **** Are both error tolerances RTOL, ATOL scalars ...
00287 C                yes - set INFO(2) = 0
00288 C                      and input scalars for both RTOL and ATOL
00289 C                 no - set INFO(2) = 1
00290 C                      and input arrays for both RTOL and ATOL ****
00291 C
00292 C       INFO(3) - The code integrates from T in the direction of TOUT
00293 C              by steps.  If you wish, it will return the computed
00294 C              solution and derivative at the next intermediate step
00295 C              (the intermediate-output mode) or TOUT, whichever comes
00296 C              first.  This is a good way to proceed if you want to
00297 C              see the behavior of the solution.  If you must have
00298 C              solutions at a great many specific TOUT points, this
00299 C              code will compute them efficiently.
00300 C
00301 C          **** Do you want the solution only at
00302 C               TOUT (and not at the next intermediate step) ...
00303 C                yes - set INFO(3) = 0
00304 C                 no - set INFO(3) = 1 ****
00305 C
00306 C       INFO(4) - To handle solutions at a great many specific
00307 C              values TOUT efficiently, this code may integrate past
00308 C              TOUT and interpolate to obtain the result at TOUT.
00309 C              Sometimes it is not possible to integrate beyond some
00310 C              point TSTOP because the equation changes there or it is
00311 C              not defined past TSTOP.  Then you must tell the code
00312 C              this stop condition.
00313 C
00314 C           **** Can the integration be carried out without any
00315 C                restrictions on the independent variable T ...
00316 C                 yes - set INFO(4) = 0
00317 C                  no - set INFO(4) = 1
00318 C                       and define the stopping point TSTOP by
00319 C                       setting RWORK(1) = TSTOP ****
00320 C
00321 C       INFO(5) - used only when INFO(12) = 0 (direct methods).
00322 C              To solve differential/algebraic systems you may wish
00323 C              to use a matrix of partial derivatives of the
00324 C              system of differential equations.  If you do not
00325 C              provide a subroutine to evaluate it analytically (see
00326 C              description of the item JAC in the call list), it will
00327 C              be approximated by numerical differencing in this code.
00328 C              Although it is less trouble for you to have the code
00329 C              compute partial derivatives by numerical differencing,
00330 C              the solution will be more reliable if you provide the
00331 C              derivatives via JAC.  Usually numerical differencing is
00332 C              more costly than evaluating derivatives in JAC, but
00333 C              sometimes it is not - this depends on your problem.
00334 C
00335 C           **** Do you want the code to evaluate the partial deriv-
00336 C                atives automatically by numerical differences ...
00337 C                 yes - set INFO(5) = 0
00338 C                  no - set INFO(5) = 1
00339 C                       and provide subroutine JAC for evaluating the
00340 C                       matrix of partial derivatives ****
00341 C
00342 C       INFO(6) - used only when INFO(12) = 0 (direct methods).
00343 C              DDASPK will perform much better if the matrix of
00344 C              partial derivatives, dG/dY + CJ*dG/dYPRIME (here CJ is
00345 C              a scalar determined by DDASPK), is banded and the code
00346 C              is told this.  In this case, the storage needed will be
00347 C              greatly reduced, numerical differencing will be performed
00348 C              much cheaper, and a number of important algorithms will
00349 C              execute much faster.  The differential equation is said 
00350 C              to have half-bandwidths ML (lower) and MU (upper) if 
00351 C              equation i involves only unknowns Y(j) with
00352 C                             i-ML .le. j .le. i+MU .
00353 C              For all i=1,2,...,NEQ.  Thus, ML and MU are the widths
00354 C              of the lower and upper parts of the band, respectively,
00355 C              with the main diagonal being excluded.  If you do not
00356 C              indicate that the equation has a banded matrix of partial
00357 C              derivatives the code works with a full matrix of NEQ**2
00358 C              elements (stored in the conventional way).  Computations
00359 C              with banded matrices cost less time and storage than with
00360 C              full matrices if  2*ML+MU .lt. NEQ.  If you tell the
00361 C              code that the matrix of partial derivatives has a banded
00362 C              structure and you want to provide subroutine JAC to
00363 C              compute the partial derivatives, then you must be careful
00364 C              to store the elements of the matrix in the special form
00365 C              indicated in the description of JAC.
00366 C
00367 C          **** Do you want to solve the problem using a full (dense)
00368 C               matrix (and not a special banded structure) ...
00369 C                yes - set INFO(6) = 0
00370 C                 no - set INFO(6) = 1
00371 C                       and provide the lower (ML) and upper (MU)
00372 C                       bandwidths by setting
00373 C                       IWORK(1)=ML
00374 C                       IWORK(2)=MU ****
00375 C
00376 C       INFO(7) - You can specify a maximum (absolute value of)
00377 C              stepsize, so that the code will avoid passing over very
00378 C              large regions.
00379 C
00380 C          ****  Do you want the code to decide on its own the maximum
00381 C                stepsize ...
00382 C                 yes - set INFO(7) = 0
00383 C                  no - set INFO(7) = 1
00384 C                       and define HMAX by setting
00385 C                       RWORK(2) = HMAX ****
00386 C
00387 C       INFO(8) -  Differential/algebraic problems may occasionally
00388 C              suffer from severe scaling difficulties on the first
00389 C              step.  If you know a great deal about the scaling of 
00390 C              your problem, you can help to alleviate this problem 
00391 C              by specifying an initial stepsize H0.
00392 C
00393 C          ****  Do you want the code to define its own initial
00394 C                stepsize ...
00395 C                 yes - set INFO(8) = 0
00396 C                  no - set INFO(8) = 1
00397 C                       and define H0 by setting
00398 C                       RWORK(3) = H0 ****
00399 C
00400 C       INFO(9) -  If storage is a severe problem, you can save some
00401 C              storage by restricting the maximum method order MAXORD.
00402 C              The default value is 5.  For each order decrease below 5,
00403 C              the code requires NEQ fewer locations, but it is likely 
00404 C              to be slower.  In any case, you must have 
00405 C              1 .le. MAXORD .le. 5.
00406 C          ****  Do you want the maximum order to default to 5 ...
00407 C                 yes - set INFO(9) = 0
00408 C                  no - set INFO(9) = 1
00409 C                       and define MAXORD by setting
00410 C                       IWORK(3) = MAXORD ****
00411 C
00412 C       INFO(10) - If you know that certain components of the
00413 C              solutions to your equations are always nonnegative
00414 C              (or nonpositive), it may help to set this
00415 C              parameter.  There are three options that are
00416 C              available:
00417 C              1.  To have constraint checking only in the initial
00418 C                  condition calculation.
00419 C              2.  To enforce nonnegativity in Y during the integration.
00420 C              3.  To enforce both options 1 and 2.
00421 C
00422 C              When selecting option 2 or 3, it is probably best to try the
00423 C              code without using this option first, and only use
00424 C              this option if that does not work very well.
00425 C
00426 C          ****  Do you want the code to solve the problem without
00427 C                invoking any special inequality constraints ...
00428 C                 yes - set INFO(10) = 0
00429 C                  no - set INFO(10) = 1 to have option 1 enforced 
00430 C                  no - set INFO(10) = 2 to have option 2 enforced
00431 C                  no - set INFO(10) = 3 to have option 3 enforced ****
00432 C
00433 C                  If you have specified INFO(10) = 1 or 3, then you
00434 C                  will also need to identify how each component of Y
00435 C                  in the initial condition calculation is constrained.
00436 C                  You must set:
00437 C                  IWORK(40+I) = +1 if Y(I) must be .GE. 0,
00438 C                  IWORK(40+I) = +2 if Y(I) must be .GT. 0,
00439 C                  IWORK(40+I) = -1 if Y(I) must be .LE. 0, while
00440 C                  IWORK(40+I) = -2 if Y(I) must be .LT. 0, while
00441 C                  IWORK(40+I) =  0 if Y(I) is not constrained.
00442 C
00443 C       INFO(11) - DDASPK normally requires the initial T, Y, and
00444 C              YPRIME to be consistent.  That is, you must have
00445 C              G(T,Y,YPRIME) = 0 at the initial T.  If you do not know
00446 C              the initial conditions precisely, in some cases
00447 C              DDASPK may be able to compute it.
00448 C
00449 C              Denoting the differential variables in Y by Y_d
00450 C              and the algebraic variables by Y_a, DDASPK can solve
00451 C              one of two initialization problems:
00452 C              1.  Given Y_d, calculate Y_a and Y'_d, or
00453 C              2.  Given Y', calculate Y.
00454 C              In either case, initial values for the given
00455 C              components are input, and initial guesses for
00456 C              the unknown components must also be provided as input.
00457 C
00458 C          ****  Are the initial T, Y, YPRIME consistent ...
00459 C
00460 C                 yes - set INFO(11) = 0
00461 C                  no - set INFO(11) = 1 to calculate option 1 above,
00462 C                    or set INFO(11) = 2 to calculate option 2 ****
00463 C
00464 C                  If you have specified INFO(11) = 1, then you
00465 C                  will also need to identify  which are the
00466 C                  differential and which are the algebraic
00467 C                  components (algebraic components are components
00468 C                  whose derivatives do not appear explicitly
00469 C                  in the function G(T,Y,YPRIME)).  You must set:
00470 C                  IWORK(LID+I) = +1 if Y(I) is a differential variable
00471 C                  IWORK(LID+I) = -1 if Y(I) is an algebraic variable,
00472 C                  where LID = 40 if INFO(10) = 0 or 2 and LID = 40+NEQ
00473 C                  if INFO(10) = 1 or 3.
00474 C
00475 C       INFO(12) - Except for the addition of the RES argument CJ,
00476 C              DDASPK by default is downward-compatible with DDASSL,
00477 C              which uses only direct (dense or band) methods to solve 
00478 C              the linear systems involved.  You must set INFO(12) to
00479 C              indicate whether you want the direct methods or the
00480 C              Krylov iterative method.
00481 C          ****   Do you want DDASPK to use standard direct methods
00482 C                 (dense or band) or the Krylov (iterative) method ...
00483 C                   direct methods - set INFO(12) = 0.
00484 C                   Krylov method  - set INFO(12) = 1,
00485 C                       and check the settings of INFO(13) and INFO(15).
00486 C
00487 C       INFO(13) - used when INFO(12) = 1 (Krylov methods).  
00488 C              DDASPK uses scalars MAXL, KMP, NRMAX, and EPLI for the
00489 C              iterative solution of linear systems.  INFO(13) allows 
00490 C              you to override the default values of these parameters.  
00491 C              These parameters and their defaults are as follows:
00492 C              MAXL = maximum number of iterations in the SPIGMR 
00493 C                 algorithm (MAXL .le. NEQ).  The default is 
00494 C                 MAXL = MIN(5,NEQ).
00495 C              KMP = number of vectors on which orthogonalization is 
00496 C                 done in the SPIGMR algorithm.  The default is 
00497 C                 KMP = MAXL, which corresponds to complete GMRES 
00498 C                 iteration, as opposed to the incomplete form.  
00499 C              NRMAX = maximum number of restarts of the SPIGMR 
00500 C                 algorithm per nonlinear iteration.  The default is
00501 C                 NRMAX = 5.
00502 C              EPLI = convergence test constant in SPIGMR algorithm.
00503 C                 The default is EPLI = 0.05.
00504 C              Note that the length of RWORK depends on both MAXL 
00505 C              and KMP.  See the definition of LRW below.
00506 C          ****   Are MAXL, KMP, and EPLI to be given their
00507 C                 default values ...
00508 C                  yes - set INFO(13) = 0
00509 C                   no - set INFO(13) = 1,
00510 C                        and set all of the following:
00511 C                        IWORK(24) = MAXL (1 .le. MAXL .le. NEQ)
00512 C                        IWORK(25) = KMP  (1 .le. KMP .le. MAXL)
00513 C                        IWORK(26) = NRMAX  (NRMAX .ge. 0)
00514 C                        RWORK(10) = EPLI (0 .lt. EPLI .lt. 1.0) ****
00515 C
00516 C        INFO(14) - used with INFO(11) > 0 (initial condition 
00517 C               calculation is requested).  In this case, you may
00518 C               request control to be returned to the calling program
00519 C               immediately after the initial condition calculation,
00520 C               before proceeding to the integration of the system
00521 C               (e.g. to examine the computed Y and YPRIME).
00522 C               If this is done, and if the initialization succeeded
00523 C               (IDID = 4), you should reset INFO(11) to 0 for the
00524 C               next call, to prevent the solver from repeating the 
00525 C               initialization (and to avoid an infinite loop). 
00526 C          ****   Do you want to proceed to the integration after
00527 C                 the initial condition calculation is done ...
00528 C                 yes - set INFO(14) = 0
00529 C                  no - set INFO(14) = 1                        ****
00530 C
00531 C        INFO(15) - used when INFO(12) = 1 (Krylov methods).
00532 C               When using preconditioning in the Krylov method,
00533 C               you must supply a subroutine, PSOL, which solves the
00534 C               associated linear systems using P.
00535 C               The usage of DDASPK is simpler if PSOL can carry out
00536 C               the solution without any prior calculation of data.
00537 C               However, if some partial derivative data is to be
00538 C               calculated in advance and used repeatedly in PSOL,
00539 C               then you must supply a JAC routine to do this,
00540 C               and set INFO(15) to indicate that JAC is to be called
00541 C               for this purpose.  For example, P might be an
00542 C               approximation to a part of the matrix A which can be
00543 C               calculated and LU-factored for repeated solutions of
00544 C               the preconditioner system.  The arrays WP and IWP
00545 C               (described under JAC and PSOL) can be used to
00546 C               communicate data between JAC and PSOL.
00547 C          ****   Does PSOL operate with no prior preparation ...
00548 C                 yes - set INFO(15) = 0 (no JAC routine)
00549 C                  no - set INFO(15) = 1
00550 C                       and supply a JAC routine to evaluate and
00551 C                       preprocess any required Jacobian data.  ****
00552 C
00553 C         INFO(16) - option to exclude algebraic variables from
00554 C               the error test.  
00555 C          ****   Do you wish to control errors locally on
00556 C                 all the variables...
00557 C                 yes - set INFO(16) = 0
00558 C                  no - set INFO(16) = 1
00559 C                       If you have specified INFO(16) = 1, then you
00560 C                       will also need to identify  which are the
00561 C                       differential and which are the algebraic
00562 C                       components (algebraic components are components
00563 C                       whose derivatives do not appear explicitly
00564 C                       in the function G(T,Y,YPRIME)).  You must set:
00565 C                       IWORK(LID+I) = +1 if Y(I) is a differential 
00566 C                                      variable, and
00567 C                       IWORK(LID+I) = -1 if Y(I) is an algebraic
00568 C                                      variable,
00569 C                       where LID = 40 if INFO(10) = 0 or 2 and 
00570 C                       LID = 40 + NEQ if INFO(10) = 1 or 3.
00571 C
00572 C       INFO(17) - used when INFO(11) > 0 (DDASPK is to do an 
00573 C              initial condition calculation).
00574 C              DDASPK uses several heuristic control quantities in the
00575 C              initial condition calculation.  They have default values,
00576 C              but can  also be set by the user using INFO(17).
00577 C              These parameters and their defaults are as follows:
00578 C              MXNIT  = maximum number of Newton iterations
00579 C                 per Jacobian or preconditioner evaluation.
00580 C                 The default is:
00581 C                 MXNIT =  5 in the direct case (INFO(12) = 0), and
00582 C                 MXNIT = 15 in the Krylov case (INFO(12) = 1).
00583 C              MXNJ   = maximum number of Jacobian or preconditioner
00584 C                 evaluations.  The default is:
00585 C                 MXNJ = 6 in the direct case (INFO(12) = 0), and
00586 C                 MXNJ = 2 in the Krylov case (INFO(12) = 1).
00587 C              MXNH   = maximum number of values of the artificial
00588 C                 stepsize parameter H to be tried if INFO(11) = 1.
00589 C                 The default is MXNH = 5.
00590 C                 NOTE: the maximum number of Newton iterations
00591 C                 allowed in all is MXNIT*MXNJ*MXNH if INFO(11) = 1,
00592 C                 and MXNIT*MXNJ if INFO(11) = 2.
00593 C              LSOFF  = flag to turn off the linesearch algorithm
00594 C                 (LSOFF = 0 means linesearch is on, LSOFF = 1 means
00595 C                 it is turned off).  The default is LSOFF = 0.
00596 C              STPTOL = minimum scaled step in linesearch algorithm.
00597 C                 The default is STPTOL = (unit roundoff)**(2/3).
00598 C              EPINIT = swing factor in the Newton iteration convergence
00599 C                 test.  The test is applied to the residual vector,
00600 C                 premultiplied by the approximate Jacobian (in the
00601 C                 direct case) or the preconditioner (in the Krylov
00602 C                 case).  For convergence, the weighted RMS norm of
00603 C                 this vector (scaled by the error weights) must be
00604 C                 less than EPINIT*EPCON, where EPCON = .33 is the
00605 C                 analogous test constant used in the time steps.
00606 C                 The default is EPINIT = .01.
00607 C          ****   Are the initial condition heuristic controls to be 
00608 C                 given their default values...
00609 C                  yes - set INFO(17) = 0
00610 C                   no - set INFO(17) = 1,
00611 C                        and set all of the following:
00612 C                        IWORK(32) = MXNIT (.GT. 0)
00613 C                        IWORK(33) = MXNJ (.GT. 0)
00614 C                        IWORK(34) = MXNH (.GT. 0)
00615 C                        IWORK(35) = LSOFF ( = 0 or 1)
00616 C                        RWORK(14) = STPTOL (.GT. 0.0)
00617 C                        RWORK(15) = EPINIT (.GT. 0.0)  ****
00618 C
00619 C         INFO(18) - option to get extra printing in initial condition 
00620 C                calculation.
00621 C          ****   Do you wish to have extra printing...
00622 C                 no  - set INFO(18) = 0
00623 C                 yes - set INFO(18) = 1 for minimal printing, or
00624 C                       set INFO(18) = 2 for full printing.
00625 C                       If you have specified INFO(18) .ge. 1, data
00626 C                       will be printed with the error handler routines.
00627 C                       To print to a non-default unit number L, include
00628 C                       the line  CALL XSETUN(L)  in your program.  ****
00629 C
00630 C   RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL)
00631 C               error tolerances to tell the code how accurately you
00632 C               want the solution to be computed.  They must be defined
00633 C               as variables because the code may change them.
00634 C               you have two choices --
00635 C                     Both RTOL and ATOL are scalars (INFO(2) = 0), or
00636 C                     both RTOL and ATOL are vectors (INFO(2) = 1).
00637 C               In either case all components must be non-negative.
00638 C
00639 C               The tolerances are used by the code in a local error
00640 C               test at each step which requires roughly that
00641 C                        abs(local error in Y(i)) .le. EWT(i) ,
00642 C               where EWT(i) = RTOL*abs(Y(i)) + ATOL is an error weight 
00643 C               quantity, for each vector component.
00644 C               (More specifically, a root-mean-square norm is used to
00645 C               measure the size of vectors, and the error test uses the
00646 C               magnitude of the solution at the beginning of the step.)
00647 C
00648 C               The true (global) error is the difference between the
00649 C               true solution of the initial value problem and the
00650 C               computed approximation.  Practically all present day
00651 C               codes, including this one, control the local error at
00652 C               each step and do not even attempt to control the global
00653 C               error directly.
00654 C
00655 C               Usually, but not always, the true accuracy of
00656 C               the computed Y is comparable to the error tolerances.
00657 C               This code will usually, but not always, deliver a more
00658 C               accurate solution if you reduce the tolerances and
00659 C               integrate again.  By comparing two such solutions you 
00660 C               can get a fairly reliable idea of the true error in the
00661 C               solution at the larger tolerances.
00662 C
00663 C               Setting ATOL = 0. results in a pure relative error test
00664 C               on that component.  Setting RTOL = 0. results in a pure
00665 C               absolute error test on that component.  A mixed test
00666 C               with non-zero RTOL and ATOL corresponds roughly to a
00667 C               relative error test when the solution component is
00668 C               much bigger than ATOL and to an absolute error test
00669 C               when the solution component is smaller than the
00670 C               threshold ATOL.
00671 C
00672 C               The code will not attempt to compute a solution at an
00673 C               accuracy unreasonable for the machine being used.  It
00674 C               will advise you if you ask for too much accuracy and
00675 C               inform you as to the maximum accuracy it believes
00676 C               possible.
00677 C
00678 C  RWORK(*) -- a real work array, which should be dimensioned in your
00679 C               calling program with a length equal to the value of
00680 C               LRW (or greater).
00681 C
00682 C  LRW -- Set it to the declared length of the RWORK array.  The
00683 C               minimum length depends on the options you have selected,
00684 C               given by a base value plus additional storage as described
00685 C               below.
00686 C
00687 C               If INFO(12) = 0 (standard direct method), the base value is
00688 C               base = 50 + max(MAXORD+4,7)*NEQ.
00689 C               The default value is MAXORD = 5 (see INFO(9)).  With the
00690 C               default MAXORD, base = 50 + 9*NEQ.
00691 C               Additional storage must be added to the base value for
00692 C               any or all of the following options:
00693 C                 if INFO(6) = 0 (dense matrix), add NEQ**2
00694 C                 if INFO(6) = 1 (banded matrix), then
00695 C                    if INFO(5) = 0, add (2*ML+MU+1)*NEQ + 2*(NEQ/(ML+MU+1)+1),
00696 C                    if INFO(5) = 1, add (2*ML+MU+1)*NEQ,
00697 C                 if INFO(16) = 1, add NEQ.
00698 C
00699 C              If INFO(12) = 1 (Krylov method), the base value is
00700 C              base = 50 + (MAXORD+5)*NEQ + (MAXL+3+MIN0(1,MAXL-KMP))*NEQ +
00701 C                      + (MAXL+3)*MAXL + 1 + LENWP.
00702 C              See PSOL for description of LENWP.  The default values are:
00703 C              MAXORD = 5 (see INFO(9)), MAXL = min(5,NEQ) and KMP = MAXL 
00704 C              (see INFO(13)).
00705 C              With the default values for MAXORD, MAXL and KMP,
00706 C              base = 91 + 18*NEQ + LENWP.
00707 C              Additional storage must be added to the base value for
00708 C              any or all of the following options:
00709 C                if INFO(16) = 1, add NEQ.
00710 C
00711 C
00712 C  IWORK(*) -- an integer work array, which should be dimensioned in
00713 C              your calling program with a length equal to the value
00714 C              of LIW (or greater).
00715 C
00716 C  LIW -- Set it to the declared length of the IWORK array.  The
00717 C             minimum length depends on the options you have selected,
00718 C             given by a base value plus additional storage as described
00719 C             below.
00720 C
00721 C             If INFO(12) = 0 (standard direct method), the base value is
00722 C             base = 40 + NEQ.
00723 C             IF INFO(10) = 1 or 3, add NEQ to the base value.
00724 C             If INFO(11) = 1 or INFO(16) =1, add NEQ to the base value.
00725 C
00726 C             If INFO(12) = 1 (Krylov method), the base value is
00727 C             base = 40 + LENIWP.
00728 C             See PSOL for description of LENIWP.
00729 C             IF INFO(10) = 1 or 3, add NEQ to the base value.
00730 C             If INFO(11) = 1 or INFO(16) = 1, add NEQ to the base value.
00731 C
00732 C
00733 C  RPAR, IPAR -- These are arrays of double precision and integer type,
00734 C             respectively, which are available for you to use
00735 C             for communication between your program that calls
00736 C             DDASPK and the RES subroutine (and the JAC and PSOL
00737 C             subroutines).  They are not altered by DDASPK.
00738 C             If you do not need RPAR or IPAR, ignore these
00739 C             parameters by treating them as dummy arguments.
00740 C             If you do choose to use them, dimension them in
00741 C             your calling program and in RES (and in JAC and PSOL)
00742 C             as arrays of appropriate length.
00743 C
00744 C  JAC -- This is the name of a routine that you may supply
00745 C         (optionally) that relates to the Jacobian matrix of the
00746 C         nonlinear system that the code must solve at each T step.
00747 C         The role of JAC (and its call sequence) depends on whether
00748 C         a direct (INFO(12) = 0) or Krylov (INFO(12) = 1) method 
00749 C         is selected.
00750 C
00751 C         **** INFO(12) = 0 (direct methods):
00752 C           If you are letting the code generate partial derivatives
00753 C           numerically (INFO(5) = 0), then JAC can be absent
00754 C           (or perhaps a dummy routine to satisfy the loader).
00755 C           Otherwise you must supply a JAC routine to compute
00756 C           the matrix A = dG/dY + CJ*dG/dYPRIME.  It must have
00757 C           the form
00758 C
00759 C           SUBROUTINE JAC (T, Y, YPRIME, PD, CJ, RPAR, IPAR)
00760 C
00761 C           The JAC routine must dimension Y, YPRIME, and PD (and RPAR
00762 C           and IPAR if used).  CJ is a scalar which is input to JAC.
00763 C           For the given values of T, Y, and YPRIME, the JAC routine
00764 C           must evaluate the nonzero elements of the matrix A, and 
00765 C           store these values in the array PD.  The elements of PD are 
00766 C           set to zero before each call to JAC, so that only nonzero
00767 C           elements need to be defined.
00768 C           The way you store the elements into the PD array depends
00769 C           on the structure of the matrix indicated by INFO(6).
00770 C           *** INFO(6) = 0 (full or dense matrix) ***
00771 C               Give PD a first dimension of NEQ.  When you evaluate the
00772 C               nonzero partial derivatives of equation i (i.e. of G(i))
00773 C               with respect to component j (of Y and YPRIME), you must
00774 C               store the element in PD according to
00775 C                  PD(i,j) = dG(i)/dY(j) + CJ*dG(i)/dYPRIME(j).
00776 C           *** INFO(6) = 1 (banded matrix with half-bandwidths ML, MU
00777 C                            as described under INFO(6)) ***
00778 C               Give PD a first dimension of 2*ML+MU+1.  When you 
00779 C               evaluate the nonzero partial derivatives of equation i 
00780 C               (i.e. of G(i)) with respect to component j (of Y and 
00781 C               YPRIME), you must store the element in PD according to 
00782 C                  IROW = i - j + ML + MU + 1
00783 C                  PD(IROW,j) = dG(i)/dY(j) + CJ*dG(i)/dYPRIME(j).
00784 C
00785 C          **** INFO(12) = 1 (Krylov method):
00786 C            If you are not calculating Jacobian data in advance for use
00787 C            in PSOL (INFO(15) = 0), JAC can be absent (or perhaps a
00788 C            dummy routine to satisfy the loader).  Otherwise, you may
00789 C            supply a JAC routine to compute and preprocess any parts of
00790 C            of the Jacobian matrix  A = dG/dY + CJ*dG/dYPRIME that are
00791 C            involved in the preconditioner matrix P.
00792 C            It is to have the form
00793 C
00794 C            SUBROUTINE JAC (RES, IRES, NEQ, T, Y, YPRIME, REWT, SAVR,
00795 C                            WK, H, CJ, WP, IWP, IER, RPAR, IPAR)
00796 C
00797 C           The JAC routine must dimension Y, YPRIME, REWT, SAVR, WK,
00798 C           and (if used) WP, IWP, RPAR, and IPAR.
00799 C           The Y, YPRIME, and SAVR arrays contain the current values
00800 C           of Y, YPRIME, and the residual G, respectively.  
00801 C           The array WK is work space of length NEQ.  
00802 C           H is the step size.  CJ is a scalar, input to JAC, that is
00803 C           normally proportional to 1/H.  REWT is an array of 
00804 C           reciprocal error weights, 1/EWT(i), where EWT(i) is
00805 C           RTOL*abs(Y(i)) + ATOL (unless you supplied routine DDAWTS
00806 C           instead), for use in JAC if needed.  For example, if JAC
00807 C           computes difference quotient approximations to partial
00808 C           derivatives, the REWT array may be useful in setting the
00809 C           increments used.  The JAC routine should do any
00810 C           factorization operations called for, in preparation for
00811 C           solving linear systems in PSOL.  The matrix P should
00812 C           be an approximation to the Jacobian,
00813 C           A = dG/dY + CJ*dG/dYPRIME.
00814 C
00815 C           WP and IWP are real and integer work arrays which you may
00816 C           use for communication between your JAC routine and your
00817 C           PSOL routine.  These may be used to store elements of the 
00818 C           preconditioner P, or related matrix data (such as factored
00819 C           forms).  They are not altered by DDASPK.
00820 C           If you do not need WP or IWP, ignore these parameters by
00821 C           treating them as dummy arguments.  If you do use them,
00822 C           dimension them appropriately in your JAC and PSOL routines.
00823 C           See the PSOL description for instructions on setting 
00824 C           the lengths of WP and IWP.
00825 C
00826 C           On return, JAC should set the error flag IER as follows..
00827 C             IER = 0    if JAC was successful,
00828 C             IER .ne. 0 if JAC was unsuccessful (e.g. if Y or YPRIME
00829 C                        was illegal, or a singular matrix is found).
00830 C           (If IER .ne. 0, a smaller stepsize will be tried.)
00831 C           IER = 0 on entry to JAC, so need be reset only on a failure.
00832 C           If RES is used within JAC, then a nonzero value of IRES will
00833 C           override any nonzero value of IER (see the RES description).
00834 C
00835 C         Regardless of the method type, subroutine JAC must not
00836 C         alter T, Y(*), YPRIME(*), H, CJ, or REWT(*).
00837 C         You must declare the name JAC in an EXTERNAL statement in
00838 C         your program that calls DDASPK.
00839 C
00840 C PSOL --  This is the name of a routine you must supply if you have
00841 C         selected a Krylov method (INFO(12) = 1) with preconditioning.
00842 C         In the direct case (INFO(12) = 0), PSOL can be absent 
00843 C         (a dummy routine may have to be supplied to satisfy the 
00844 C         loader).  Otherwise, you must provide a PSOL routine to 
00845 C         solve linear systems arising from preconditioning.
00846 C         When supplied with INFO(12) = 1, the PSOL routine is to 
00847 C         have the form
00848 C
00849 C         SUBROUTINE PSOL (NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT,
00850 C                          WP, IWP, B, EPLIN, IER, RPAR, IPAR)
00851 C
00852 C         The PSOL routine must solve linear systems of the form 
00853 C         P*x = b where P is the left preconditioner matrix.
00854 C
00855 C         The right-hand side vector b is in the B array on input, and
00856 C         PSOL must return the solution vector x in B.
00857 C         The Y, YPRIME, and SAVR arrays contain the current values
00858 C         of Y, YPRIME, and the residual G, respectively.  
00859 C
00860 C         Work space required by JAC and/or PSOL, and space for data to
00861 C         be communicated from JAC to PSOL is made available in the form
00862 C         of arrays WP and IWP, which are parts of the RWORK and IWORK
00863 C         arrays, respectively.  The lengths of these real and integer
00864 C         work spaces WP and IWP must be supplied in LENWP and LENIWP,
00865 C         respectively, as follows..
00866 C           IWORK(27) = LENWP = length of real work space WP
00867 C           IWORK(28) = LENIWP = length of integer work space IWP.
00868 C
00869 C         WK is a work array of length NEQ for use by PSOL.
00870 C         CJ is a scalar, input to PSOL, that is normally proportional
00871 C         to 1/H (H = stepsize).  If the old value of CJ
00872 C         (at the time of the last JAC call) is needed, it must have
00873 C         been saved by JAC in WP.
00874 C
00875 C         WGHT is an array of weights, to be used if PSOL uses an
00876 C         iterative method and performs a convergence test.  (In terms
00877 C         of the argument REWT to JAC, WGHT is REWT/sqrt(NEQ).)
00878 C         If PSOL uses an iterative method, it should use EPLIN
00879 C         (a heuristic parameter) as the bound on the weighted norm of
00880 C         the residual for the computed solution.  Specifically, the
00881 C         residual vector R should satisfy
00882 C              SQRT (SUM ( (R(i)*WGHT(i))**2 ) ) .le. EPLIN
00883 C
00884 C         PSOL must not alter NEQ, T, Y, YPRIME, SAVR, CJ, WGHT, EPLIN.
00885 C
00886 C         On return, PSOL should set the error flag IER as follows..
00887 C           IER = 0 if PSOL was successful,
00888 C           IER .lt. 0 if an unrecoverable error occurred, meaning
00889 C                 control will be passed to the calling routine,
00890 C           IER .gt. 0 if a recoverable error occurred, meaning that
00891 C                 the step will be retried with the same step size
00892 C                 but with a call to JAC to update necessary data,
00893 C                 unless the Jacobian data is current, in which case
00894 C                 the step will be retried with a smaller step size.
00895 C           IER = 0 on entry to PSOL so need be reset only on a failure.
00896 C
00897 C         You must declare the name PSOL in an EXTERNAL statement in
00898 C         your program that calls DDASPK.
00899 C
00900 C
00901 C  OPTIONALLY REPLACEABLE SUBROUTINE:
00902 C
00903 C  DDASPK uses a weighted root-mean-square norm to measure the 
00904 C  size of various error vectors.  The weights used in this norm
00905 C  are set in the following subroutine:
00906 C
00907 C    SUBROUTINE DDAWTS (NEQ, IWT, RTOL, ATOL, Y, EWT, RPAR, IPAR)
00908 C    DIMENSION RTOL(*), ATOL(*), Y(*), EWT(*), RPAR(*), IPAR(*)
00909 C
00910 C  A DDAWTS routine has been included with DDASPK which sets the
00911 C  weights according to
00912 C    EWT(I) = RTOL*ABS(Y(I)) + ATOL
00913 C  in the case of scalar tolerances (IWT = 0) or
00914 C    EWT(I) = RTOL(I)*ABS(Y(I)) + ATOL(I)
00915 C  in the case of array tolerances (IWT = 1).  (IWT is INFO(2).)
00916 C  In some special cases, it may be appropriate for you to define
00917 C  your own error weights by writing a subroutine DDAWTS to be 
00918 C  called instead of the version supplied.  However, this should 
00919 C  be attempted only after careful thought and consideration. 
00920 C  If you supply this routine, you may use the tolerances and Y 
00921 C  as appropriate, but do not overwrite these variables.  You
00922 C  may also use RPAR and IPAR to communicate data as appropriate.
00923 C  ***Note: Aside from the values of the weights, the choice of 
00924 C  norm used in DDASPK (weighted root-mean-square) is not subject
00925 C  to replacement by the user.  In this respect, DDASPK is not
00926 C  downward-compatible with the original DDASSL solver (in which
00927 C  the norm routine was optionally user-replaceable).
00928 C
00929 C
00930 C------OUTPUT - AFTER ANY RETURN FROM DDASPK----------------------------
00931 C
00932 C  The principal aim of the code is to return a computed solution at
00933 C  T = TOUT, although it is also possible to obtain intermediate
00934 C  results along the way.  To find out whether the code achieved its
00935 C  goal or if the integration process was interrupted before the task
00936 C  was completed, you must check the IDID parameter.
00937 C
00938 C
00939 C   T -- The output value of T is the point to which the solution
00940 C        was successfully advanced.
00941 C
00942 C   Y(*) -- contains the computed solution approximation at T.
00943 C
00944 C   YPRIME(*) -- contains the computed derivative approximation at T.
00945 C
00946 C   IDID -- reports what the code did, described as follows:
00947 C
00948 C                     *** TASK COMPLETED ***
00949 C                Reported by positive values of IDID
00950 C
00951 C           IDID = 1 -- a step was successfully taken in the
00952 C                   intermediate-output mode.  The code has not
00953 C                   yet reached TOUT.
00954 C
00955 C           IDID = 2 -- the integration to TSTOP was successfully
00956 C                   completed (T = TSTOP) by stepping exactly to TSTOP.
00957 C
00958 C           IDID = 3 -- the integration to TOUT was successfully
00959 C                   completed (T = TOUT) by stepping past TOUT.
00960 C                   Y(*) and YPRIME(*) are obtained by interpolation.
00961 C
00962 C           IDID = 4 -- the initial condition calculation, with
00963 C                   INFO(11) > 0, was successful, and INFO(14) = 1.
00964 C                   No integration steps were taken, and the solution
00965 C                   is not considered to have been started.
00966 C
00967 C                    *** TASK INTERRUPTED ***
00968 C                Reported by negative values of IDID
00969 C
00970 C           IDID = -1 -- a large amount of work has been expended
00971 C                     (about 500 steps).
00972 C
00973 C           IDID = -2 -- the error tolerances are too stringent.
00974 C
00975 C           IDID = -3 -- the local error test cannot be satisfied
00976 C                     because you specified a zero component in ATOL
00977 C                     and the corresponding computed solution component
00978 C                     is zero.  Thus, a pure relative error test is
00979 C                     impossible for this component.
00980 C
00981 C           IDID = -5 -- there were repeated failures in the evaluation
00982 C                     or processing of the preconditioner (in JAC).
00983 C
00984 C           IDID = -6 -- DDASPK had repeated error test failures on the
00985 C                     last attempted step.
00986 C
00987 C           IDID = -7 -- the nonlinear system solver in the time integration
00988 C                     could not converge.
00989 C
00990 C           IDID = -8 -- the matrix of partial derivatives appears
00991 C                     to be singular (direct method).
00992 C
00993 C           IDID = -9 -- the nonlinear system solver in the time integration
00994 C                     failed to achieve convergence, and there were repeated 
00995 C                     error test failures in this step.
00996 C
00997 C           IDID =-10 -- the nonlinear system solver in the time integration 
00998 C                     failed to achieve convergence because IRES was equal 
00999 C                     to -1.
01000 C
01001 C           IDID =-11 -- IRES = -2 was encountered and control is
01002 C                     being returned to the calling program.
01003 C
01004 C           IDID =-12 -- DDASPK failed to compute the initial Y, YPRIME.
01005 C
01006 C           IDID =-13 -- unrecoverable error encountered inside user's
01007 C                     PSOL routine, and control is being returned to
01008 C                     the calling program.
01009 C
01010 C           IDID =-14 -- the Krylov linear system solver could not 
01011 C                     achieve convergence.
01012 C
01013 C           IDID =-15,..,-32 -- Not applicable for this code.
01014 C
01015 C                    *** TASK TERMINATED ***
01016 C                reported by the value of IDID=-33
01017 C
01018 C           IDID = -33 -- the code has encountered trouble from which
01019 C                   it cannot recover.  A message is printed
01020 C                   explaining the trouble and control is returned
01021 C                   to the calling program.  For example, this occurs
01022 C                   when invalid input is detected.
01023 C
01024 C   RTOL, ATOL -- these quantities remain unchanged except when
01025 C               IDID = -2.  In this case, the error tolerances have been
01026 C               increased by the code to values which are estimated to
01027 C               be appropriate for continuing the integration.  However,
01028 C               the reported solution at T was obtained using the input
01029 C               values of RTOL and ATOL.
01030 C
01031 C   RWORK, IWORK -- contain information which is usually of no interest
01032 C               to the user but necessary for subsequent calls. 
01033 C               However, you may be interested in the performance data
01034 C               listed below.  These quantities are accessed in RWORK 
01035 C               and IWORK but have internal mnemonic names, as follows..
01036 C
01037 C               RWORK(3)--contains H, the step size h to be attempted
01038 C                        on the next step.
01039 C
01040 C               RWORK(4)--contains TN, the current value of the
01041 C                        independent variable, i.e. the farthest point
01042 C                        integration has reached.  This will differ 
01043 C                        from T if interpolation has been performed 
01044 C                        (IDID = 3).
01045 C
01046 C               RWORK(7)--contains HOLD, the stepsize used on the last
01047 C                        successful step.  If INFO(11) = INFO(14) = 1,
01048 C                        this contains the value of H used in the
01049 C                        initial condition calculation.
01050 C
01051 C               IWORK(7)--contains K, the order of the method to be 
01052 C                        attempted on the next step.
01053 C
01054 C               IWORK(8)--contains KOLD, the order of the method used
01055 C                        on the last step.
01056 C
01057 C               IWORK(11)--contains NST, the number of steps (in T) 
01058 C                        taken so far.
01059 C
01060 C               IWORK(12)--contains NRE, the number of calls to RES 
01061 C                        so far.
01062 C
01063 C               IWORK(13)--contains NJE, the number of calls to JAC so
01064 C                        far (Jacobian or preconditioner evaluations).
01065 C
01066 C               IWORK(14)--contains NETF, the total number of error test
01067 C                        failures so far.
01068 C
01069 C               IWORK(15)--contains NCFN, the total number of nonlinear
01070 C                        convergence failures so far (includes counts
01071 C                        of singular iteration matrix or singular
01072 C                        preconditioners).
01073 C
01074 C               IWORK(16)--contains NCFL, the number of convergence
01075 C                        failures of the linear iteration so far.
01076 C
01077 C               IWORK(17)--contains LENIW, the length of IWORK actually
01078 C                        required.  This is defined on normal returns 
01079 C                        and on an illegal input return for
01080 C                        insufficient storage.
01081 C
01082 C               IWORK(18)--contains LENRW, the length of RWORK actually
01083 C                        required.  This is defined on normal returns 
01084 C                        and on an illegal input return for
01085 C                        insufficient storage.
01086 C
01087 C               IWORK(19)--contains NNI, the total number of nonlinear
01088 C                        iterations so far (each of which calls a
01089 C                        linear solver).
01090 C
01091 C               IWORK(20)--contains NLI, the total number of linear
01092 C                        (Krylov) iterations so far.
01093 C
01094 C               IWORK(21)--contains NPS, the number of PSOL calls so
01095 C                        far, for preconditioning solve operations or
01096 C                        for solutions with the user-supplied method.
01097 C
01098 C               Note: The various counters in IWORK do not include 
01099 C               counts during a call made with INFO(11) > 0 and
01100 C               INFO(14) = 1.
01101 C
01102 C
01103 C------INPUT - WHAT TO DO TO CONTINUE THE INTEGRATION  -----------------
01104 C              (CALLS AFTER THE FIRST)
01105 C
01106 C     This code is organized so that subsequent calls to continue the
01107 C     integration involve little (if any) additional effort on your
01108 C     part.  You must monitor the IDID parameter in order to determine
01109 C     what to do next.
01110 C
01111 C     Recalling that the principal task of the code is to integrate
01112 C     from T to TOUT (the interval mode), usually all you will need
01113 C     to do is specify a new TOUT upon reaching the current TOUT.
01114 C
01115 C     Do not alter any quantity not specifically permitted below.  In
01116 C     particular do not alter NEQ, T, Y(*), YPRIME(*), RWORK(*), 
01117 C     IWORK(*), or the differential equation in subroutine RES.  Any 
01118 C     such alteration constitutes a new problem and must be treated 
01119 C     as such, i.e. you must start afresh.
01120 C
01121 C     You cannot change from array to scalar error control or vice
01122 C     versa (INFO(2)), but you can change the size of the entries of
01123 C     RTOL or ATOL.  Increasing a tolerance makes the equation easier
01124 C     to integrate.  Decreasing a tolerance will make the equation
01125 C     harder to integrate and should generally be avoided.
01126 C
01127 C     You can switch from the intermediate-output mode to the
01128 C     interval mode (INFO(3)) or vice versa at any time.
01129 C
01130 C     If it has been necessary to prevent the integration from going
01131 C     past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
01132 C     code will not integrate to any TOUT beyond the currently
01133 C     specified TSTOP.  Once TSTOP has been reached, you must change
01134 C     the value of TSTOP or set INFO(4) = 0.  You may change INFO(4)
01135 C     or TSTOP at any time but you must supply the value of TSTOP in
01136 C     RWORK(1) whenever you set INFO(4) = 1.
01137 C
01138 C     Do not change INFO(5), INFO(6), INFO(12-17) or their associated
01139 C     IWORK/RWORK locations unless you are going to restart the code.
01140 C
01141 C                    *** FOLLOWING A COMPLETED TASK ***
01142 C
01143 C     If..
01144 C     IDID = 1, call the code again to continue the integration
01145 C                  another step in the direction of TOUT.
01146 C
01147 C     IDID = 2 or 3, define a new TOUT and call the code again.
01148 C                  TOUT must be different from T.  You cannot change
01149 C                  the direction of integration without restarting.
01150 C
01151 C     IDID = 4, reset INFO(11) = 0 and call the code again to begin
01152 C                  the integration.  (If you leave INFO(11) > 0 and
01153 C                  INFO(14) = 1, you may generate an infinite loop.)
01154 C                  In this situation, the next call to DASPK is 
01155 C                  considered to be the first call for the problem,
01156 C                  in that all initializations are done.
01157 C
01158 C                    *** FOLLOWING AN INTERRUPTED TASK ***
01159 C
01160 C     To show the code that you realize the task was interrupted and
01161 C     that you want to continue, you must take appropriate action and
01162 C     set INFO(1) = 1.
01163 C
01164 C     If..
01165 C     IDID = -1, the code has taken about 500 steps.  If you want to
01166 C                  continue, set INFO(1) = 1 and call the code again.
01167 C                  An additional 500 steps will be allowed.
01168 C
01169 C
01170 C     IDID = -2, the error tolerances RTOL, ATOL have been increased
01171 C                  to values the code estimates appropriate for
01172 C                  continuing.  You may want to change them yourself.
01173 C                  If you are sure you want to continue with relaxed
01174 C                  error tolerances, set INFO(1) = 1 and call the code
01175 C                  again.
01176 C
01177 C     IDID = -3, a solution component is zero and you set the
01178 C                  corresponding component of ATOL to zero.  If you
01179 C                  are sure you want to continue, you must first alter
01180 C                  the error criterion to use positive values of ATOL 
01181 C                  for those components corresponding to zero solution
01182 C                  components, then set INFO(1) = 1 and call the code
01183 C                  again.
01184 C
01185 C     IDID = -4  --- cannot occur with this code.
01186 C
01187 C     IDID = -5, your JAC routine failed with the Krylov method.  Check
01188 C                  for errors in JAC and restart the integration.
01189 C
01190 C     IDID = -6, repeated error test failures occurred on the last
01191 C                  attempted step in DDASPK.  A singularity in the
01192 C                  solution may be present.  If you are absolutely
01193 C                  certain you want to continue, you should restart
01194 C                  the integration.  (Provide initial values of Y and
01195 C                  YPRIME which are consistent.)
01196 C
01197 C     IDID = -7, repeated convergence test failures occurred on the last
01198 C                  attempted step in DDASPK.  An inaccurate or ill-
01199 C                  conditioned Jacobian or preconditioner may be the
01200 C                  problem.  If you are absolutely certain you want
01201 C                  to continue, you should restart the integration.
01202 C
01203 C
01204 C     IDID = -8, the matrix of partial derivatives is singular, with
01205 C                  the use of direct methods.  Some of your equations
01206 C                  may be redundant.  DDASPK cannot solve the problem
01207 C                  as stated.  It is possible that the redundant
01208 C                  equations could be removed, and then DDASPK could
01209 C                  solve the problem.  It is also possible that a
01210 C                  solution to your problem either does not exist
01211 C                  or is not unique.
01212 C
01213 C     IDID = -9, DDASPK had multiple convergence test failures, preceded
01214 C                  by multiple error test failures, on the last
01215 C                  attempted step.  It is possible that your problem is
01216 C                  ill-posed and cannot be solved using this code.  Or,
01217 C                  there may be a discontinuity or a singularity in the
01218 C                  solution.  If you are absolutely certain you want to
01219 C                  continue, you should restart the integration.
01220 C
01221 C     IDID = -10, DDASPK had multiple convergence test failures
01222 C                  because IRES was equal to -1.  If you are
01223 C                  absolutely certain you want to continue, you
01224 C                  should restart the integration.
01225 C
01226 C     IDID = -11, there was an unrecoverable error (IRES = -2) from RES
01227 C                  inside the nonlinear system solver.  Determine the
01228 C                  cause before trying again.
01229 C
01230 C     IDID = -12, DDASPK failed to compute the initial Y and YPRIME
01231 C                  vectors.  This could happen because the initial 
01232 C                  approximation to Y or YPRIME was not very good, or
01233 C                  because no consistent values of these vectors exist.
01234 C                  The problem could also be caused by an inaccurate or
01235 C                  singular iteration matrix, or a poor preconditioner.
01236 C
01237 C     IDID = -13, there was an unrecoverable error encountered inside 
01238 C                  your PSOL routine.  Determine the cause before 
01239 C                  trying again.
01240 C
01241 C     IDID = -14, the Krylov linear system solver failed to achieve
01242 C                  convergence.  This may be due to ill-conditioning
01243 C                  in the iteration matrix, or a singularity in the
01244 C                  preconditioner (if one is being used).
01245 C                  Another possibility is that there is a better
01246 C                  choice of Krylov parameters (see INFO(13)).
01247 C                  Possibly the failure is caused by redundant equations
01248 C                  in the system, or by inconsistent equations.
01249 C                  In that case, reformulate the system to make it
01250 C                  consistent and non-redundant.
01251 C
01252 C     IDID = -15,..,-32 --- Cannot occur with this code.
01253 C
01254 C                       *** FOLLOWING A TERMINATED TASK ***
01255 C
01256 C     If IDID = -33, you cannot continue the solution of this problem.
01257 C                  An attempt to do so will result in your run being
01258 C                  terminated.
01259 C
01260 C  ---------------------------------------------------------------------
01261 C
01262 C***REFERENCES
01263 C  1.  L. R. Petzold, A Description of DASSL: A Differential/Algebraic
01264 C      System Solver, in Scientific Computing, R. S. Stepleman et al.
01265 C      (Eds.), North-Holland, Amsterdam, 1983, pp. 65-68.
01266 C  2.  K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical 
01267 C      Solution of Initial-Value Problems in Differential-Algebraic
01268 C      Equations, Elsevier, New York, 1989.
01269 C  3.  P. N. Brown and A. C. Hindmarsh, Reduced Storage Matrix Methods
01270 C      in Stiff ODE Systems, J. Applied Mathematics and Computation,
01271 C      31 (1989), pp. 40-91.
01272 C  4.  P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Using Krylov
01273 C      Methods in the Solution of Large-Scale Differential-Algebraic
01274 C      Systems, SIAM J. Sci. Comp., 15 (1994), pp. 1467-1488.
01275 C  5.  P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Consistent
01276 C      Initial Condition Calculation for Differential-Algebraic
01277 C      Systems, LLNL Report UCRL-JC-122175, August 1995; submitted to
01278 C      SIAM J. Sci. Comp.
01279 C
01280 C***ROUTINES CALLED
01281 C
01282 C   The following are all the subordinate routines used by DDASPK.
01283 C
01284 C   DDASIC computes consistent initial conditions.
01285 C   DYYPNW updates Y and YPRIME in linesearch for initial condition
01286 C          calculation.
01287 C   DDSTP  carries out one step of the integration.
01288 C   DCNSTR/DCNST0 check the current solution for constraint violations.
01289 C   DDAWTS sets error weight quantities.
01290 C   DINVWT tests and inverts the error weights.
01291 C   DDATRP performs interpolation to get an output solution.
01292 C   DDWNRM computes the weighted root-mean-square norm of a vector.
01293 C   D1MACH provides the unit roundoff of the computer.
01294 C   XERRWD/XSETF/XSETUN/IXSAV is a package to handle error messages. 
01295 C   DDASID nonlinear equation driver to initialize Y and YPRIME using
01296 C          direct linear system solver methods.  Interfaces to Newton
01297 C          solver (direct case).
01298 C   DNSID  solves the nonlinear system for unknown initial values by
01299 C          modified Newton iteration and direct linear system methods.
01300 C   DLINSD carries out linesearch algorithm for initial condition
01301 C          calculation (direct case).
01302 C   DFNRMD calculates weighted norm of preconditioned residual in
01303 C          initial condition calculation (direct case).
01304 C   DNEDD  nonlinear equation driver for direct linear system solver
01305 C          methods.  Interfaces to Newton solver (direct case).
01306 C   DMATD  assembles the iteration matrix (direct case).
01307 C   DNSD   solves the associated nonlinear system by modified
01308 C          Newton iteration and direct linear system methods.
01309 C   DSLVD  interfaces to linear system solver (direct case).
01310 C   DDASIK nonlinear equation driver to initialize Y and YPRIME using
01311 C          Krylov iterative linear system methods.  Interfaces to
01312 C          Newton solver (Krylov case).
01313 C   DNSIK  solves the nonlinear system for unknown initial values by
01314 C          Newton iteration and Krylov iterative linear system methods.
01315 C   DLINSK carries out linesearch algorithm for initial condition
01316 C          calculation (Krylov case).
01317 C   DFNRMK calculates weighted norm of preconditioned residual in
01318 C          initial condition calculation (Krylov case).
01319 C   DNEDK  nonlinear equation driver for iterative linear system solver
01320 C          methods.  Interfaces to Newton solver (Krylov case).
01321 C   DNSK   solves the associated nonlinear system by Inexact Newton
01322 C          iteration and (linear) Krylov iteration.
01323 C   DSLVK  interfaces to linear system solver (Krylov case).
01324 C   DSPIGM solves a linear system by SPIGMR algorithm.
01325 C   DATV   computes matrix-vector product in Krylov algorithm.
01326 C   DORTH  performs orthogonalization of Krylov basis vectors.
01327 C   DHEQR  performs QR factorization of Hessenberg matrix.
01328 C   DHELS  finds least-squares solution of Hessenberg linear system.
01329 C   DGETRF, DGETRS, DGBTRF, DGBTRS are LAPACK routines for solving 
01330 C          linear systems (dense or band direct methods).
01331 C   DAXPY, DCOPY, DDOT, DNRM2, DSCAL are Basic Linear Algebra (BLAS)
01332 C          routines.
01333 C
01334 C The routines called directly by DDASPK are:
01335 C   DCNST0, DDAWTS, DINVWT, D1MACH, DDWNRM, DDASIC, DDATRP, DDSTP,
01336 C   XERRWD
01337 C
01338 C***END PROLOGUE DDASPK
01339 C
01340 C
01341       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
01342       LOGICAL DONE, LAVL, LCFN, LCFL, LWARN
01343       DIMENSION Y(*),YPRIME(*)
01344       DIMENSION INFO(20)
01345       DIMENSION RWORK(LRW),IWORK(LIW)
01346       DIMENSION RTOL(*),ATOL(*)
01347       DIMENSION RPAR(*),IPAR(*)
01348       CHARACTER MSG*80
01349       EXTERNAL  RES, JAC, PSOL, DDASID, DDASIK, DNEDD, DNEDK
01350 C
01351 C     Set pointers into IWORK.
01352 C
01353       PARAMETER (LML=1, LMU=2, LMTYPE=4, 
01354      *   LIWM=1, LMXORD=3, LJCALC=5, LPHASE=6, LK=7, LKOLD=8,
01355      *   LNS=9, LNSTL=10, LNST=11, LNRE=12, LNJE=13, LETF=14, LNCFN=15,
01356      *   LNCFL=16, LNIW=17, LNRW=18, LNNI=19, LNLI=20, LNPS=21,
01357      *   LNPD=22, LMITER=23, LMAXL=24, LKMP=25, LNRMAX=26, LLNWP=27,
01358      *   LLNIWP=28, LLOCWP=29, LLCIWP=30, LKPRIN=31,
01359      *   LMXNIT=32, LMXNJ=33, LMXNH=34, LLSOFF=35, LICNS=41)
01360 C
01361 C     Set pointers into RWORK.
01362 C
01363       PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4, LCJ=5, LCJOLD=6,
01364      *   LHOLD=7, LS=8, LROUND=9, LEPLI=10, LSQRN=11, LRSQRN=12,
01365      *   LEPCON=13, LSTOL=14, LEPIN=15,
01366      *   LALPHA=21, LBETA=27, LGAMMA=33, LPSI=39, LSIGMA=45, LDELTA=51)
01367 C
01368       SAVE LID, LENID, NONNEG
01369 C
01370 C
01371 C***FIRST EXECUTABLE STATEMENT  DDASPK
01372 C
01373 C
01374       IF(INFO(1).NE.0) GO TO 100
01375 C
01376 C-----------------------------------------------------------------------
01377 C     This block is executed for the initial call only.
01378 C     It contains checking of inputs and initializations.
01379 C-----------------------------------------------------------------------
01380 C
01381 C     First check INFO array to make sure all elements of INFO
01382 C     Are within the proper range.  (INFO(1) is checked later, because
01383 C     it must be tested on every call.) ITEMP holds the location
01384 C     within INFO which may be out of range.
01385 C
01386       DO 10 I=2,9
01387          ITEMP = I
01388          IF (INFO(I) .NE. 0 .AND. INFO(I) .NE. 1) GO TO 701
01389  10      CONTINUE
01390       ITEMP = 10
01391       IF(INFO(10).LT.0 .OR. INFO(10).GT.3) GO TO 701
01392       ITEMP = 11
01393       IF(INFO(11).LT.0 .OR. INFO(11).GT.2) GO TO 701
01394       DO 15 I=12,17
01395          ITEMP = I
01396          IF (INFO(I) .NE. 0 .AND. INFO(I) .NE. 1) GO TO 701
01397  15      CONTINUE
01398       ITEMP = 18
01399       IF(INFO(18).LT.0 .OR. INFO(18).GT.2) GO TO 701
01400 
01401 C
01402 C     Check NEQ to see if it is positive.
01403 C
01404       IF (NEQ .LE. 0) GO TO 702
01405 C
01406 C     Check and compute maximum order.
01407 C
01408       MXORD=5
01409       IF (INFO(9) .NE. 0) THEN
01410          MXORD=IWORK(LMXORD)
01411          IF (MXORD .LT. 1 .OR. MXORD .GT. 5) GO TO 703
01412          ENDIF
01413       IWORK(LMXORD)=MXORD
01414 C
01415 C     Set and/or check inputs for constraint checking (INFO(10) .NE. 0).
01416 C     Set values for ICNFLG, NONNEG, and pointer LID.
01417 C
01418       ICNFLG = 0
01419       NONNEG = 0
01420       LID = LICNS
01421       IF (INFO(10) .EQ. 0) GO TO 20
01422       IF (INFO(10) .EQ. 1) THEN
01423          ICNFLG = 1
01424          NONNEG = 0
01425          LID = LICNS + NEQ
01426       ELSEIF (INFO(10) .EQ. 2) THEN
01427          ICNFLG = 0
01428          NONNEG = 1
01429       ELSE
01430          ICNFLG = 1
01431          NONNEG = 1
01432          LID = LICNS + NEQ
01433       ENDIF
01434 C
01435  20   CONTINUE
01436 C
01437 C     Set and/or check inputs for Krylov solver (INFO(12) .NE. 0).
01438 C     If indicated, set default values for MAXL, KMP, NRMAX, and EPLI.
01439 C     Otherwise, verify inputs required for iterative solver.
01440 C
01441       IF (INFO(12) .EQ. 0) GO TO 25
01442 C
01443       IWORK(LMITER) = INFO(12)
01444       IF (INFO(13) .EQ. 0) THEN
01445          IWORK(LMAXL) = MIN(5,NEQ)
01446          IWORK(LKMP) = IWORK(LMAXL)
01447          IWORK(LNRMAX) = 5
01448          RWORK(LEPLI) = 0.05D0
01449       ELSE
01450          IF(IWORK(LMAXL) .LT. 1 .OR. IWORK(LMAXL) .GT. NEQ) GO TO 720
01451          IF(IWORK(LKMP) .LT. 1 .OR. IWORK(LKMP) .GT. IWORK(LMAXL))
01452      1      GO TO 721
01453          IF(IWORK(LNRMAX) .LT. 0) GO TO 722
01454          IF(RWORK(LEPLI).LE.0.0D0 .OR. RWORK(LEPLI).GE.1.0D0)GO TO 723
01455          ENDIF
01456 C
01457  25   CONTINUE
01458 C
01459 C     Set and/or check controls for the initial condition calculation
01460 C     (INFO(11) .GT. 0).  If indicated, set default values.
01461 C     Otherwise, verify inputs required for iterative solver.
01462 C
01463       IF (INFO(11) .EQ. 0) GO TO 30
01464       IF (INFO(17) .EQ. 0) THEN
01465         IWORK(LMXNIT) = 5
01466         IF (INFO(12) .GT. 0) IWORK(LMXNIT) = 15
01467         IWORK(LMXNJ) = 6
01468         IF (INFO(12) .GT. 0) IWORK(LMXNJ) = 2
01469         IWORK(LMXNH) = 5
01470         IWORK(LLSOFF) = 0
01471         RWORK(LEPIN) = 0.01D0
01472       ELSE
01473         IF (IWORK(LMXNIT) .LE. 0) GO TO 725
01474         IF (IWORK(LMXNJ) .LE. 0) GO TO 725
01475         IF (IWORK(LMXNH) .LE. 0) GO TO 725
01476         LSOFF = IWORK(LLSOFF)
01477         IF (LSOFF .LT. 0 .OR. LSOFF .GT. 1) GO TO 725
01478         IF (RWORK(LEPIN) .LE. 0.0D0) GO TO 725
01479         ENDIF
01480 C
01481  30   CONTINUE
01482 C
01483 C     Below is the computation and checking of the work array lengths
01484 C     LENIW and LENRW, using direct methods (INFO(12) = 0) or
01485 C     the Krylov methods (INFO(12) = 1).
01486 C
01487       LENIC = 0
01488       IF (INFO(10) .EQ. 1 .OR. INFO(10) .EQ. 3) LENIC = NEQ
01489       LENID = 0
01490       IF (INFO(11) .EQ. 1 .OR. INFO(16) .EQ. 1) LENID = NEQ
01491       IF (INFO(12) .EQ. 0) THEN
01492 C
01493 C        Compute MTYPE, etc.  Check ML and MU.
01494 C
01495          NCPHI = MAX(MXORD + 1, 4)
01496          IF(INFO(6).EQ.0) THEN 
01497             LENPD = NEQ**2
01498             LENRW = 50 + (NCPHI+3)*NEQ + LENPD
01499             IF(INFO(5).EQ.0) THEN
01500                IWORK(LMTYPE)=2
01501             ELSE
01502                IWORK(LMTYPE)=1
01503             ENDIF
01504          ELSE
01505             IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717
01506             IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718
01507             LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ
01508             IF(INFO(5).EQ.0) THEN
01509                IWORK(LMTYPE)=5
01510                MBAND=IWORK(LML)+IWORK(LMU)+1
01511                MSAVE=(NEQ/MBAND)+1
01512                LENRW = 50 + (NCPHI+3)*NEQ + LENPD + 2*MSAVE
01513             ELSE
01514                IWORK(LMTYPE)=4
01515                LENRW = 50 + (NCPHI+3)*NEQ + LENPD
01516             ENDIF
01517          ENDIF
01518 C
01519 C        Compute LENIW, LENWP, LENIWP.
01520 C
01521          LENIW = 40 + LENIC + LENID + NEQ
01522          LENWP = 0
01523          LENIWP = 0
01524 C
01525       ELSE IF (INFO(12) .EQ. 1)  THEN
01526          MAXL = IWORK(LMAXL)
01527          LENWP = IWORK(LLNWP)
01528          LENIWP = IWORK(LLNIWP)
01529          LENPD = (MAXL+3+MIN0(1,MAXL-IWORK(LKMP)))*NEQ
01530      1         + (MAXL+3)*MAXL + 1 + LENWP
01531          LENRW = 50 + (IWORK(LMXORD)+5)*NEQ + LENPD
01532          LENIW = 40 + LENIC + LENID + LENIWP
01533 C
01534       ENDIF
01535       IF(INFO(16) .NE. 0) LENRW = LENRW + NEQ
01536 C
01537 C     Check lengths of RWORK and IWORK.
01538 C
01539       IWORK(LNIW)=LENIW
01540       IWORK(LNRW)=LENRW
01541       IWORK(LNPD)=LENPD
01542       IWORK(LLOCWP) = LENPD-LENWP+1
01543       IF(LRW.LT.LENRW)GO TO 704
01544       IF(LIW.LT.LENIW)GO TO 705
01545 C
01546 C     Check ICNSTR for legality.
01547 C
01548       IF (LENIC .GT. 0) THEN
01549         DO 40 I = 1,NEQ
01550           ICI = IWORK(LICNS-1+I)
01551           IF (ICI .LT. -2 .OR. ICI .GT. 2) GO TO 726
01552  40       CONTINUE
01553         ENDIF
01554 C
01555 C     Check Y for consistency with constraints.
01556 C
01557       IF (LENIC .GT. 0) THEN
01558         CALL DCNST0(NEQ,Y,IWORK(LICNS),IRET)
01559         IF (IRET .NE. 0) GO TO 727
01560         ENDIF
01561 C
01562 C     Check ID for legality.
01563 C
01564       IF (LENID .GT. 0) THEN
01565         DO 50 I = 1,NEQ
01566           IDI = IWORK(LID-1+I)
01567           IF (IDI .NE. 1 .AND. IDI .NE. -1) GO TO 724
01568  50       CONTINUE
01569         ENDIF
01570 C
01571 C     Check to see that TOUT is different from T.
01572 C
01573       IF(TOUT .EQ. T)GO TO 719
01574 C
01575 C     Check HMAX.
01576 C
01577       IF(INFO(7) .NE. 0) THEN
01578          HMAX = RWORK(LHMAX)
01579          IF (HMAX .LE. 0.0D0) GO TO 710
01580          ENDIF
01581 C
01582 C     Initialize counters and other flags.
01583 C
01584       IWORK(LNST)=0
01585       IWORK(LNRE)=0
01586       IWORK(LNJE)=0
01587       IWORK(LETF)=0
01588       IWORK(LNCFN)=0
01589       IWORK(LNNI)=0
01590       IWORK(LNLI)=0
01591       IWORK(LNPS)=0
01592       IWORK(LNCFL)=0
01593       IWORK(LKPRIN)=INFO(18)
01594       IDID=1
01595       GO TO 200
01596 C
01597 C-----------------------------------------------------------------------
01598 C     This block is for continuation calls only.
01599 C     Here we check INFO(1), and if the last step was interrupted,
01600 C     we check whether appropriate action was taken.
01601 C-----------------------------------------------------------------------
01602 C
01603 100   CONTINUE
01604       IF(INFO(1).EQ.1)GO TO 110
01605       ITEMP = 1
01606       IF(INFO(1).NE.-1)GO TO 701
01607 C
01608 C     If we are here, the last step was interrupted by an error
01609 C     condition from DDSTP, and appropriate action was not taken.
01610 C     This is a fatal error.
01611 C
01612       MSG = 'DASPK--  THE LAST STEP TERMINATED WITH A NEGATIVE'
01613       CALL XERRWD(MSG,49,201,0,0,0,0,0,0.0D0,0.0D0)
01614       MSG = 'DASPK--  VALUE (=I1) OF IDID AND NO APPROPRIATE'
01615       CALL XERRWD(MSG,47,202,0,1,IDID,0,0,0.0D0,0.0D0)
01616       MSG = 'DASPK--  ACTION WAS TAKEN. RUN TERMINATED'
01617       CALL XERRWD(MSG,41,203,1,0,0,0,0,0.0D0,0.0D0)
01618       RETURN
01619 110   CONTINUE
01620 C
01621 C-----------------------------------------------------------------------
01622 C     This block is executed on all calls.
01623 C
01624 C     Counters are saved for later checks of performance.
01625 C     Then the error tolerance parameters are checked, and the
01626 C     work array pointers are set.
01627 C-----------------------------------------------------------------------
01628 C
01629 200   CONTINUE
01630 C
01631 C     Save counters for use later.
01632 C
01633       IWORK(LNSTL)=IWORK(LNST)
01634       NLI0 = IWORK(LNLI)
01635       NNI0 = IWORK(LNNI)
01636       NCFN0 = IWORK(LNCFN)
01637       NCFL0 = IWORK(LNCFL)
01638       NWARN = 0
01639 C
01640 C     Check RTOL and ATOL.
01641 C
01642       NZFLG = 0
01643       RTOLI = RTOL(1)
01644       ATOLI = ATOL(1)
01645       DO 210 I=1,NEQ
01646          IF (INFO(2) .EQ. 1) RTOLI = RTOL(I)
01647          IF (INFO(2) .EQ. 1) ATOLI = ATOL(I)
01648          IF (RTOLI .GT. 0.0D0 .OR. ATOLI .GT. 0.0D0) NZFLG = 1
01649          IF (RTOLI .LT. 0.0D0) GO TO 706
01650          IF (ATOLI .LT. 0.0D0) GO TO 707
01651 210      CONTINUE
01652       IF (NZFLG .EQ. 0) GO TO 708
01653 C
01654 C     Set pointers to RWORK and IWORK segments.
01655 C     For direct methods, SAVR is not used.
01656 C
01657       IWORK(LLCIWP) = LID + LENID
01658       LSAVR = LDELTA
01659       IF (INFO(12) .NE. 0) LSAVR = LDELTA + NEQ
01660       LE = LSAVR + NEQ
01661       LWT = LE + NEQ
01662       LVT = LWT
01663       IF (INFO(16) .NE. 0) LVT = LWT + NEQ
01664       LPHI = LVT + NEQ
01665       LWM = LPHI + (IWORK(LMXORD)+1)*NEQ
01666       IF (INFO(1) .EQ. 1) GO TO 400
01667 C
01668 C-----------------------------------------------------------------------
01669 C     This block is executed on the initial call only.
01670 C     Set the initial step size, the error weight vector, and PHI.
01671 C     Compute unknown initial components of Y and YPRIME, if requested.
01672 C-----------------------------------------------------------------------
01673 C
01674 300   CONTINUE
01675       TN=T
01676       IDID=1
01677 C
01678 C     Set error weight array WT and altered weight array VT.
01679 C
01680       CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
01681       CALL DINVWT(NEQ,RWORK(LWT),IER)
01682       IF (IER .NE. 0) GO TO 713
01683       IF (INFO(16) .NE. 0) THEN
01684         DO 305 I = 1, NEQ
01685  305      RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1)
01686         ENDIF
01687 C
01688 C     Compute unit roundoff and HMIN.
01689 C
01690       UROUND = D1MACH(4)
01691       RWORK(LROUND) = UROUND
01692       HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT))
01693 C
01694 C     Set/check STPTOL control for initial condition calculation.
01695 C     
01696       IF (INFO(11) .NE. 0) THEN
01697         IF( INFO(17) .EQ. 0) THEN
01698           RWORK(LSTOL) = UROUND**.6667D0
01699         ELSE
01700           IF (RWORK(LSTOL) .LE. 0.0D0) GO TO 725
01701           ENDIF
01702         ENDIF
01703 C
01704 C     Compute EPCON and square root of NEQ and its reciprocal, used
01705 C     inside iterative solver.
01706 C
01707       RWORK(LEPCON) = 0.33D0
01708       FLOATN = NEQ
01709       RWORK(LSQRN) = SQRT(FLOATN)
01710       RWORK(LRSQRN) = 1.D0/RWORK(LSQRN)
01711 C
01712 C     Check initial interval to see that it is long enough.
01713 C
01714       TDIST = ABS(TOUT - T)
01715       IF(TDIST .LT. HMIN) GO TO 714
01716 C
01717 C     Check H0, if this was input.
01718 C
01719       IF (INFO(8) .EQ. 0) GO TO 310
01720          H0 = RWORK(LH)
01721          IF ((TOUT - T)*H0 .LT. 0.0D0) GO TO 711
01722          IF (H0 .EQ. 0.0D0) GO TO 712
01723          GO TO 320
01724 310    CONTINUE
01725 C
01726 C     Compute initial stepsize, to be used by either
01727 C     DDSTP or DDASIC, depending on INFO(11).
01728 C
01729       H0 = 0.001D0*TDIST
01730       YPNORM = DDWNRM(NEQ,YPRIME,RWORK(LVT),RPAR,IPAR)
01731       IF (YPNORM .GT. 0.5D0/H0) H0 = 0.5D0/YPNORM
01732       H0 = SIGN(H0,TOUT-T)
01733 C
01734 C     Adjust H0 if necessary to meet HMAX bound.
01735 C
01736 320   IF (INFO(7) .EQ. 0) GO TO 330
01737          RH = ABS(H0)/RWORK(LHMAX)
01738          IF (RH .GT. 1.0D0) H0 = H0/RH
01739 C
01740 C     Check against TSTOP, if applicable.
01741 C
01742 330   IF (INFO(4) .EQ. 0) GO TO 340
01743          TSTOP = RWORK(LTSTOP)
01744          IF ((TSTOP - T)*H0 .LT. 0.0D0) GO TO 715
01745          IF ((T + H0 - TSTOP)*H0 .GT. 0.0D0) H0 = TSTOP - T
01746          IF ((TSTOP - TOUT)*H0 .LT. 0.0D0) GO TO 709
01747 C
01748 340   IF (INFO(11) .EQ. 0) GO TO 370
01749 C
01750 C     Compute unknown components of initial Y and YPRIME, depending
01751 C     on INFO(11) and INFO(12).  INFO(12) represents the nonlinear
01752 C     solver type (direct/Krylov).  Pass the name of the specific 
01753 C     nonlinear solver, depending on INFO(12).  The location of the work
01754 C     arrays SAVR, YIC, YPIC, PWK also differ in the two cases.
01755 C
01756       NWT = 1
01757       EPCONI = RWORK(LEPIN)*RWORK(LEPCON)
01758 350   IF (INFO(12) .EQ. 0) THEN
01759          LYIC = LPHI + 2*NEQ
01760          LYPIC = LYIC + NEQ
01761          LPWK = LYPIC
01762          CALL DDASIC(TN,Y,YPRIME,NEQ,INFO(11),IWORK(LID),
01763      *     RES,JAC,PSOL,H0,RWORK(LWT),NWT,IDID,RPAR,IPAR,
01764      *     RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE),
01765      *     RWORK(LYIC),RWORK(LYPIC),RWORK(LPWK),RWORK(LWM),IWORK(LIWM),
01766      *     HMIN,RWORK(LROUND),RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
01767      *     EPCONI,RWORK(LSTOL),INFO(15),ICNFLG,IWORK(LICNS),DDASID)
01768       ELSE IF (INFO(12) .EQ. 1) THEN
01769          LYIC = LWM
01770          LYPIC = LYIC + NEQ
01771          LPWK = LYPIC + NEQ
01772          CALL DDASIC(TN,Y,YPRIME,NEQ,INFO(11),IWORK(LID),
01773      *     RES,JAC,PSOL,H0,RWORK(LWT),NWT,IDID,RPAR,IPAR,
01774      *     RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE),
01775      *     RWORK(LYIC),RWORK(LYPIC),RWORK(LPWK),RWORK(LWM),IWORK(LIWM),
01776      *     HMIN,RWORK(LROUND),RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
01777      *     EPCONI,RWORK(LSTOL),INFO(15),ICNFLG,IWORK(LICNS),DDASIK)
01778       ENDIF
01779 C
01780       IF (IDID .LT. 0) GO TO 600
01781 C
01782 C     DDASIC was successful.  If this was the first call to DDASIC,
01783 C     update the WT array (with the current Y) and call it again.
01784 C
01785       IF (NWT .EQ. 2) GO TO 355
01786       NWT = 2
01787       CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
01788       CALL DINVWT(NEQ,RWORK(LWT),IER)
01789       IF (IER .NE. 0) GO TO 713
01790       GO TO 350
01791 C
01792 C     If INFO(14) = 1, return now with IDID = 4.
01793 C
01794 355   IF (INFO(14) .EQ. 1) THEN
01795         IDID = 4
01796         H = H0
01797         IF (INFO(11) .EQ. 1) RWORK(LHOLD) = H0
01798         GO TO 590
01799       ENDIF
01800 C
01801 C     Update the WT and VT arrays one more time, with the new Y.
01802 C
01803       CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
01804       CALL DINVWT(NEQ,RWORK(LWT),IER)
01805       IF (IER .NE. 0) GO TO 713
01806       IF (INFO(16) .NE. 0) THEN
01807         DO 357 I = 1, NEQ
01808  357      RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1)
01809         ENDIF
01810 C
01811 C     Reset the initial stepsize to be used by DDSTP.
01812 C     Use H0, if this was input.  Otherwise, recompute H0,
01813 C     and adjust it if necessary to meet HMAX bound.
01814 C
01815       IF (INFO(8) .NE. 0) THEN
01816          H0 = RWORK(LH)
01817          GO TO 360
01818          ENDIF
01819 C
01820       H0 = 0.001D0*TDIST
01821       YPNORM = DDWNRM(NEQ,YPRIME,RWORK(LVT),RPAR,IPAR)
01822       IF (YPNORM .GT. 0.5D0/H0) H0 = 0.5D0/YPNORM
01823       H0 = SIGN(H0,TOUT-T)
01824 C
01825 360   IF (INFO(7) .NE. 0) THEN
01826          RH = ABS(H0)/RWORK(LHMAX)
01827          IF (RH .GT. 1.0D0) H0 = H0/RH
01828          ENDIF
01829 C
01830 C     Check against TSTOP, if applicable.
01831 C
01832       IF (INFO(4) .NE. 0) THEN
01833          TSTOP = RWORK(LTSTOP)
01834          IF ((T + H0 - TSTOP)*H0 .GT. 0.0D0) H0 = TSTOP - T
01835          ENDIF
01836 C
01837 C     Load H and RWORK(LH) with H0.
01838 C
01839 370   H = H0
01840       RWORK(LH) = H
01841 C
01842 C     Load Y and H*YPRIME into PHI(*,1) and PHI(*,2).
01843 C
01844       ITEMP = LPHI + NEQ
01845       DO 380 I = 1,NEQ
01846          RWORK(LPHI + I - 1) = Y(I)
01847 380      RWORK(ITEMP + I - 1) = H*YPRIME(I)
01848 C
01849       GO TO 500
01850 C
01851 C-----------------------------------------------------------------------
01852 C     This block is for continuation calls only.
01853 C     Its purpose is to check stop conditions before taking a step.
01854 C     Adjust H if necessary to meet HMAX bound.
01855 C-----------------------------------------------------------------------
01856 C
01857 400   CONTINUE
01858       UROUND=RWORK(LROUND)
01859       DONE = .FALSE.
01860       TN=RWORK(LTN)
01861       H=RWORK(LH)
01862       IF(INFO(7) .EQ. 0) GO TO 410
01863          RH = ABS(H)/RWORK(LHMAX)
01864          IF(RH .GT. 1.0D0) H = H/RH
01865 410   CONTINUE
01866       IF(T .EQ. TOUT) GO TO 719
01867       IF((T - TOUT)*H .GT. 0.0D0) GO TO 711
01868       IF(INFO(4) .EQ. 1) GO TO 430
01869       IF(INFO(3) .EQ. 1) GO TO 420
01870       IF((TN-TOUT)*H.LT.0.0D0)GO TO 490
01871       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
01872      *  RWORK(LPHI),RWORK(LPSI))
01873       T=TOUT
01874       IDID = 3
01875       DONE = .TRUE.
01876       GO TO 490
01877 420   IF((TN-T)*H .LE. 0.0D0) GO TO 490
01878       IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425
01879       CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
01880      *  RWORK(LPHI),RWORK(LPSI))
01881       T = TN
01882       IDID = 1
01883       DONE = .TRUE.
01884       GO TO 490
01885 425   CONTINUE
01886       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
01887      *  RWORK(LPHI),RWORK(LPSI))
01888       T = TOUT
01889       IDID = 3
01890       DONE = .TRUE.
01891       GO TO 490
01892 430   IF(INFO(3) .EQ. 1) GO TO 440
01893       TSTOP=RWORK(LTSTOP)
01894       IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715
01895       IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709
01896       IF((TN-TOUT)*H.LT.0.0D0)GO TO 450
01897       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
01898      *   RWORK(LPHI),RWORK(LPSI))
01899       T=TOUT
01900       IDID = 3
01901       DONE = .TRUE.
01902       GO TO 490
01903 440   TSTOP = RWORK(LTSTOP)
01904       IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715
01905       IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709
01906       IF((TN-T)*H .LE. 0.0D0) GO TO 450
01907       IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445
01908       CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
01909      *  RWORK(LPHI),RWORK(LPSI))
01910       T = TN
01911       IDID = 1
01912       DONE = .TRUE.
01913       GO TO 490
01914 445   CONTINUE
01915       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
01916      *  RWORK(LPHI),RWORK(LPSI))
01917       T = TOUT
01918       IDID = 3
01919       DONE = .TRUE.
01920       GO TO 490
01921 450   CONTINUE
01922 C
01923 C     Check whether we are within roundoff of TSTOP.
01924 C
01925       IF(ABS(TN-TSTOP).GT.100.0D0*UROUND*
01926      *   (ABS(TN)+ABS(H)))GO TO 460
01927       CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD),
01928      *  RWORK(LPHI),RWORK(LPSI))
01929       IDID=2
01930       T=TSTOP
01931       DONE = .TRUE.
01932       GO TO 490
01933 460   TNEXT=TN+H
01934       IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490
01935       H=TSTOP-TN
01936       RWORK(LH)=H
01937 C
01938 490   IF (DONE) GO TO 590
01939 C
01940 C-----------------------------------------------------------------------
01941 C     The next block contains the call to the one-step integrator DDSTP.
01942 C     This is a looping point for the integration steps.
01943 C     Check for too many steps.
01944 C     Check for poor Newton/Krylov performance.
01945 C     Update WT.  Check for too much accuracy requested.
01946 C     Compute minimum stepsize.
01947 C-----------------------------------------------------------------------
01948 C
01949 500   CONTINUE
01950 C
01951 C     Check for too many steps.
01952 C
01953       IF((IWORK(LNST)-IWORK(LNSTL)).LT.500) GO TO 505
01954            IDID=-1
01955            GO TO 527
01956 C
01957 C Check for poor Newton/Krylov performance.
01958 C
01959 505   IF (INFO(12) .EQ. 0) GO TO 510
01960       NSTD = IWORK(LNST) - IWORK(LNSTL)
01961       NNID = IWORK(LNNI) - NNI0
01962       IF (NSTD .LT. 10 .OR. NNID .EQ. 0) GO TO 510
01963       AVLIN = REAL(IWORK(LNLI) - NLI0)/REAL(NNID)
01964       RCFN = REAL(IWORK(LNCFN) - NCFN0)/REAL(NSTD)
01965       RCFL = REAL(IWORK(LNCFL) - NCFL0)/REAL(NNID)
01966       FMAXL = IWORK(LMAXL)
01967       LAVL = AVLIN .GT. FMAXL
01968       LCFN = RCFN .GT. 0.9D0
01969       LCFL = RCFL .GT. 0.9D0
01970       LWARN = LAVL .OR. LCFN .OR. LCFL
01971       IF (.NOT.LWARN) GO TO 510
01972       NWARN = NWARN + 1
01973       IF (NWARN .GT. 10) GO TO 510
01974       IF (LAVL) THEN
01975         MSG = 'DASPK-- Warning. Poor iterative algorithm performance   '
01976         CALL XERRWD (MSG, 56, 501, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01977         MSG = '      at T = R1. Average no. of linear iterations = R2  '
01978         CALL XERRWD (MSG, 56, 501, 0, 0, 0, 0, 2, TN, AVLIN)
01979         ENDIF
01980       IF (LCFN) THEN
01981         MSG = 'DASPK-- Warning. Poor iterative algorithm performance   '
01982         CALL XERRWD (MSG, 56, 502, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01983         MSG = '      at T = R1. Nonlinear convergence failure rate = R2'
01984         CALL XERRWD (MSG, 56, 502, 0, 0, 0, 0, 2, TN, RCFN)
01985         ENDIF
01986       IF (LCFL) THEN
01987         MSG = 'DASPK-- Warning. Poor iterative algorithm performance   '
01988         CALL XERRWD (MSG, 56, 503, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
01989         MSG = '      at T = R1. Linear convergence failure rate = R2   '
01990         CALL XERRWD (MSG, 56, 503, 0, 0, 0, 0, 2, TN, RCFL)
01991         ENDIF
01992 C
01993 C     Update WT and VT, if this is not the first call.
01994 C
01995 510   CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),RWORK(LWT),
01996      *            RPAR,IPAR)
01997       CALL DINVWT(NEQ,RWORK(LWT),IER)
01998       IF (IER .NE. 0) THEN
01999         IDID = -3
02000         GO TO 527
02001         ENDIF
02002       IF (INFO(16) .NE. 0) THEN
02003         DO 515 I = 1, NEQ
02004  515      RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1)
02005         ENDIF
02006 C
02007 C     Test for too much accuracy requested.
02008 C
02009       R = DDWNRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*100.0D0*UROUND
02010       IF (R .LE. 1.0D0) GO TO 525
02011 C
02012 C     Multiply RTOL and ATOL by R and return.
02013 C
02014       IF(INFO(2).EQ.1)GO TO 523
02015            RTOL(1)=R*RTOL(1)
02016            ATOL(1)=R*ATOL(1)
02017            IDID=-2
02018            GO TO 527
02019 523   DO 524 I=1,NEQ
02020            RTOL(I)=R*RTOL(I)
02021 524        ATOL(I)=R*ATOL(I)
02022       IDID=-2
02023       GO TO 527
02024 525   CONTINUE
02025 C
02026 C     Compute minimum stepsize.
02027 C
02028       HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT))
02029 C
02030 C     Test H vs. HMAX
02031       IF (INFO(7) .NE. 0) THEN
02032          RH = ABS(H)/RWORK(LHMAX)
02033          IF (RH .GT. 1.0D0) H = H/RH
02034          ENDIF
02035 C
02036 C     Call the one-step integrator.
02037 C     Note that INFO(12) represents the nonlinear solver type.
02038 C     Pass the required nonlinear solver, depending upon INFO(12).
02039 C
02040       IF (INFO(12) .EQ. 0) THEN
02041          CALL DDSTP(TN,Y,YPRIME,NEQ,
02042      *      RES,JAC,PSOL,H,RWORK(LWT),RWORK(LVT),INFO(1),IDID,RPAR,IPAR,
02043      *      RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE),
02044      *      RWORK(LWM),IWORK(LIWM),
02045      *      RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),
02046      *      RWORK(LPSI),RWORK(LSIGMA),
02047      *      RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),RWORK(LS),HMIN,
02048      *      RWORK(LROUND), RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
02049      *      RWORK(LEPCON), IWORK(LPHASE),IWORK(LJCALC),INFO(15),
02050      *      IWORK(LK), IWORK(LKOLD),IWORK(LNS),NONNEG,INFO(12),
02051      *      DNEDD)
02052       ELSE IF (INFO(12) .EQ. 1) THEN
02053          CALL DDSTP(TN,Y,YPRIME,NEQ,
02054      *      RES,JAC,PSOL,H,RWORK(LWT),RWORK(LVT),INFO(1),IDID,RPAR,IPAR,
02055      *      RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE),
02056      *      RWORK(LWM),IWORK(LIWM),
02057      *      RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),
02058      *      RWORK(LPSI),RWORK(LSIGMA),
02059      *      RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),RWORK(LS),HMIN,
02060      *      RWORK(LROUND), RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
02061      *      RWORK(LEPCON), IWORK(LPHASE),IWORK(LJCALC),INFO(15),
02062      *      IWORK(LK), IWORK(LKOLD),IWORK(LNS),NONNEG,INFO(12),
02063      *      DNEDK)
02064       ENDIF
02065 C
02066 527   IF(IDID.LT.0)GO TO 600
02067 C
02068 C-----------------------------------------------------------------------
02069 C     This block handles the case of a successful return from DDSTP
02070 C     (IDID=1).  Test for stop conditions.
02071 C-----------------------------------------------------------------------
02072 C
02073       IF(INFO(4).NE.0)GO TO 540
02074            IF(INFO(3).NE.0)GO TO 530
02075              IF((TN-TOUT)*H.LT.0.0D0)GO TO 500
02076              CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
02077      *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
02078              IDID=3
02079              T=TOUT
02080              GO TO 580
02081 530          IF((TN-TOUT)*H.GE.0.0D0)GO TO 535
02082              T=TN
02083              IDID=1
02084              GO TO 580
02085 535          CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
02086      *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
02087              IDID=3
02088              T=TOUT
02089              GO TO 580
02090 540   IF(INFO(3).NE.0)GO TO 550
02091       IF((TN-TOUT)*H.LT.0.0D0)GO TO 542
02092          CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
02093      *     IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
02094          T=TOUT
02095          IDID=3
02096          GO TO 580
02097 542   IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*
02098      *   (ABS(TN)+ABS(H)))GO TO 545
02099       TNEXT=TN+H
02100       IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500
02101       H=TSTOP-TN
02102       GO TO 500
02103 545   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
02104      *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
02105       IDID=2
02106       T=TSTOP
02107       GO TO 580
02108 550   IF((TN-TOUT)*H.GE.0.0D0)GO TO 555
02109       IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552
02110       T=TN
02111       IDID=1
02112       GO TO 580
02113 552   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
02114      *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
02115       IDID=2
02116       T=TSTOP
02117       GO TO 580
02118 555   CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
02119      *   IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
02120       T=TOUT
02121       IDID=3
02122 580   CONTINUE
02123 C
02124 C-----------------------------------------------------------------------
02125 C     All successful returns from DDASPK are made from this block.
02126 C-----------------------------------------------------------------------
02127 C
02128 590   CONTINUE
02129       RWORK(LTN)=TN
02130       RWORK(LH)=H
02131       RETURN
02132 C
02133 C-----------------------------------------------------------------------
02134 C     This block handles all unsuccessful returns other than for
02135 C     illegal input.
02136 C-----------------------------------------------------------------------
02137 C
02138 600   CONTINUE
02139       ITEMP = -IDID
02140       GO TO (610,620,630,700,655,640,650,660,670,675,
02141      *  680,685,690,695), ITEMP
02142 C
02143 C     The maximum number of steps was taken before
02144 C     reaching tout.
02145 C
02146 610   MSG = 'DASPK--  AT CURRENT T (=R1)  500 STEPS'
02147       CALL XERRWD(MSG,38,610,0,0,0,0,1,TN,0.0D0)
02148       MSG = 'DASPK--  TAKEN ON THIS CALL BEFORE REACHING TOUT'
02149       CALL XERRWD(MSG,48,611,0,0,0,0,0,0.0D0,0.0D0)
02150       GO TO 700
02151 C
02152 C     Too much accuracy for machine precision.
02153 C
02154 620   MSG = 'DASPK--  AT T (=R1) TOO MUCH ACCURACY REQUESTED'
02155       CALL XERRWD(MSG,47,620,0,0,0,0,1,TN,0.0D0)
02156       MSG = 'DASPK--  FOR PRECISION OF MACHINE. RTOL AND ATOL'
02157       CALL XERRWD(MSG,48,621,0,0,0,0,0,0.0D0,0.0D0)
02158       MSG = 'DASPK--  WERE INCREASED TO APPROPRIATE VALUES'
02159       CALL XERRWD(MSG,45,622,0,0,0,0,0,0.0D0,0.0D0)
02160       GO TO 700
02161 C
02162 C     WT(I) .LE. 0.0D0 for some I (not at start of problem).
02163 C
02164 630   MSG = 'DASPK--  AT T (=R1) SOME ELEMENT OF WT'
02165       CALL XERRWD(MSG,38,630,0,0,0,0,1,TN,0.0D0)
02166       MSG = 'DASPK--  HAS BECOME .LE. 0.0'
02167       CALL XERRWD(MSG,28,631,0,0,0,0,0,0.0D0,0.0D0)
02168       GO TO 700
02169 C
02170 C     Error test failed repeatedly or with H=HMIN.
02171 C
02172 640   MSG = 'DASPK--  AT T (=R1) AND STEPSIZE H (=R2) THE'
02173       CALL XERRWD(MSG,44,640,0,0,0,0,2,TN,H)
02174       MSG='DASPK--  ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN'
02175       CALL XERRWD(MSG,57,641,0,0,0,0,0,0.0D0,0.0D0)
02176       GO TO 700
02177 C
02178 C     Nonlinear solver failed to converge repeatedly or with H=HMIN.
02179 C
02180 650   MSG = 'DASPK--  AT T (=R1) AND STEPSIZE H (=R2) THE'
02181       CALL XERRWD(MSG,44,650,0,0,0,0,2,TN,H)
02182       MSG = 'DASPK--  NONLINEAR SOLVER FAILED TO CONVERGE'
02183       CALL XERRWD(MSG,44,651,0,0,0,0,0,0.0D0,0.0D0)
02184       MSG = 'DASPK--  REPEATEDLY OR WITH ABS(H)=HMIN'
02185       CALL XERRWD(MSG,40,652,0,0,0,0,0,0.0D0,0.0D0)
02186       GO TO 700
02187 C
02188 C     The preconditioner had repeated failures.
02189 C
02190 655   MSG = 'DASPK--  AT T (=R1) AND STEPSIZE H (=R2) THE'
02191       CALL XERRWD(MSG,44,655,0,0,0,0,2,TN,H)
02192       MSG = 'DASPK--  PRECONDITIONER HAD REPEATED FAILURES.'
02193       CALL XERRWD(MSG,46,656,0,0,0,0,0,0.0D0,0.0D0)
02194       GO TO 700
02195 C
02196 C     The iteration matrix is singular.
02197 C
02198 660   MSG = 'DASPK--  AT T (=R1) AND STEPSIZE H (=R2) THE'
02199       CALL XERRWD(MSG,44,660,0,0,0,0,2,TN,H)
02200       MSG = 'DASPK--  ITERATION MATRIX IS SINGULAR.'
02201       CALL XERRWD(MSG,38,661,0,0,0,0,0,0.0D0,0.0D0)
02202       GO TO 700
02203 C
02204 C     Nonlinear system failure preceded by error test failures.
02205 C
02206 670   MSG = 'DASPK--  AT T (=R1) AND STEPSIZE H (=R2) THE'
02207       CALL XERRWD(MSG,44,670,0,0,0,0,2,TN,H)
02208       MSG = 'DASPK--  NONLINEAR SOLVER COULD NOT CONVERGE.'
02209       CALL XERRWD(MSG,45,671,0,0,0,0,0,0.0D0,0.0D0)
02210       MSG = 'DASPK--  ALSO, THE ERROR TEST FAILED REPEATEDLY.'
02211       CALL XERRWD(MSG,49,672,0,0,0,0,0,0.0D0,0.0D0)
02212       GO TO 700
02213 C
02214 C     Nonlinear system failure because IRES = -1.
02215 C
02216 675   MSG = 'DASPK--  AT T (=R1) AND STEPSIZE H (=R2) THE'
02217       CALL XERRWD(MSG,44,675,0,0,0,0,2,TN,H)
02218       MSG = 'DASPK--  NONLINEAR SYSTEM SOLVER COULD NOT CONVERGE'
02219       CALL XERRWD(MSG,51,676,0,0,0,0,0,0.0D0,0.0D0)
02220       MSG = 'DASPK--  BECAUSE IRES WAS EQUAL TO MINUS ONE'
02221       CALL XERRWD(MSG,44,677,0,0,0,0,0,0.0D0,0.0D0)
02222       GO TO 700
02223 C
02224 C     Failure because IRES = -2.
02225 C
02226 680   MSG = 'DASPK--  AT T (=R1) AND STEPSIZE H (=R2)'
02227       CALL XERRWD(MSG,40,680,0,0,0,0,2,TN,H)
02228       MSG = 'DASPK--  IRES WAS EQUAL TO MINUS TWO'
02229       CALL XERRWD(MSG,36,681,0,0,0,0,0,0.0D0,0.0D0)
02230       GO TO 700
02231 C
02232 C     Failed to compute initial YPRIME.
02233 C
02234 685   MSG = 'DASPK--  AT T (=R1) AND STEPSIZE H (=R2) THE'
02235       CALL XERRWD(MSG,44,685,0,0,0,0,0,0.0D0,0.0D0)
02236       MSG = 'DASPK--  INITIAL (Y,YPRIME) COULD NOT BE COMPUTED'
02237       CALL XERRWD(MSG,49,686,0,0,0,0,2,TN,H0)
02238       GO TO 700
02239 C
02240 C     Failure because IER was negative from PSOL.
02241 C
02242 690   MSG = 'DASPK--  AT T (=R1) AND STEPSIZE H (=R2)'
02243       CALL XERRWD(MSG,40,690,0,0,0,0,2,TN,H)
02244       MSG = 'DASPK--  IER WAS NEGATIVE FROM PSOL'
02245       CALL XERRWD(MSG,35,691,0,0,0,0,0,0.0D0,0.0D0)
02246       GO TO 700
02247 C
02248 C     Failure because the linear system solver could not converge.
02249 C
02250 695   MSG = 'DASPK--  AT T (=R1) AND STEPSIZE H (=R2) THE'
02251       CALL XERRWD(MSG,44,695,0,0,0,0,2,TN,H)
02252       MSG = 'DASPK--  LINEAR SYSTEM SOLVER COULD NOT CONVERGE.'
02253       CALL XERRWD(MSG,50,696,0,0,0,0,0,0.0D0,0.0D0)
02254       GO TO 700
02255 C
02256 C
02257 700   CONTINUE
02258       INFO(1)=-1
02259       T=TN
02260       RWORK(LTN)=TN
02261       RWORK(LH)=H
02262       RETURN
02263 C
02264 C-----------------------------------------------------------------------
02265 C     This block handles all error returns due to illegal input,
02266 C     as detected before calling DDSTP.
02267 C     First the error message routine is called.  If this happens
02268 C     twice in succession, execution is terminated.
02269 C-----------------------------------------------------------------------
02270 C
02271 701   MSG = 'DASPK--  ELEMENT (=I1) OF INFO VECTOR IS NOT VALID'
02272       CALL XERRWD(MSG,50,1,0,1,ITEMP,0,0,0.0D0,0.0D0)
02273       GO TO 750
02274 702   MSG = 'DASPK--  NEQ (=I1) .LE. 0'
02275       CALL XERRWD(MSG,25,2,0,1,NEQ,0,0,0.0D0,0.0D0)
02276       GO TO 750
02277 703   MSG = 'DASPK--  MAXORD (=I1) NOT IN RANGE'
02278       CALL XERRWD(MSG,34,3,0,1,MXORD,0,0,0.0D0,0.0D0)
02279       GO TO 750
02280 704   MSG='DASPK--  RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)'
02281       CALL XERRWD(MSG,60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0)
02282       GO TO 750
02283 705   MSG='DASPK--  IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)'
02284       CALL XERRWD(MSG,60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0)
02285       GO TO 750
02286 706   MSG = 'DASPK--  SOME ELEMENT OF RTOL IS .LT. 0'
02287       CALL XERRWD(MSG,39,6,0,0,0,0,0,0.0D0,0.0D0)
02288       GO TO 750
02289 707   MSG = 'DASPK--  SOME ELEMENT OF ATOL IS .LT. 0'
02290       CALL XERRWD(MSG,39,7,0,0,0,0,0,0.0D0,0.0D0)
02291       GO TO 750
02292 708   MSG = 'DASPK--  ALL ELEMENTS OF RTOL AND ATOL ARE ZERO'
02293       CALL XERRWD(MSG,47,8,0,0,0,0,0,0.0D0,0.0D0)
02294       GO TO 750
02295 709   MSG='DASPK--  INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)'
02296       CALL XERRWD(MSG,54,9,0,0,0,0,2,TSTOP,TOUT)
02297       GO TO 750
02298 710   MSG = 'DASPK--  HMAX (=R1) .LT. 0.0'
02299       CALL XERRWD(MSG,28,10,0,0,0,0,1,HMAX,0.0D0)
02300       GO TO 750
02301 711   MSG = 'DASPK--  TOUT (=R1) BEHIND T (=R2)'
02302       CALL XERRWD(MSG,34,11,0,0,0,0,2,TOUT,T)
02303       GO TO 750
02304 712   MSG = 'DASPK--  INFO(8)=1 AND H0=0.0'
02305       CALL XERRWD(MSG,29,12,0,0,0,0,0,0.0D0,0.0D0)
02306       GO TO 750
02307 713   MSG = 'DASPK--  SOME ELEMENT OF WT IS .LE. 0.0'
02308       CALL XERRWD(MSG,39,13,0,0,0,0,0,0.0D0,0.0D0)
02309       GO TO 750
02310 714   MSG='DASPK-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION'
02311       CALL XERRWD(MSG,60,14,0,0,0,0,2,TOUT,T)
02312       GO TO 750
02313 715   MSG = 'DASPK--  INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)'
02314       CALL XERRWD(MSG,49,15,0,0,0,0,2,TSTOP,T)
02315       GO TO 750
02316 717   MSG = 'DASPK--  ML (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ'
02317       CALL XERRWD(MSG,52,17,0,1,IWORK(LML),0,0,0.0D0,0.0D0)
02318       GO TO 750
02319 718   MSG = 'DASPK--  MU (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ'
02320       CALL XERRWD(MSG,52,18,0,1,IWORK(LMU),0,0,0.0D0,0.0D0)
02321       GO TO 750
02322 719   MSG = 'DASPK--  TOUT (=R1) IS EQUAL TO T (=R2)'
02323       CALL XERRWD(MSG,39,19,0,0,0,0,2,TOUT,T)
02324       GO TO 750
02325 720   MSG = 'DASPK--  MAXL (=I1) ILLEGAL. EITHER .LT. 1 OR .GT. NEQ'
02326       CALL XERRWD(MSG,54,20,0,1,IWORK(LMAXL),0,0,0.0D0,0.0D0)
02327       GO TO 750
02328 721   MSG = 'DASPK--  KMP (=I1) ILLEGAL. EITHER .LT. 1 OR .GT. MAXL'
02329       CALL XERRWD(MSG,54,21,0,1,IWORK(LKMP),0,0,0.0D0,0.0D0)
02330       GO TO 750
02331 722   MSG = 'DASPK--  NRMAX (=I1) ILLEGAL. .LT. 0'
02332       CALL XERRWD(MSG,36,22,0,1,IWORK(LNRMAX),0,0,0.0D0,0.0D0)
02333       GO TO 750
02334 723   MSG = 'DASPK--  EPLI (=R1) ILLEGAL. EITHER .LE. 0.D0 OR .GE. 1.D0'
02335       CALL XERRWD(MSG,58,23,0,0,0,0,1,RWORK(LEPLI),0.0D0)
02336       GO TO 750
02337 724   MSG = 'DASPK--  ILLEGAL IWORK VALUE FOR INFO(11) .NE. 0'
02338       CALL XERRWD(MSG,48,24,0,0,0,0,0,0.0D0,0.0D0)
02339       GO TO 750
02340 725   MSG = 'DASPK--  ONE OF THE INPUTS FOR INFO(17) = 1 IS ILLEGAL'
02341       CALL XERRWD(MSG,54,25,0,0,0,0,0,0.0D0,0.0D0)
02342       GO TO 750
02343 726   MSG = 'DASPK--  ILLEGAL IWORK VALUE FOR INFO(10) .NE. 0'
02344       CALL XERRWD(MSG,48,26,0,0,0,0,0,0.0D0,0.0D0)
02345       GO TO 750
02346 727   MSG = 'DASPK--  Y(I) AND IWORK(40+I) (I=I1) INCONSISTENT'
02347       CALL XERRWD(MSG,49,27,0,1,IRET,0,0,0.0D0,0.0D0)
02348       GO TO 750
02349 750   IF(INFO(1).EQ.-1) GO TO 760
02350       INFO(1)=-1
02351       IDID=-33
02352       RETURN
02353 760   MSG = 'DASPK--  REPEATED OCCURRENCES OF ILLEGAL INPUT'
02354       CALL XERRWD(MSG,46,701,0,0,0,0,0,0.0D0,0.0D0)
02355 770   MSG = 'DASPK--  RUN TERMINATED. APPARENT INFINITE LOOP'
02356       CALL XERRWD(MSG,47,702,1,0,0,0,0,0.0D0,0.0D0)
02357       RETURN
02358 C
02359 C------END OF SUBROUTINE DDASPK-----------------------------------------
02360       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines