Logo Search packages:      
Sourcecode: octave3.2 version File versions


     *  G,NG,JROOT)
C***DATE WRITTEN   821001   (YYMMDD)
C             L - 316, P.O. Box 808,
C             LIVERMORE, CA.    94550
C***PURPOSE  This code solves a system of differential/algebraic
C            equations of the form F(T,Y,YPRIME) = 0.
C *Usage:
C     *   JROOT(NG)
C     *   RWORK(LRW), RPAR
C *Arguments:
C  RES:EXT  This is a subroutine which you provide to define the
C           differential/algebraic system.
C  NEQ:IN  This is the number of equations to be solved.
C  T:INOUT  This is the current value of the independent variable.
C  Y(*):INOUT  This array contains the solution components at T.
C  YPRIME(*):INOUT  This array contains the derivatives of the solution
C                   components at T.
C  TOUT:IN  This is a point at which a solution is desired.
C  INFO(N):IN  The basic task of the code is to solve the system from T
C              to TOUT and return an answer at TOUT.  INFO is an integer
C              array which is used to communicate exactly how you want
C              this task to be carried out.  N must be greater than or
C              equal to 15.
C  RTOL,ATOL:INOUT  These quantities represent absolute and relative
C                   error tolerances which you provide to indicate how
C                   accurately you wish the solution to be computed.
C                   You may choose them to be both scalars or else
C                   both vectors.
C  IDID:OUT  This scalar quantity is an indicator reporting what the
C            code did.  You must monitor this integer variable to decide
C            what action to take next.
C  RWORK:WORK  A real work array of length LRW which provides the
C               code with needed storage space.
C  LRW:IN  The length of RWORK.
C  IWORK:WORK  An integer work array of length LIW which probides the
C               code with needed storage space.
C  LIW:IN  The length of IWORK.
C  RPAR,IPAR:IN  These are real and integer parameter arrays which
C                you can use for communication between your calling
C                program and the RES subroutine (and the JAC subroutine)
C  JAC:EXT  This is the name of a subroutine which you may choose to
C           provide for defining a matrix of partial derivatives
C           described below.
C  G  This is the name of the subroutine for defining
C     constraint functions, G(T,Y), whose roots are desired
C     during the integration.  This name must be declared
C     external in the calling program.
C  NG  This is the number of constraint functions G(I).
C      If there are none, set NG=0, and pass a dummy name
C      for G.
C  JROOT  This is an integer array of length NG for output
C         of root information.
C *Description
C  Subroutine DDASRT uses the backward differentiation formulas of
C  orders one through five to solve a system of the above form for Y and
C  YPRIME.  Values for Y and YPRIME at the initial time must be given as
C  input.  These values must be consistent, (that is, if T,Y,YPRIME are
C  the given initial values, they must satisfy F(T,Y,YPRIME) = 0.).  The
C  subroutine solves the system from T to TOUT.
C  It is easy to continue the solution to get results at additional
C  TOUT.  This is the interval mode of operation.  Intermediate results
C  can also be obtained easily by using the intermediate-output
C  capability.  If DDASRT detects a sign-change in G(T,Y), then
C  it will return the intermediate value of T and Y for which
C  G(T,Y) = 0.
C  The first call of the code is defined to be the start of each new
C  problem. Read through the descriptions of all the following items,
C  provide sufficient storage space for designated arrays, set
C  appropriate variables for the initialization of the problem, and
C  give information about how you want the problem to be solved.
C  RES -- Provide a subroutine of the form
C         to define the system of differential/algebraic
C         equations which is to be solved. For the given values
C         of T,Y and YPRIME, the subroutine should
C         return the residual of the defferential/algebraic
C         system
C             DELTA = F(T,Y,YPRIME)
C         (DELTA(*) is a vector of length NEQ which is
C         output for RES.)
C         Subroutine RES must not alter T,Y or YPRIME.
C         You must declare the name RES in an external
C         statement in your program that calls DDASRT.
C         You must dimension Y,YPRIME and DELTA in RES.
C         IRES is an integer flag which is always equal to
C         zero on input. Subroutine RES should alter IRES
C         only if it encounters an illegal value of Y or
C         a stop condition. Set IRES = -1 if an input value
C         is illegal, and DDASRT will try to solve the problem
C         without getting IRES = -1. If IRES = -2, DDASRT
C         will return control to the calling program
C         with IDID = -11.
C         RPAR and IPAR are real and integer parameter arrays which
C         you can use for communication between your calling program
C         and subroutine RES. They are not altered by DDASRT. If you
C         do not need RPAR or IPAR, ignore these parameters by treat-
C         ing them as dummy arguments. If you do choose to use them,
C         dimension them in your calling program and in RES as arrays
C         of appropriate length.
C  NEQ -- Set it to the number of differential equations.
C         (NEQ .GE. 1)
C  T -- Set it to the initial point of the integration.
C       T must be defined as a variable.
C  Y(*) -- Set this vector to the initial values of the NEQ solution
C          components at the initial point. You must dimension Y of
C          length at least NEQ in your calling program.
C  YPRIME(*) -- Set this vector to the initial values of
C               the NEQ first derivatives of the solution
C               components at the initial point. You
C               must dimension YPRIME at least NEQ
C               in your calling program. If you do not
C               know initial values of some of the solution
C               components, see the explanation of INFO(11).
C  TOUT - Set it to the first point at which a solution
C         is desired. You can not take TOUT = T.
C         integration either forward in T (TOUT .GT. T) or
C         backward in T (TOUT .LT. T) is permitted.
C         The code advances the solution from T to TOUT using
C         step sizes which are automatically selected so as to
C         achieve the desired accuracy. If you wish, the code will
C         return with the solution and its derivative at
C         intermediate steps (intermediate-output mode) so that
C         you can monitor them, but you still must provide TOUT in
C         accord with the basic aim of the code.
C         the first step taken by the code is a critical one
C         because it must reflect how fast the solution changes near
C         the initial point. The code automatically selects an
C         initial step size which is practically always suitable for
C         the problem. By using the fact that the code will not step
C         past TOUT in the first step, you could, if necessary,
C         restrict the length of the initial step size.
C         For some problems it may not be permissable to integrate
C         past a point TSTOP because a discontinuity occurs there
C         or the solution or its derivative is not defined beyond
C         TSTOP. When you have declared a TSTOP point (SEE INFO(4)
C         and RWORK(1)), you have told the code not to integrate
C         past TSTOP. In this case any TOUT beyond TSTOP is invalid
C         input.
C  INFO(*) - Use the INFO array to give the code more details about
C            how you want your problem solved. This array should be
C            dimensioned of length 15, though DDASRT uses
C            only the first twelve entries. You must respond to all of
C            the following items which are arranged as questions. The
C            simplest use of the code corresponds to answering all
C            questions as yes, i.e. setting all entries of INFO to 0.
C       INFO(1) - This parameter enables the code to initialize
C              itself. You must set it to indicate the start of every
C              new problem.
C          **** Is this the first call for this problem ...
C                Yes - Set INFO(1) = 0
C                 No - Not applicable here.
C                      See below for continuation calls.  ****
C       INFO(2) - How much accuracy you want of your solution
C              is specified by the error tolerances RTOL and ATOL.
C              The simplest use is to take them both to be scalars.
C              To obtain more flexibility, they can both be vectors.
C              The code must be told your choice.
C          **** Are both error tolerances RTOL, ATOL scalars ...
C                Yes - Set INFO(2) = 0
C                      and input scalars for both RTOL and ATOL
C                 No - Set INFO(2) = 1
C                      and input arrays for both RTOL and ATOL ****
C       INFO(3) - The code integrates from T in the direction
C              of TOUT by steps. If you wish, it will return the
C              computed solution and derivative at the next
C              intermediate step (the intermediate-output mode) or
C              TOUT, whichever comes first. This is a good way to
C              proceed if you want to see the behavior of the solution.
C              If you must have solutions at a great many specific
C              TOUT points, this code will compute them efficiently.
C          **** Do you want the solution only at
C                TOUT (and not at the next intermediate step) ...
C                 Yes - Set INFO(3) = 0
C                  No - Set INFO(3) = 1 ****
C       INFO(4) - To handle solutions at a great many specific
C              values TOUT efficiently, this code may integrate past
C              TOUT and interpolate to obtain the result at TOUT.
C              Sometimes it is not possible to integrate beyond some
C              point TSTOP because the equation changes there or it is
C              not defined past TSTOP. Then you must tell the code
C              not to go past.
C           **** Can the integration be carried out without any
C                restrictions on the independent variable T ...
C                 Yes - Set INFO(4)=0
C                  No - Set INFO(4)=1
C                       and define the stopping point TSTOP by
C                       setting RWORK(1)=TSTOP ****
C       INFO(5) - To solve differential/algebraic problems it is
C              necessary to use a matrix of partial derivatives of the
C              system of differential equations. If you do not
C              provide a subroutine to evaluate it analytically (see
C              description of the item JAC in the call list), it will
C              be approximated by numerical differencing in this code.
C              although it is less trouble for you to have the code
C              compute partial derivatives by numerical differencing,
C              the solution will be more reliable if you provide the
C              derivatives via JAC. Sometimes numerical differencing
C              is cheaper than evaluating derivatives in JAC and
C              sometimes it is not - this depends on your problem.
C           **** Do you want the code to evaluate the partial
C                derivatives automatically by numerical differences ...
C                   Yes - Set INFO(5)=0
C                    No - Set INFO(5)=1
C                  and provide subroutine JAC for evaluating the
C                  matrix of partial derivatives ****
C       INFO(6) - DDASRT will perform much better if the matrix of
C              partial derivatives, DG/DY + CJ*DG/DYPRIME,
C              (here CJ is a scalar determined by DDASRT)
C              is banded and the code is told this. In this
C              case, the storage needed will be greatly reduced,
C              numerical differencing will be performed much cheaper,
C              and a number of important algorithms will execute much
C              faster. The differential equation is said to have
C              half-bandwidths ML (lower) and MU (upper) if equation i
C              involves only unknowns Y(J) with
C                             I-ML .LE. J .LE. I+MU
C              for all I=1,2,...,NEQ. Thus, ML and MU are the widths
C              of the lower and upper parts of the band, respectively,
C              with the main diagonal being excluded. If you do not
C              indicate that the equation has a banded matrix of partial
C              derivatives, the code works with a full matrix of NEQ**2
C              elements (stored in the conventional way). Computations
C              with banded matrices cost less time and storage than with
C              full matrices if 2*ML+MU .LT. NEQ. If you tell the
C              code that the matrix of partial derivatives has a banded
C              structure and you want to provide subroutine JAC to
C              compute the partial derivatives, then you must be careful
C              to store the elements of the matrix in the special form
C              indicated in the description of JAC.
C          **** Do you want to solve the problem using a full
C               (dense) matrix (and not a special banded
C               structure) ...
C                Yes - Set INFO(6)=0
C                 No - Set INFO(6)=1
C                       and provide the lower (ML) and upper (MU)
C                       bandwidths by setting
C                       IWORK(1)=ML
C                       IWORK(2)=MU ****
C        INFO(7) -- You can specify a maximum (absolute value of)
C              stepsize, so that the code
C              will avoid passing over very
C              large regions.
C          ****  Do you want the code to decide
C                on its own maximum stepsize?
C                Yes - Set INFO(7)=0
C                 No - Set INFO(7)=1
C                      and define HMAX by setting
C                      RWORK(2)=HMAX ****
C        INFO(8) -- Differential/algebraic problems
C              may occaisionally suffer from
C              severe scaling difficulties on the
C              first step. If you know a great deal
C              about the scaling of your problem, you can
C              help to alleviate this problem by
C              specifying an initial stepsize H0.
C          ****  Do you want the code to define
C                its own initial stepsize?
C                Yes - Set INFO(8)=0
C                 No - Set INFO(8)=1
C                      and define H0 by setting
C                      RWORK(3)=H0 ****
C        INFO(9) -- If storage is a severe problem,
C              you can save some locations by
C              restricting the maximum order MAXORD.
C              the default value is 5. for each
C              order decrease below 5, the code
C              requires NEQ fewer locations, however
C              it is likely to be slower. In any
C              case, you must have 1 .LE. MAXORD .LE. 5
C          ****  Do you want the maximum order to
C                default to 5?
C                Yes - Set INFO(9)=0
C                 No - Set INFO(9)=1
C                      and define MAXORD by setting
C                      IWORK(3)=MAXORD ****
C        INFO(10) --If you know that the solutions to your equations
C               will always be nonnegative, it may help to set this
C               parameter. However, it is probably best to
C               try the code without using this option first,
C               and only to use this option if that doesn

'tC               work very well.C           ****  Do you want the code to solve the problem withoutC                 invoking any special nonnegativity constraints?C                  Yes - Set INFO(10)=0C                   No - Set INFO(10)=1CC        INFO(11) --DDASRT normally requires the initial T,C               Y, and YPRIME to be consistent. That is,C               you must have F(T,Y,YPRIME) = 0 at the initialC               time. If you do not know the initialC               derivative precisely, you can let DDASRT tryC               to compute it.C          ****   Are the initial T, Y, YPRIME consistent?C                 Yes - Set INFO(11) = 0C                  No - Set INFO(11) = 1,C                       and set YPRIME to an initial approximationC                       to YPRIME.  (If you have no idea whatC                       YPRIME should be, set it to zero. NoteC                       that the initial Y should be suchC                       that there must exist a YPRIME so thatC                       F(T,Y,YPRIME) = 0.)CC        INFO(12) --Maximum number of steps.C          ****   Do you want to let DDASRT use the default limit forC                 the number of steps?C                 Yes - Set INFO(12) = 0C                  No - Set INFO(12) = 1,C                       and define the maximum number of stepsC                       by setting IWORK(21)=MXSTEPCC   RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOLC               error tolerances to tell the code how accurately youC               want the solution to be computed. They must be definedC               as variables because the code may change them. YouC               have two choices --C                     Both RTOL and ATOL are scalars. (INFO(2)=0)C                     Both RTOL and ATOL are vectors. (INFO(2)=1)C               in either case all components must be non-negative.CC               The tolerances are used by the code in a local errorC               test at each step which requires roughly thatC                     ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOLC               for each vector component.C               (More specifically, a root-mean-square norm is used toC               measure the size of vectors, and the error test uses theC               magnitude of the solution at the beginning of the step.)CC               The true (global) error is the difference between theC               true solution of the initial value problem and theC               computed approximation. Practically all present dayC               codes, including this one, control the local error atC               each step and do not even attempt to control the globalC               error directly.C               Usually, but not always, the true accuracy of theC               computed Y is comparable to the error tolerances. ThisC               code will usually, but not always, deliver a moreC               accurate solution if you reduce the tolerances andC               integrate again. By comparing two such solutions youC               can get a fairly reliable idea of the true error in theC               solution at the bigger tolerances.CC               Setting ATOL=0. results in a pure relative error test onC               that component. Setting RTOL=0. results in a pureC               absolute error test on that component. A mixed testC               with non-zero RTOL and ATOL corresponds roughly to aC               relative error test when the solution component is muchC               bigger than ATOL and to an absolute error test when theC               solution component is smaller than the threshhold ATOL.CC               The code will not attempt to compute a solution at anC               accuracy unreasonable for the machine being used. ItC               will advise you if you ask for too much accuracy andC               inform you as to the maximum accuracy it believesC               possible.CC  RWORK(*) --  Dimension this real work array of length LRW in yourC               calling program.CC  LRW -- Set it to the declared length of the RWORK array.C               You must haveC                    LRW .GE. 50+(MAXORD+4)*NEQ+NEQ**2+3*NGC               for the full (dense) JACOBIAN case (when INFO(6)=0), orC                    LRW .GE. 50+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ+3*NGC               for the banded user-defined JACOBIAN caseC               (when INFO(5)=1 and INFO(6)=1), orC                     LRW .GE. 50+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQC                           +2*(NEQ/(ML+MU+1)+1)+3*NGC               for the banded finite-difference-generated JACOBIAN caseC               (when INFO(5)=0 and INFO(6)=1)CC  IWORK(*) --  Dimension this integer work array of length LIW inC               your calling program.CC  LIW -- Set it to the declared length of the IWORK array.C               you must have LIW .GE. 21+NEQCC  RPAR, IPAR -- These are parameter arrays, of real and integerC               type, respectively. You can use them for communicationC               between your program that calls DDASRT and theC               RES subroutine (and the JAC subroutine). They are notC               altered by DDASRT. If you do not need RPAR or IPAR,C               ignore these parameters by treating them as dummyC               arguments. If you do choose to use them, dimensionC               them in your calling program and in RES (and in JAC)C               as arrays of appropriate length.CC  JAC -- If you have set INFO(5)=0, you can ignore this parameterC               by treating it as a dummy argument. Otherwise, you mustC               provide a subroutine of the formC               JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR)C               to define the matrix of partial derivativesC               PD=DG/DY+CJ*DG/DYPRIMEC               CJ is a scalar which is input to JAC.C               For the given values of T,Y,YPRIME, theC               subroutine must evaluate the non-zero partialC               derivatives for each equation and each solutionC               component, and store these values in theC               matrix PD. The elements of PD are set to zeroC               before each call to JAC so only non-zero elementsC               need to be defined.CC               Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ.C               You must declare the name JAC in anC               EXTERNAL STATEMENT in your program that callsC               DDASRT. You must dimension Y, YPRIME and PDC               in JAC.CC               The way you must store the elements into the PD matrixC               depends on the structure of the matrix which youC               indicated by INFO(6).C               *** INFO(6)=0 -- Full (dense) matrix ***C                   Give PD a first dimension of NEQ.C                   When you evaluate the (non-zero) partial derivativeC                   of equation I with respect to variable J, you mustC                   store it in PD according toC                   PD(I,J) = * DF(I)/DY(J)+CJ*DF(I)/DYPRIME(J)*C               *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MUC                   upper diagonal bands (refer to INFO(6) descriptionC                   of ML and MU) ***C                   Give PD a first dimension of 2*ML+MU+1.C                   when you evaluate the (non-zero) partial derivativeC                   of equation I with respect to variable J, you mustC                   store it in PD according toC                   IROW = I - J + ML + MU + 1C                   PD(IROW,J) = *DF(I)/DY(J)+CJ*DF(I)/DYPRIME(J)*C               RPAR and IPAR are real and integer parameter arraysC               which you can use for communication between your callingC               program and your JACOBIAN subroutine JAC. They are notC               altered by DDASRT. If you do not need RPAR or IPAR,C               ignore these parameters by treating them as dummyC               arguments. If you do choose to use them, dimensionC               them in your calling program and in JAC as arrays ofC               appropriate length.CC  G -- This is the name of the subroutine for defining constraintC               functions, whose roots are desired during theC               integration.  It is to have the formC                   SUBROUTINE G(NEQ,T,Y,NG,GOUT,RPAR,IPAR)C                   DIMENSION Y(NEQ),GOUT(NG),C               where NEQ, T, Y and NG are INPUT, and the array GOUT isC               output.  NEQ, T, and Y have the same meaning as in theC               RES routine, and GOUT is an array of length NG.C               For I=1,...,NG, this routine is to load into GOUT(I)C               the value at (T,Y) of the I-th constraint function G(I).C               DDASRT will find roots of the G(I) of odd multiplicityC               (that is, sign changes) as they occur duringC               the integration.  G must be declared EXTERNAL in theC               calling program.CC               CAUTION..because of numerical errors in the functionsC               G(I) due to roundoff and integration error, DDASRTC               may return false roots, or return the same root at twoC               or more nearly equal values of T.  If such false rootsC               are suspected, the user should consider smaller errorC               tolerances and/or higher precision in the evaluation ofC               the G(I).CC               If a root of some G(I) defines the end of the problem,C               the input to DDASRT should nevertheless allowC               integration to a point slightly past that ROOT, soC               that DDASRT can locate the root by interpolation.CC  NG -- The number of constraint functions G(I).  If there are none,C               set NG = 0, and pass a dummy name for G.CC JROOT -- This is an integer array of length NG.  It is used only forC               output.  On a return where one or more roots have beenC               found, JROOT(I)=1 If G(I) has a root at T,C               or JROOT(I)=0 if not.CCCC  OPTIONALLY REPLACEABLE NORM ROUTINE:C  DDASRT uses a weighted norm DDANRM to measure the sizeC  of vectors such as the estimated error in each step.C  A FUNCTION subprogramC    DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR)C    DIMENSION V(NEQ),WT(NEQ)C  is used to define this norm. Here, V is the vectorC  whose norm is to be computed, and WT is a vector ofC  weights.  A DDANRM routine has been included with DDASRTC  which computes the weighted root-mean-square normC  given byC    DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2)C  this norm is suitable for most problems. In someC  special cases, it may be more convenient and/orC  efficient to define your own norm by writing a functionC  subprogram to be called instead of DDANRM. This shouldC  ,however, be attempted only after careful thought andC  consideration.CCC------OUTPUT-AFTER ANY RETURN FROM DDASRT----CC  The principal aim of the code is to return a computed solution atC  TOUT, although it is also possible to obtain intermediate resultsC  along the way. To find out whether the code achieved its goalC  or if the integration process was interrupted before the task wasC  completed, you must check the IDID parameter.CCC   T -- The solution was successfully advanced to theC               output value of T.CC   Y(*) -- Contains the computed solution approximation at T.CC   YPRIME(*) -- Contains the computed derivativeC               approximation at T.CC   IDID -- Reports what the code did.CC                     *** Task completed ***C                Reported by positive values of IDIDCC           IDID = 1 -- A step was successfully taken in theC                   intermediate-output mode. The code has notC                   yet reached TOUT.CC           IDID = 2 -- The integration to TSTOP was successfullyC                   completed (T=TSTOP) by stepping exactly to TSTOP.CC           IDID = 3 -- The integration to TOUT was successfullyC                   completed (T=TOUT) by stepping past TOUT.C                   Y(*) is obtained by interpolation.C                   YPRIME(*) is obtained by interpolation.CC           IDID = 4 -- The integration was successfully completedC                   by finding one or more roots of G at T.CC                    *** Task interrupted ***C                Reported by negative values of IDIDCC           IDID = -1 -- A large amount of work has been expended.C                   (About INFO(12) steps)CC           IDID = -2 -- The error tolerances are too stringent.CC           IDID = -3 -- The local error test cannot be satisfiedC                   because you specified a zero component in ATOLC                   and the corresponding computed solutionC                   component is zero. Thus, a pure relative errorC                   test is impossible for this component.CC           IDID = -6 -- DDASRT had repeated error testC                   failures on the last attempted step.CC           IDID = -7 -- The corrector could not converge.CC           IDID = -8 -- The matrix of partial derivativesC                   is singular.CC           IDID = -9 -- The corrector could not converge.C                   there were repeated error test failuresC                   in this step.CC           IDID =-10 -- The corrector could not convergeC                   because IRES was equal to minus one.CC           IDID =-11 -- IRES equal to -2 was encounteredC                   and control is being returned to theC                   calling program.CC           IDID =-12 -- DDASRT failed to compute the initialC                   YPRIME.CCCC           IDID = -13,..,-32 -- Not applicable for this codeCC                    *** Task terminated ***C                Reported by the value of IDID=-33CC           IDID = -33 -- The code has encountered trouble from whichC                   it cannot recover. A message is printedC                   explaining the trouble and control is returnedC                   to the calling program. For example, this occursC                   when invalid input is detected.CC   RTOL, ATOL -- These quantities remain unchanged except whenC               IDID = -2. In this case, the error tolerances have beenC               increased by the code to values which are estimated toC               be appropriate for continuing the integration. However,C               the reported solution at T was obtained using the inputC               values of RTOL and ATOL.CC   RWORK, IWORK -- Contain information which is usually of noC               interest to the user but necessary for subsequent calls.C               However, you may find use forCC               RWORK(3)--Which contains the step size H to beC                       attempted on the next step.CC               RWORK(4)--Which contains the current value of theC                       independent variable, i.e., the farthest pointC                       integration has reached. This will be differentC                       from T only when interpolation has beenC                       performed (IDID=3).CC               RWORK(7)--Which contains the stepsize usedC                       on the last successful step.CC               IWORK(7)--Which contains the order of the method toC                       be attempted on the next step.CC               IWORK(8)--Which contains the order of the method usedC                       on the last step.CC               IWORK(11)--Which contains the number of steps taken soC                        far.CC               IWORK(12)--Which contains the number of calls to RESC                        so far.CC               IWORK(13)--Which contains the number of evaluations ofC                        the matrix of partial derivatives needed soC                        far.CC               IWORK(14)--Which contains the total numberC                        of error test failures so far.CC               IWORK(15)--Which contains the total numberC                        of convergence test failures so far.C                        (includes singular iteration matrixC                        failures.)CC               IWORK(16)--Which contains the total number of callsC                        to the constraint function g so farCCCC   INPUT -- What to do to continue the integrationC            (calls after the first)                **CC     This code is organized so that subsequent calls to continue theC     integration involve little (if any) additional effort on yourC     part. You must monitor the IDID parameter in order to determineC     what to do next.CC     Recalling that the principal task of the code is to integrateC     from T to TOUT (the interval mode), usually all you will needC     to do is specify a new TOUT upon reaching the current TOUT.CC     Do not alter any quantity not specifically permitted below,C     in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*)C     or the differential equation in subroutine RES. Any suchC     alteration constitutes a new problem and must be treated as such,C     i.e., you must start afresh.CC     You cannot change from vector to scalar error control or viceC     versa (INFO(2)), but you can change the size of the entries ofC     RTOL, ATOL. Increasing a tolerance makes the equation easierC     to integrate. Decreasing a tolerance will make the equationC     harder to integrate and should generally be avoided.CC     You can switch from the intermediate-output mode to theC     interval mode (INFO(3)) or vice versa at any time.CC     If it has been necessary to prevent the integration from goingC     past a point TSTOP (INFO(4), RWORK(1)), keep in mind that theC     code will not integrate to any TOUT beyond the currentlyC     specified TSTOP. Once TSTOP has been reached you must changeC     the value of TSTOP or set INFO(4)=0. You may change INFO(4)C     or TSTOP at any time but you must supply the value of TSTOP inC     RWORK(1) whenever you set INFO(4)=1.CC     Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2)C     unless you are going to restart the code.CC                    *** Following a completed task ***C     IfC     IDID = 1, call the code again to continue the integrationC                  another step in the direction of TOUT.CC     IDID = 2 or 3, define a new TOUT and call the code again.C                  TOUT must be different from T. You cannot changeC                  the direction of integration without restarting.CC     IDID = 4, call the code again to continue the integrationC                  another step in the direction of TOUT.  You mayC                  change the functions in G after a return with IDID=4,C                  but the number of constraint functions NG must remainC                  the same.  If you wish to changeC                  the functions in RES or in G, then youC                  must restart the code.CC                    *** Following an interrupted task ***C                  To show the code that you realize the task wasC                  interrupted and that you want to continue, youC                  must take appropriate action and set INFO(1) = 1C     IfC     IDID = -1, The code has reached the step iteration.C                  If you want to continue, set INFO(1) = 1 andC                  call the code again.  See also INFO(12).CC     IDID = -2, The error tolerances RTOL, ATOL have beenC                  increased to values the code estimates appropriateC                  for continuing. You may want to change themC                  yourself. If you are sure you want to continueC                  with relaxed error tolerances, set INFO(1)=1 andC                  call the code again.CC     IDID = -3, A solution component is zero and you set theC                  corresponding component of ATOL to zero. If youC                  are sure you want to continue, you must firstC                  alter the error criterion to use positive valuesC                  for those components of ATOL corresponding to zeroC                  solution components, then set INFO(1)=1 and callC                  the code again.CC     IDID = -4,-5  --- Cannot occur with this code.CC     IDID = -6, Repeated error test failures occurred on theC                  last attempted step in DDASRT. A singularity in theC                  solution may be present. If you are absolutelyC                  certain you want to continue, you should restartC                  the integration. (Provide initial values of Y andC                  YPRIME which are consistent)CC     IDID = -7, Repeated convergence test failures occurredC                  on the last attempted step in DDASRT. An inaccurateC                  or ill-conditioned JACOBIAN may be the problem. IfC                  you are absolutely certain you want to continue, youC                  should restart the integration.CC     IDID = -8, The matrix of partial derivatives is singular.C                  Some of your equations may be redundant.C                  DDASRT cannot solve the problem as stated.C                  It is possible that the redundant equationsC                  could be removed, and then DDASRT couldC                  solve the problem. It is also possibleC                  that a solution to your problem eitherC                  does not exist or is not unique.CC     IDID = -9, DDASRT had multiple convergence testC                  failures, preceeded by multiple errorC                  test failures, on the last attempted step.C                  It is possible that your problemC                  is ill-posed, and cannot be solvedC                  using this code. Or, there may be aC                  discontinuity or a singularity in theC                  solution. If you are absolutely certainC                  you want to continue, you should restartC                  the integration.CC    IDID =-10, DDASRT had multiple convergence test failuresC                  because IRES was equal to minus one.C                  If you are absolutely certain you wantC                  to continue, you should restart theC                  integration.CC    IDID =-11, IRES=-2 was encountered, and control is beingC                  returned to the calling program.CC    IDID =-12, DDASRT failed to compute the initial YPRIME.C               This could happen because the initialC               approximation to YPRIME was not very good, orC               if a YPRIME consistent with the initial YC               does not exist. The problem could also be causedC               by an inaccurate or singular iteration matrix.CCCC     IDID = -13,..,-32 --- Cannot occur with this code.CC                       *** Following a terminated task ***C     If IDID= -33, you cannot continue the solution of thisC                  problem. An attempt to do so will result in yourC                  run being terminated.CC  ---------------------------------------------------------------------CC***REFERENCEC      K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical C      Solution of Initial-Value Problems in Differential-AlgebraicC      Equations, Elsevier, New York, 1989.CC***ROUTINES CALLED  DDASTP,DDAINI,DDANRM,DDAWTS,DDATRP,DRCHEK,DROOTS,C                    XERRWD,D1MACHC***END PROLOGUE  DDASRTCC**EndC      IMPLICIT DOUBLE PRECISION(A-H,O-Z)      LOGICAL DONE      EXTERNAL RES, JAC, G      DIMENSION Y(*),YPRIME(*)      DIMENSION INFO(15)      DIMENSION RWORK(*),IWORK(*)      DIMENSION RTOL(*),ATOL(*)      DIMENSION RPAR(*),IPAR(*)      CHARACTER MSG*80CC     SET POINTERS INTO IWORK      PARAMETER (LML=1, LMU=2, LMXORD=3, LMTYPE=4, LNST=11,     *  LNRE=12, LNJE=13, LETF=14, LCTF=15, LNGE=16, LNPD=17,     *  LIRFND=18, LMXSTP=21, LIPVT=22, LJCALC=5, LPHASE=6, LK=7,     *  LKOLD=8, LNS=9, LNSTL=10, LIWM=1)CC     SET RELATIVE OFFSET INTO RWORK      PARAMETER (NPD=1)CC     SET POINTERS INTO RWORK      PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4,     *  LCJ=5, LCJOLD=6, LHOLD=7, LS=8, LROUND=9,     *  LALPHA=11, LBETA=17, LGAMMA=23,     *  LPSI=29, LSIGMA=35, LT0=41, LTLAST=42, LALPHR=43, LX2=44,     *  LDELTA=51)CC***FIRST EXECUTABLE STATEMENT  DDASRT      IF(INFO(1).NE.0)GO TO 100CC-----------------------------------------------------------------------C     THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY.C     IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS.C-----------------------------------------------------------------------CC     FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFOC     ARE EITHER ZERO OR ONE.      DO 10 I=2,12         IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 70110       CONTINUEC      IF(NEQ.LE.0)GO TO 702CC     CHECK AND COMPUTE MAXIMUM ORDER      MXORD=5      IF(INFO(9).EQ.0)GO TO 20         MXORD=IWORK(LMXORD)         IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 70320       IWORK(LMXORD)=MXORDCC     COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU.      IF(INFO(6).NE.0)GO TO 40         LENPD=NEQ**2         LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD+3*NG         IF(INFO(5).NE.0)GO TO 30            IWORK(LMTYPE)=2            GO TO 6030          IWORK(LMTYPE)=1            GO TO 6040    IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717      IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718      LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ      IF(INFO(5).NE.0)GO TO 50         IWORK(LMTYPE)=5         MBAND=IWORK(LML)+IWORK(LMU)+1         MSAVE=(NEQ/MBAND)+1         LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE+3*NG         GO TO 6050       IWORK(LMTYPE)=4         LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD+3*NGCC     CHECK LENGTHS OF RWORK AND IWORK60    LENIW=21+NEQ      IWORK(LNPD)=LENPD      IF(LRW.LT.LENRW)GO TO 704      IF(LIW.LT.LENIW)GO TO 705CC     CHECK TO SEE THAT TOUT IS DIFFERENT FROM TC     Also check to see that NG is larger than 0.      IF(TOUT .EQ. T)GO TO 719      IF(NG .LT. 0) GO TO 730CC     CHECK HMAX      IF(INFO(7).EQ.0)GO TO 70         HMAX=RWORK(LHMAX)         IF(HMAX.LE.0.0D0)GO TO 71070    CONTINUECC     CHECK AND COMPUTE MAXIMUM STEPS      MXSTP=500      IF(INFO(12).EQ.0)GO TO 80        MXSTP=IWORK(LMXSTP)        IF(MXSTP.LT.0)GO TO 71680      IWORK(LMXSTP)=MXSTPCC     INITIALIZE COUNTERS      IWORK(LNST)=0      IWORK(LNRE)=0      IWORK(LNJE)=0      IWORK(LNGE)=0C      IWORK(LNSTL)=0      IDID=1      GO TO 200CC-----------------------------------------------------------------------C     THIS BLOCK IS FOR CONTINUATION CALLSC     ONLY. HERE WE CHECK INFO(1),AND IF THEC     LAST STEP WAS INTERRUPTED WE CHECK WHETHERC     APPROPRIATE ACTION WAS TAKEN.C-----------------------------------------------------------------------C100   CONTINUE      IF(INFO(1).EQ.1)GO TO 110      IF(INFO(1).NE.-1)GO TO 701C     IF WE ARE HERE, THE LAST STEP WAS INTERRUPTEDC     BY AN ERROR CONDITION FROM DDASTP,ANDC     APPROPRIATE ACTION WAS NOT TAKEN. THISC     IS A FATAL ERROR.      MSG = 'DASRT--  THE LAST STEP TERMINATED WITH A NEGATIVE

'      CALL XERRWD(MSG,49,201,0,0,0,0,0,0.0D0,0.0D0)      MSG = 'DASRT--  VALUE (=I1) OF IDID AND NO APPROPRIATE

'      CALL XERRWD(MSG,47,202,0,1,IDID,0,0,0.0D0,0.0D0)      MSG = 'DASRT--  ACTION WAS TAKEN. RUN TERMINATED

'      CALL XERRWD(MSG,41,203,1,0,0,0,0,0.0D0,0.0D0)      RETURN110   CONTINUE      IWORK(LNSTL)=IWORK(LNST)CC-----------------------------------------------------------------------C     THIS BLOCK IS EXECUTED ON ALL CALLS.C     THE ERROR TOLERANCE PARAMETERS AREC     CHECKED, AND THE WORK ARRAY POINTERSC     ARE SET.C-----------------------------------------------------------------------C200   CONTINUEC     CHECK RTOL,ATOL      NZFLG=0      RTOLI=RTOL(1)      ATOLI=ATOL(1)      DO 210 I=1,NEQ         IF(INFO(2).EQ.1)RTOLI=RTOL(I)         IF(INFO(2).EQ.1)ATOLI=ATOL(I)         IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1         IF(RTOLI.LT.0.0D0)GO TO 706         IF(ATOLI.LT.0.0D0)GO TO 707210      CONTINUE      IF(NZFLG.EQ.0)GO TO 708CC     SET UP RWORK STORAGE.IWORK STORAGE IS FIXEDC     IN DATA STATEMENT.      LG0=LDELTA+NEQ      LG1=LG0+NG      LGX=LG1+NG      LE=LGX+NG      LWT=LE+NEQ      LPHI=LWT+NEQ      LPD=LPHI+(IWORK(LMXORD)+1)*NEQ      LWM=LPD      NTEMP=NPD+IWORK(LNPD)      IF(INFO(1).EQ.1)GO TO 400CC-----------------------------------------------------------------------C     THIS BLOCK IS EXECUTED ON THE INITIAL CALLC     ONLY. SET THE INITIAL STEP SIZE, ANDC     THE ERROR WEIGHT VECTOR, AND PHI.C     COMPUTE INITIAL YPRIME, IF NECESSARY.C-----------------------------------------------------------------------C300   CONTINUE      TN=T      IDID=1CC     SET ERROR WEIGHT VECTOR WT      CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)      DO 305 I = 1,NEQ         IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713305      CONTINUECC     COMPUTE UNIT ROUNDOFF AND HMIN      UROUND = D1MACH(4)      RWORK(LROUND) = UROUND      HMIN = 4.0D0*UROUND*DMAX1(DABS(T),DABS(TOUT))CC     CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH      TDIST = DABS(TOUT - T)      IF(TDIST .LT. HMIN) GO TO 714CC     CHECK H0, IF THIS WAS INPUT      IF (INFO(8) .EQ. 0) GO TO 310         HO = RWORK(LH)         IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711         IF (HO .EQ. 0.0D0) GO TO 712         GO TO 320310    CONTINUECC     COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHERC     DDASTP OR DDAINI, DEPENDING ON INFO(11)      HO = 0.001D0*TDIST      YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR)      IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM      HO = DSIGN(HO,TOUT-T)C     ADJUST HO IF NECESSARY TO MEET HMAX BOUND320   IF (INFO(7) .EQ. 0) GO TO 330         RH = DABS(HO)/RWORK(LHMAX)         IF (RH .GT. 1.0D0) HO = HO/RHC     COMPUTE TSTOP, IF APPLICABLE330   IF (INFO(4) .EQ. 0) GO TO 340         TSTOP = RWORK(LTSTOP)         IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715         IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T         IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709CC     COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE340   IF (INFO(11) .EQ. 0) GO TO 350      CALL DDAINI(TN,Y,YPRIME,NEQ,     *  RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR,     *  RWORK(LPHI),RWORK(LDELTA),RWORK(LE),     *  RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND),     *  INFO(10),NTEMP)      IF (IDID .LT. 0) GO TO 390CC     LOAD H WITH H0.  STORE H IN RWORK(LH)350   H = HO      RWORK(LH) = HCC     LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2)360   ITEMP = LPHI + NEQ      DO 370 I = 1,NEQ         RWORK(LPHI + I - 1) = Y(I)370      RWORK(ITEMP + I - 1) = H*YPRIME(I)CC     INITIALIZE T0 IN RWORK AND CHECK FOR A ZERO OF G NEAR THEC     INITIAL T.C      RWORK(LT0) = T      IWORK(LIRFND) = 0      RWORK(LPSI)=H      RWORK(LPSI+1)=2.0D0*H      IWORK(LKOLD)=1      IF(NG .EQ. 0) GO TO 390      CALL DRCHEK(1,G,NG,NEQ,T,TOUT,Y,RWORK(LE),RWORK(LPHI),     *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),     *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),     *  RWORK,IWORK,RPAR,IPAR)      IF(IRT .NE. 0) GO TO 732CC     Check for a root in the interval (T0,TN], unless DDASRTC     did not have to initialize YPRIME.C      IF(NG .EQ. 0 .OR. INFO(11) .EQ. 0) GO TO 390      CALL DRCHEK(3,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI),     *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),     *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),     *  RWORK,IWORK,RPAR,IPAR)      IF(IRT .NE. 1) GO TO 390      IWORK(LIRFND) = 1      IDID = 4      T = RWORK(LT0)      GO TO 580C390   GO TO 500CC-------------------------------------------------------C     THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITSC     PURPOSE IS TO CHECK STOP CONDITIONS BEFOREC     TAKING A STEP.C     ADJUST H IF NECESSARY TO MEET HMAX BOUNDC-------------------------------------------------------C400   CONTINUE      UROUND=RWORK(LROUND)      DONE = .FALSE.      TN=RWORK(LTN)      H=RWORK(LH)      IF(NG .EQ. 0) GO TO 405CC     Check for a zero of G near TN.C      CALL DRCHEK(2,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI),     *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),     *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),     *  RWORK,IWORK,RPAR,IPAR)      IF(IRT .NE. 1) GO TO 405      IWORK(LIRFND) = 1      IDID = 4      T = RWORK(LT0)      DONE = .TRUE.      GO TO 490C405   CONTINUE      IF(INFO(7) .EQ. 0) GO TO 410         RH = DABS(H)/RWORK(LHMAX)         IF(RH .GT. 1.0D0) H = H/RH410   CONTINUE      IF(T .EQ. TOUT) GO TO 719      IF((T - TOUT)*H .GT. 0.0D0) GO TO 711      IF(INFO(4) .EQ. 1) GO TO 430      IF(INFO(3) .EQ. 1) GO TO 420      IF((TN-TOUT)*H.LT.0.0D0)GO TO 490      CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),     *  RWORK(LPHI),RWORK(LPSI))      T=TOUT      IDID = 3      DONE = .TRUE.      GO TO 490420   IF((TN-T)*H .LE. 0.0D0) GO TO 490      IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425      CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),     *  RWORK(LPHI),RWORK(LPSI))      T = TN      IDID = 1      DONE = .TRUE.      GO TO 490425   CONTINUE      CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),     *  RWORK(LPHI),RWORK(LPSI))      T = TOUT      IDID = 3      DONE = .TRUE.      GO TO 490430   IF(INFO(3) .EQ. 1) GO TO 440      TSTOP=RWORK(LTSTOP)      IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715      IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709      IF((TN-TOUT)*H.LT.0.0D0)GO TO 450      CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),     *   RWORK(LPHI),RWORK(LPSI))      T=TOUT      IDID = 3      DONE = .TRUE.      GO TO 490440   TSTOP = RWORK(LTSTOP)      IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715      IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709      IF((TN-T)*H .LE. 0.0D0) GO TO 450      IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445      CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),     *  RWORK(LPHI),RWORK(LPSI))      T = TN      IDID = 1      DONE = .TRUE.      GO TO 490445   CONTINUE      CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),     *  RWORK(LPHI),RWORK(LPSI))      T = TOUT      IDID = 3      DONE = .TRUE.      GO TO 490450   CONTINUEC     CHECK WHETHER WE ARE WITH IN ROUNDOFF OF TSTOP      IF(DABS(TN-TSTOP).GT.100.0D0*UROUND*     *   (DABS(TN)+DABS(H)))GO TO 460      CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD),     *  RWORK(LPHI),RWORK(LPSI))      IDID=2      T=TSTOP      DONE = .TRUE.      GO TO 490460   TNEXT=TN+H      IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490      H=TSTOP-TN      RWORK(LH)=HC490   IF (DONE) GO TO 590CC-------------------------------------------------------C     THE NEXT BLOCK CONTAINS THE CALL TO THEC     ONE-STEP INTEGRATOR DDASTP.C     THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS.C     CHECK FOR TOO MANY STEPS.C     UPDATE WT.C     CHECK FOR TOO MUCH ACCURACY REQUESTED.C     COMPUTE MINIMUM STEPSIZE.C-------------------------------------------------------C500   CONTINUEC     CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME      IF (IDID .EQ. -12) GO TO 527CC     CHECK FOR TOO MANY STEPS      IF((IWORK(LNST)-IWORK(LNSTL)).LT.IWORK(LMXSTP))     *   GO TO 510           IDID=-1           GO TO 527CC     UPDATE WT510   CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),     *  RWORK(LWT),RPAR,IPAR)      DO 520 I=1,NEQ         IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520           IDID=-3           GO TO 527520   CONTINUECC     TEST FOR TOO MUCH ACCURACY REQUESTED.      R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*     *   100.0D0*UROUND      IF(R.LE.1.0D0)GO TO 525C     MULTIPLY RTOL AND ATOL BY R AND RETURN      IF(INFO(2).EQ.1)GO TO 523           RTOL(1)=R*RTOL(1)           ATOL(1)=R*ATOL(1)           IDID=-2           GO TO 527523   DO 524 I=1,NEQ           RTOL(I)=R*RTOL(I)524        ATOL(I)=R*ATOL(I)      IDID=-2      GO TO 527525   CONTINUECC     COMPUTE MINIMUM STEPSIZE      HMIN=4.0D0*UROUND*DMAX1(DABS(TN),DABS(TOUT))CC     TEST H VS. HMAX      IF (INFO(7) .EQ. 0) GO TO 526         RH = ABS(H)/RWORK(LHMAX)         IF (RH .GT. 1.0D0) H = H/RH526   CONTINUE     C      CALL DDASTP(TN,Y,YPRIME,NEQ,     *   RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR,     *   RWORK(LPHI),RWORK(LDELTA),RWORK(LE),     *   RWORK(LWM),IWORK(LIWM),     *   RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),     *   RWORK(LPSI),RWORK(LSIGMA),     *   RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),     *   RWORK(LS),HMIN,RWORK(LROUND),     *   IWORK(LPHASE),IWORK(LJCALC),IWORK(LK),     *   IWORK(LKOLD),IWORK(LNS),INFO(10),NTEMP)527   IF(IDID.LT.0)GO TO 600CC--------------------------------------------------------C     THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURNC     FROM DDASTP (IDID=1).  TEST FOR STOP CONDITIONS.C--------------------------------------------------------C      IF(NG .EQ. 0) GO TO 529CC     Check for a zero of G near TN.C      CALL DRCHEK(3,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI),     *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),     *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),     *  RWORK,IWORK,RPAR,IPAR)      IF(IRT .NE. 1) GO TO 529      IWORK(LIRFND) = 1      IDID = 4      T = RWORK(LT0)      GO TO 580C529   CONTINUE      IF(INFO(4).NE.0)GO TO 540           IF(INFO(3).NE.0)GO TO 530             IF((TN-TOUT)*H.LT.0.0D0)GO TO 500             CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,     *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))             IDID=3             T=TOUT             GO TO 580530          IF((TN-TOUT)*H.GE.0.0D0)GO TO 535             T=TN             IDID=1             GO TO 580535          CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,     *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))             IDID=3             T=TOUT             GO TO 580540   IF(INFO(3).NE.0)GO TO 550      IF((TN-TOUT)*H.LT.0.0D0)GO TO 542         CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,     *     IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))         T=TOUT         IDID=3         GO TO 580542   IF(DABS(TN-TSTOP).LE.100.0D0*UROUND*     *   (DABS(TN)+DABS(H)))GO TO 545      TNEXT=TN+H      IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500      H=TSTOP-TN      GO TO 500545   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,     *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))      IDID=2      T=TSTOP      GO TO 580550   IF((TN-TOUT)*H.GE.0.0D0)GO TO 555      IF(DABS(TN-TSTOP).LE.100.0D0*UROUND*(DABS(TN)+DABS(H)))GO TO 552      T=TN      IDID=1      GO TO 580552   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,     *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))      IDID=2      T=TSTOP      GO TO 580555   CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,     *   IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))      T=TOUT      IDID=3580   CONTINUECC--------------------------------------------------------C     ALL SUCCESSFUL RETURNS FROM DDASRT ARE MADE FROMC     THIS BLOCK.C--------------------------------------------------------C590   CONTINUE      RWORK(LTN)=TN      RWORK(LH)=H      RWORK(LTLAST) = T      RETURNCC-----------------------------------------------------------------------C     THIS BLOCK HANDLES ALL UNSUCCESSFULC     RETURNS OTHER THAN FOR ILLEGAL INPUT.C-----------------------------------------------------------------------C600   CONTINUE      ITEMP=-IDID      GO TO (610,620,630,690,690,640,650,660,670,675,     *  680,685), ITEMPCC     THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFOREC     REACHING TOUT610   MSG = 'DASRT--  AT CURRENT T (=R1)  500 STEPS




'      CALL XERRWD(MSG,48,621,0,0,0,0,0,0.0D0,0.0D0)      MSG = 'DASRT--  WERE INCREASED TO APPROPRIATE VALUES

'      CALL XERRWD(MSG,45,622,0,0,0,0,0,0.0D0,0.0D0)C      GO TO 690C     WT(I) .LE. 0.0D0 FOR SOME I (NOT AT START OF PROBLEM)630   MSG = 'DASRT--  AT T (=R1) SOME ELEMENT OF WT

'      CALL XERRWD(MSG,38,630,0,0,0,0,1,TN,0.0D0)      MSG = 'DASRT--  HAS BECOME .LE. 0.0

'      CALL XERRWD(MSG,28,631,0,0,0,0,0,0.0D0,0.0D0)      GO TO 690CC     ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN640   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2) THE




'      CALL XERRWD(MSG,48,651,0,0,0,0,0,0.0D0,0.0D0)      MSG = 'DASRT--  OR WITH ABS(H)=HMIN

'      CALL XERRWD(MSG,28,652,0,0,0,0,0,0.0D0,0.0D0)      GO TO 690CC     THE ITERATION MATRIX IS SINGULAR660   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2) THE




'      CALL XERRWD(MSG,49,671,0,0,0,0,0,0.0D0,0.0D0)      MSG = 'DASRT--  ERROR TEST FAILED REPEATEDLY.

'      CALL XERRWD(MSG,38,672,0,0,0,0,0,0.0D0,0.0D0)      GO TO 690CC     CORRECTOR FAILURE BECAUSE IRES = -1675   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2) THE


'      CALL XERRWD(MSG,45,676,0,0,0,0,0,0.0D0,0.0D0)      MSG = 'DASRT--  IRES WAS EQUAL TO MINUS ONE

'      CALL XERRWD(MSG,36,677,0,0,0,0,0,0.0D0,0.0D0)      GO TO 690CC     FAILURE BECAUSE IRES = -2680   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2)

'      CALL XERRWD(MSG,40,680,0,0,0,0,2,TN,H)      MSG = 'DASRT--  IRES WAS EQUAL TO MINUS TWO

'      CALL XERRWD(MSG,36,681,0,0,0,0,0,0.0D0,0.0D0)      GO TO 690CC     FAILED TO COMPUTE INITIAL YPRIME685   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2) THE


'      CALL XERRWD(MSG,45,686,0,0,0,0,0,0.0D0,0.0D0)      GO TO 690690   CONTINUE      INFO(1)=-1      T=TN      RWORK(LTN)=TN      RWORK(LH)=H      RETURNC-----------------------------------------------------------------------C     THIS BLOCK HANDLES ALL ERROR RETURNS DUEC     TO ILLEGAL INPUT, AS DETECTED BEFORE CALLINGC     DDASTP. FIRST THE ERROR MESSAGE ROUTINE ISC     CALLED. IF THIS HAPPENS TWICE INC     SUCCESSION, EXECUTION IS TERMINATEDCC-----------------------------------------------------------------------701   MSG = 'DASRT--  SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE

'      CALL XERRWD(MSG,55,1,0,0,0,0,0,0.0D0,0.0D0)      GO TO 750702   MSG = 'DASRT--  NEQ (=I1) .LE. 0

'      CALL XERRWD(MSG,25,2,0,1,NEQ,0,0,0.0D0,0.0D0)      GO TO 750703   MSG = 'DASRT--  MAXORD (=I1) NOT IN RANGE

'      CALL XERRWD(MSG,34,3,0,1,MXORD,0,0,0.0D0,0.0D0)      GO TO 750704   MSG='DASRT--  RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)

'      CALL XERRWD(MSG,60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0)      GO TO 750705   MSG='DASRT--  IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)

'      CALL XERRWD(MSG,60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0)      GO TO 750706   MSG = 'DASRT--  SOME ELEMENT OF RTOL IS .LT. 0

'      CALL XERRWD(MSG,39,6,0,0,0,0,0,0.0D0,0.0D0)      GO TO 750707   MSG = 'DASRT--  SOME ELEMENT OF ATOL IS .LT. 0

'      CALL XERRWD(MSG,39,7,0,0,0,0,0,0.0D0,0.0D0)      GO TO 750708   MSG = 'DASRT--  ALL ELEMENTS OF RTOL AND ATOL ARE ZERO

'      CALL XERRWD(MSG,47,8,0,0,0,0,0,0.0D0,0.0D0)      GO TO 750709   MSG='DASRT--  INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)

'      CALL XERRWD(MSG,54,9,0,0,0,0,2,TSTOP,TOUT)      GO TO 750710   MSG = 'DASRT--  HMAX (=R1) .LT. 0.0

'      CALL XERRWD(MSG,28,10,0,0,0,0,1,HMAX,0.0D0)      GO TO 750711   MSG = 'DASRT--  TOUT (=R1) BEHIND T (=R2)

'      CALL XERRWD(MSG,34,11,0,0,0,0,2,TOUT,T)      GO TO 750712   MSG = 'DASRT--  INFO(8)=1 AND H0=0.0

'      CALL XERRWD(MSG,29,12,0,0,0,0,0,0.0D0,0.0D0)      GO TO 750713   MSG = 'DASRT--  SOME ELEMENT OF WT IS .LE. 0.0

'      CALL XERRWD(MSG,39,13,0,0,0,0,0,0.0D0,0.0D0)      GO TO 750714   MSG='DASRT-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION

'      CALL XERRWD(MSG,60,14,0,0,0,0,2,TOUT,T)      GO TO 750715   MSG = 'DASRT--  INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)

'      CALL XERRWD(MSG,49,15,0,0,0,0,2,TSTOP,T)      GO TO 750716   MSG = 'DASRT--  INFO(12)=1 AND MXSTP (=I1) .LT. 0

'      CALL XERRWD(MSG,42,16,0,1,IWORK(LMXSTP),0,0,0.0D0,0.0D0)      GO TO 750717   MSG = 'DASRT--  ML (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ

'      CALL XERRWD(MSG,52,17,0,1,IWORK(LML),0,0,0.0D0,0.0D0)      GO TO 750718   MSG = 'DASRT--  MU (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ

'      CALL XERRWD(MSG,52,18,0,1,IWORK(LMU),0,0,0.0D0,0.0D0)      GO TO 750719   MSG = 'DASRT--  TOUT (=R1) IS EQUAL TO T (=R2)

'      CALL XERRWD(MSG,39,19,0,0,0,0,2,TOUT,T)      GO TO 750730   MSG = 'DASRT--  NG (=I1) .LT. 0

'      CALL XERRWD(MSG,24,30,1,1,NG,0,0,0.0D0,0.0D0)      GO TO 750732   MSG = 'DASRT--  ONE OR MORE COMPONENTS OF G HAS A ROOT

'      CALL XERRWD(MSG,47,32,1,0,0,0,0,0.0D0,0.0D0)      MSG = '         TOO NEAR TO THE INITIAL POINT

'      CALL XERRWD(MSG,38,32,1,0,0,0,0,0.0D0,0.0D0)750   IF(INFO(1).EQ.-1) GO TO 760      INFO(1)=-1      IDID=-33      RETURN760   MSG = 'DASRT--  REPEATED OCCURRENCES OF ILLEGAL INPUT

'      CALL XERRWD(MSG,46,801,0,0,0,0,0,0.0D0,0.0D0)770   MSG = 'DASRT--  RUN TERMINATED. APPARENT INFINITE LOOP

Generated by  Doxygen 1.6.0   Back to index