research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
cantera/ext/math/ddaspk.f

6612 lines
245 KiB
Fortran
Raw Blame History

SUBROUTINE DDASPK (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
* IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL)
C
C***BEGIN PROLOGUE DDASPK
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 910624 (Added HMAX test at 525 in main driver.)
C***REVISION DATE 920929 (CJ in RES call, RES counter fix.)
C***REVISION DATE 921215 (Warnings on poor iteration performance)
C***REVISION DATE 921216 (NRMAX as optional input)
C***REVISION DATE 930315 (Name change: DDINI to DDINIT)
C***REVISION DATE 940822 (Replaced initial condition calculation)
C***REVISION DATE 941101 (Added linesearch in I.C. calculations)
C***REVISION DATE 941220 (Misc. corrections throughout)
C***REVISION DATE 950125 (Added DINVWT routine)
C***REVISION DATE 950714 (Misc. corrections throughout)
C***REVISION DATE 950802 (Default NRMAX = 5, based on tests.)
C***REVISION DATE 950808 (Optional error test added.)
C***REVISION DATE 950814 (Added I.C. constraints and INFO(14))
C***REVISION DATE 950828 (Various minor corrections.)
C***REVISION DATE 951006 (Corrected WT scaling in DFNRMK.)
C***REVISION DATE 951030 (Corrected history update at end of DDASTP.)
C***REVISION DATE 960129 (Corrected RL bug in DLINSD, DLINSK.)
C***REVISION DATE 960301 (Added NONNEG to SAVE statement.)
C***REVISION DATE 000512 (Removed copyright notices.)
C***REVISION DATE 000622 (Corrected LWM value using NCPHI.)
C***REVISION DATE 000628 (Corrected I.C. stopping tests when index = 0.)
C***REVISION DATE 000628 (Fixed alpha test in I.C. calc., Krylov case.)
C***REVISION DATE 000628 (Improved restart in I.C. calc., Krylov case.)
C***REVISION DATE 000628 (Minor corrections throughout.)
C***REVISION DATE 000711 (Fixed Newton convergence test in DNSD, DNSK.)
C***REVISION DATE 000712 (Fixed tests on TN - TOUT below 420 and 440.)
C***CATEGORY NO. I1A2
C***KEYWORDS DIFFERENTIAL/ALGEBRAIC, BACKWARD DIFFERENTIATION FORMULAS,
C IMPLICIT DIFFERENTIAL SYSTEMS, KRYLOV ITERATION
C***AUTHORS Linda R. Petzold, Peter N. Brown, Alan C. Hindmarsh, and
C Clement W. Ulrich
C Center for Computational Sciences & Engineering, L-316
C Lawrence Livermore National Laboratory
C P.O. Box 808,
C Livermore, CA 94551
C***PURPOSE This code solves a system of differential/algebraic
C equations of the form
C G(t,y,y') = 0 ,
C using a combination of Backward Differentiation Formula
C (BDF) methods and a choice of two linear system solution
C methods: direct (dense or band) or Krylov (iterative).
C This version is in double precision.
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C *Usage:
C
C IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR(*)
C DOUBLE PRECISION T, Y(*), YPRIME(*), TOUT, RTOL(*), ATOL(*),
C RWORK(LRW), RPAR(*)
C EXTERNAL RES, JAC, PSOL
C
C CALL DDASPK (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
C * IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC, PSOL)
C
C Quantities which may be altered by the code are:
C T, Y(*), YPRIME(*), INFO(1), RTOL, ATOL, IDID, RWORK(*), IWORK(*)
C
C
C *Arguments:
C
C RES:EXT This is the name of a subroutine which you
C provide to define the residual function G(t,y,y')
C of the differential/algebraic system.
C
C NEQ:IN This is the number of equations in the system.
C
C T:INOUT This is the current value of the independent
C variable.
C
C Y(*):INOUT This array contains the solution components at T.
C
C YPRIME(*):INOUT This array contains the derivatives of the solution
C components at T.
C
C TOUT:IN This is a point at which a solution is desired.
C
C INFO(N):IN This is an integer array used to communicate details
C of how the solution is to be carried out, such as
C tolerance type, matrix structure, step size and
C order limits, and choice of nonlinear system method.
C N must be at least 20.
C
C RTOL,ATOL:INOUT These quantities represent absolute and relative
C error tolerances (on local error) which you provide
C to indicate how accurately you wish the solution to
C be computed. You may choose them to be both scalars
C or else both arrays of length NEQ.
C
C IDID:OUT This integer scalar is an indicator reporting what
C the code did. You must monitor this variable to
C decide what action to take next.
C
C RWORK:WORK A real work array of length LRW which provides the
C code with needed storage space.
C
C LRW:IN The length of RWORK.
C
C IWORK:WORK An integer work array of length LIW which provides
C the code with needed storage space.
C
C LIW:IN The length of IWORK.
C
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, JAC, and PSOL subroutines.
C
C JAC:EXT This is the name of a subroutine which you may
C provide (optionally) for calculating Jacobian
C (partial derivative) data involved in solving linear
C systems within DDASPK.
C
C PSOL:EXT This is the name of a subroutine which you must
C provide for solving linear systems if you selected
C a Krylov method. The purpose of PSOL is to solve
C linear systems involving a left preconditioner P.
C
C *Overview
C
C The DDASPK solver uses the backward differentiation formulas of
C orders one through five to solve a system of the form G(t,y,y') = 0
C for y = Y and y' = YPRIME. Values for Y and YPRIME at the initial
C time must be given as input. These values should be consistent,
C that is, if T, Y, YPRIME are the given initial values, they should
C satisfy G(T,Y,YPRIME) = 0. However, if consistent values are not
C known, in many cases you can have DDASPK solve for them -- see INFO(11).
C (This and other options are described in more detail below.)
C
C Normally, DDASPK solves the system from T to TOUT. It is easy to
C continue the solution to get results at additional TOUT. This is
C the interval mode of operation. Intermediate results can also be
C obtained easily by specifying INFO(3).
C
C On each step taken by DDASPK, a sequence of nonlinear algebraic
C systems arises. These are solved by one of two types of
C methods:
C * a Newton iteration with a direct method for the linear
C systems involved (INFO(12) = 0), or
C * a Newton iteration with a preconditioned Krylov iterative
C method for the linear systems involved (INFO(12) = 1).
C
C The direct method choices are dense and band matrix solvers,
C with either a user-supplied or an internal difference quotient
C Jacobian matrix, as specified by INFO(5) and INFO(6).
C In the band case, INFO(6) = 1, you must supply half-bandwidths
C in IWORK(1) and IWORK(2).
C
C The Krylov method is the Generalized Minimum Residual (GMRES)
C method, in either complete or incomplete form, and with
C scaling and preconditioning. The method is implemented
C in an algorithm called SPIGMR. Certain options in the Krylov
C method case are specified by INFO(13) and INFO(15).
C
C If the Krylov method is chosen, you may supply a pair of routines,
C JAC and PSOL, to apply preconditioning to the linear system.
C If the system is A*x = b, the matrix is A = dG/dY + CJ*dG/dYPRIME
C (of order NEQ). This system can then be preconditioned in the form
C (P-inverse)*A*x = (P-inverse)*b, with left preconditioner P.
C (DDASPK does not allow right preconditioning.)
C Then the Krylov method is applied to this altered, but equivalent,
C linear system, hopefully with much better performance than without
C preconditioning. (In addition, a diagonal scaling matrix based on
C the tolerances is also introduced into the altered system.)
C
C The JAC routine evaluates any data needed for solving systems
C with coefficient matrix P, and PSOL carries out that solution.
C In any case, in order to improve convergence, you should try to
C make P approximate the matrix A as much as possible, while keeping
C the system P*x = b reasonably easy and inexpensive to solve for x,
C given a vector b.
C
C
C *Description
C
C------INPUT - WHAT TO DO ON THE FIRST CALL TO DDASPK-------------------
C
C
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
C
C RES -- Provide a subroutine of the form
C
C SUBROUTINE RES (T, Y, YPRIME, CJ, DELTA, IRES, RPAR, IPAR)
C
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 return
C the residual of the differential/algebraic system
C DELTA = G(T,Y,YPRIME)
C DELTA is a vector of length NEQ which is output from RES.
C
C Subroutine RES must not alter T, Y, YPRIME, or CJ.
C You must declare the name RES in an EXTERNAL
C statement in your program that calls DDASPK.
C You must dimension Y, YPRIME, and DELTA in RES.
C
C The input argument CJ can be ignored, or used to rescale
C constraint equations in the system (see Ref. 2, p. 145).
C Note: In this respect, DDASPK is not downward-compatible
C with DDASSL, which does not have the RES argument CJ.
C
C IRES is an integer flag which is always equal to zero
C on input. Subroutine RES should alter IRES only if it
C encounters an illegal value of Y or a stop condition.
C Set IRES = -1 if an input value is illegal, and DDASPK
C will try to solve the problem without getting IRES = -1.
C If IRES = -2, DDASPK will return control to the calling
C program with IDID = -11.
C
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 DDASPK. 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
C NEQ -- Set it to the number of equations in the system (NEQ .GE. 1).
C
C T -- Set it to the initial point of the integration. (T must be
C a variable.)
C
C Y(*) -- Set this array 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
C YPRIME(*) -- Set this array to the initial values of the NEQ first
C derivatives of the solution components at the initial
C point. You must dimension YPRIME at least NEQ in your
C calling program.
C
C TOUT - Set it to the first point at which a solution is desired.
C You cannot take TOUT = T. Integration either forward in T
C (TOUT .GT. T) or backward in T (TOUT .LT. T) is permitted.
C
C The code advances the solution from T to TOUT using step
C sizes which are automatically selected so as to achieve the
C desired accuracy. If you wish, the code will return with the
C solution and its derivative at intermediate steps (the
C intermediate-output mode) so that you can monitor them,
C but you still must provide TOUT in accord with the basic
C aim of the code.
C
C The first step taken by the code is a critical one because
C it must reflect how fast the solution changes near the
C initial point. The code automatically selects an initial
C step size which is practically always suitable for the
C problem. By using the fact that the code will not step past
C TOUT in the first step, you could, if necessary, restrict the
C length of the initial step.
C
C For some problems it may not be permissible 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 past
C TSTOP. In this case any tout beyond TSTOP is invalid input.
C
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 20, though DDASPK uses only the
C first 15 entries. You must respond to all of the following
C items, which are arranged as questions. The simplest use
C of DDASPK corresponds to setting all entries of INFO to 0.
C
C INFO(1) - This parameter enables the code to initialize itself.
C You must set it to indicate the start of every new
C problem.
C
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
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 arrays.
C The code must be told your choice.
C
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
C INFO(3) - The code integrates from T in the direction of TOUT
C by steps. If you wish, it will return the computed
C solution and derivative at the next intermediate step
C (the intermediate-output mode) or TOUT, whichever comes
C first. This is a good way to proceed if you want to
C see the behavior of the solution. If you must have
C solutions at a great many specific TOUT points, this
C code will compute them efficiently.
C
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
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 this stop condition.
C
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
C INFO(5) - used only when INFO(12) = 0 (direct methods).
C To solve differential/algebraic systems you may wish
C 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. Usually numerical differencing is
C more costly than evaluating derivatives in JAC, but
C sometimes it is not - this depends on your problem.
C
C **** Do you want the code to evaluate the partial deriv-
C atives 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
C INFO(6) - used only when INFO(12) = 0 (direct methods).
C DDASPK will perform much better if the matrix of
C partial derivatives, dG/dY + CJ*dG/dYPRIME (here CJ is
C a scalar determined by DDASPK), is banded and the code
C is told this. In this case, the storage needed will be
C greatly reduced, numerical differencing will be performed
C much cheaper, and a number of important algorithms will
C execute much faster. The differential equation is said
C to have half-bandwidths ML (lower) and MU (upper) if
C equation i 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
C **** Do you want to solve the problem using a full (dense)
C matrix (and not a special banded 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
C INFO(7) - You can specify a maximum (absolute value of)
C stepsize, so that the code will avoid passing over very
C large regions.
C
C **** Do you want the code to decide on its own the maximum
C stepsize ...
C yes - set INFO(7) = 0
C no - set INFO(7) = 1
C and define HMAX by setting
C RWORK(2) = HMAX ****
C
C INFO(8) - Differential/algebraic problems may occasionally
C suffer from severe scaling difficulties on the first
C step. If you know a great deal about the scaling of
C your problem, you can help to alleviate this problem
C by specifying an initial stepsize H0.
C
C **** Do you want the code to define its own initial
C stepsize ...
C yes - set INFO(8) = 0
C no - set INFO(8) = 1
C and define H0 by setting
C RWORK(3) = H0 ****
C
C INFO(9) - If storage is a severe problem, you can save some
C storage by restricting the maximum method order MAXORD.
C The default value is 5. For each order decrease below 5,
C the code requires NEQ fewer locations, but it is likely
C to be slower. In any case, you must have
C 1 .le. MAXORD .le. 5.
C **** Do you want the maximum order to 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
C INFO(10) - If you know that certain components of the
C solutions to your equations are always nonnegative
C (or nonpositive), it may help to set this
C parameter. There are three options that are
C available:
C 1. To have constraint checking only in the initial
C condition calculation.
C 2. To enforce nonnegativity in Y during the integration.
C 3. To enforce both options 1 and 2.
C
C When selecting option 2 or 3, it is probably best to try the
C code without using this option first, and only use
C this option if that does not work very well.
C
C **** Do you want the code to solve the problem without
C invoking any special inequality constraints ...
C yes - set INFO(10) = 0
C no - set INFO(10) = 1 to have option 1 enforced
C no - set INFO(10) = 2 to have option 2 enforced
C no - set INFO(10) = 3 to have option 3 enforced ****
C
C If you have specified INFO(10) = 1 or 3, then you
C will also need to identify how each component of Y
C in the initial condition calculation is constrained.
C You must set:
C IWORK(40+I) = +1 if Y(I) must be .GE. 0,
C IWORK(40+I) = +2 if Y(I) must be .GT. 0,
C IWORK(40+I) = -1 if Y(I) must be .LE. 0, while
C IWORK(40+I) = -2 if Y(I) must be .LT. 0, while
C IWORK(40+I) = 0 if Y(I) is not constrained.
C
C INFO(11) - DDASPK normally requires the initial T, Y, and
C YPRIME to be consistent. That is, you must have
C G(T,Y,YPRIME) = 0 at the initial T. If you do not know
C the initial conditions precisely, in some cases
C DDASPK may be able to compute it.
C
C Denoting the differential variables in Y by Y_d
C and the algebraic variables by Y_a, DDASPK can solve
C one of two initialization problems:
C 1. Given Y_d, calculate Y_a and Y'_d, or
C 2. Given Y', calculate Y.
C In either case, initial values for the given
C components are input, and initial guesses for
C the unknown components must also be provided as input.
C
C **** Are the initial T, Y, YPRIME consistent ...
C
C yes - set INFO(11) = 0
C no - set INFO(11) = 1 to calculate option 1 above,
C or set INFO(11) = 2 to calculate option 2 ****
C
C If you have specified INFO(11) = 1, then you
C will also need to identify which are the
C differential and which are the algebraic
C components (algebraic components are components
C whose derivatives do not appear explicitly
C in the function G(T,Y,YPRIME)). You must set:
C IWORK(LID+I) = +1 if Y(I) is a differential variable
C IWORK(LID+I) = -1 if Y(I) is an algebraic variable,
C where LID = 40 if INFO(10) = 0 or 2 and LID = 40+NEQ
C if INFO(10) = 1 or 3.
C
C INFO(12) - Except for the addition of the RES argument CJ,
C DDASPK by default is downward-compatible with DDASSL,
C which uses only direct (dense or band) methods to solve
C the linear systems involved. You must set INFO(12) to
C indicate whether you want the direct methods or the
C Krylov iterative method.
C **** Do you want DDASPK to use standard direct methods
C (dense or band) or the Krylov (iterative) method ...
C direct methods - set INFO(12) = 0.
C Krylov method - set INFO(12) = 1,
C and check the settings of INFO(13) and INFO(15).
C
C INFO(13) - used when INFO(12) = 1 (Krylov methods).
C DDASPK uses scalars MAXL, KMP, NRMAX, and EPLI for the
C iterative solution of linear systems. INFO(13) allows
C you to override the default values of these parameters.
C These parameters and their defaults are as follows:
C MAXL = maximum number of iterations in the SPIGMR
C algorithm (MAXL .le. NEQ). The default is
C MAXL = MIN(5,NEQ).
C KMP = number of vectors on which orthogonalization is
C done in the SPIGMR algorithm. The default is
C KMP = MAXL, which corresponds to complete GMRES
C iteration, as opposed to the incomplete form.
C NRMAX = maximum number of restarts of the SPIGMR
C algorithm per nonlinear iteration. The default is
C NRMAX = 5.
C EPLI = convergence test constant in SPIGMR algorithm.
C The default is EPLI = 0.05.
C Note that the length of RWORK depends on both MAXL
C and KMP. See the definition of LRW below.
C **** Are MAXL, KMP, and EPLI to be given their
C default values ...
C yes - set INFO(13) = 0
C no - set INFO(13) = 1,
C and set all of the following:
C IWORK(24) = MAXL (1 .le. MAXL .le. NEQ)
C IWORK(25) = KMP (1 .le. KMP .le. MAXL)
C IWORK(26) = NRMAX (NRMAX .ge. 0)
C RWORK(10) = EPLI (0 .lt. EPLI .lt. 1.0) ****
C
C INFO(14) - used with INFO(11) > 0 (initial condition
C calculation is requested). In this case, you may
C request control to be returned to the calling program
C immediately after the initial condition calculation,
C before proceeding to the integration of the system
C (e.g. to examine the computed Y and YPRIME).
C If this is done, and if the initialization succeeded
C (IDID = 4), you should reset INFO(11) to 0 for the
C next call, to prevent the solver from repeating the
C initialization (and to avoid an infinite loop).
C **** Do you want to proceed to the integration after
C the initial condition calculation is done ...
C yes - set INFO(14) = 0
C no - set INFO(14) = 1 ****
C
C INFO(15) - used when INFO(12) = 1 (Krylov methods).
C When using preconditioning in the Krylov method,
C you must supply a subroutine, PSOL, which solves the
C associated linear systems using P.
C The usage of DDASPK is simpler if PSOL can carry out
C the solution without any prior calculation of data.
C However, if some partial derivative data is to be
C calculated in advance and used repeatedly in PSOL,
C then you must supply a JAC routine to do this,
C and set INFO(15) to indicate that JAC is to be called
C for this purpose. For example, P might be an
C approximation to a part of the matrix A which can be
C calculated and LU-factored for repeated solutions of
C the preconditioner system. The arrays WP and IWP
C (described under JAC and PSOL) can be used to
C communicate data between JAC and PSOL.
C **** Does PSOL operate with no prior preparation ...
C yes - set INFO(15) = 0 (no JAC routine)
C no - set INFO(15) = 1
C and supply a JAC routine to evaluate and
C preprocess any required Jacobian data. ****
C
C INFO(16) - option to exclude algebraic variables from
C the error test.
C **** Do you wish to control errors locally on
C all the variables...
C yes - set INFO(16) = 0
C no - set INFO(16) = 1
C If you have specified INFO(16) = 1, then you
C will also need to identify which are the
C differential and which are the algebraic
C components (algebraic components are components
C whose derivatives do not appear explicitly
C in the function G(T,Y,YPRIME)). You must set:
C IWORK(LID+I) = +1 if Y(I) is a differential
C variable, and
C IWORK(LID+I) = -1 if Y(I) is an algebraic
C variable,
C where LID = 40 if INFO(10) = 0 or 2 and
C LID = 40 + NEQ if INFO(10) = 1 or 3.
C
C INFO(17) - used when INFO(11) > 0 (DDASPK is to do an
C initial condition calculation).
C DDASPK uses several heuristic control quantities in the
C initial condition calculation. They have default values,
C but can also be set by the user using INFO(17).
C These parameters and their defaults are as follows:
C MXNIT = maximum number of Newton iterations
C per Jacobian or preconditioner evaluation.
C The default is:
C MXNIT = 5 in the direct case (INFO(12) = 0), and
C MXNIT = 15 in the Krylov case (INFO(12) = 1).
C MXNJ = maximum number of Jacobian or preconditioner
C evaluations. The default is:
C MXNJ = 6 in the direct case (INFO(12) = 0), and
C MXNJ = 2 in the Krylov case (INFO(12) = 1).
C MXNH = maximum number of values of the artificial
C stepsize parameter H to be tried if INFO(11) = 1.
C The default is MXNH = 5.
C NOTE: the maximum number of Newton iterations
C allowed in all is MXNIT*MXNJ*MXNH if INFO(11) = 1,
C and MXNIT*MXNJ if INFO(11) = 2.
C LSOFF = flag to turn off the linesearch algorithm
C (LSOFF = 0 means linesearch is on, LSOFF = 1 means
C it is turned off). The default is LSOFF = 0.
C STPTOL = minimum scaled step in linesearch algorithm.
C The default is STPTOL = (unit roundoff)**(2/3).
C EPINIT = swing factor in the Newton iteration convergence
C test. The test is applied to the residual vector,
C premultiplied by the approximate Jacobian (in the
C direct case) or the preconditioner (in the Krylov
C case). For convergence, the weighted RMS norm of
C this vector (scaled by the error weights) must be
C less than EPINIT*EPCON, where EPCON = .33 is the
C analogous test constant used in the time steps.
C The default is EPINIT = .01.
C **** Are the initial condition heuristic controls to be
C given their default values...
C yes - set INFO(17) = 0
C no - set INFO(17) = 1,
C and set all of the following:
C IWORK(32) = MXNIT (.GT. 0)
C IWORK(33) = MXNJ (.GT. 0)
C IWORK(34) = MXNH (.GT. 0)
C IWORK(35) = LSOFF ( = 0 or 1)
C RWORK(14) = STPTOL (.GT. 0.0)
C RWORK(15) = EPINIT (.GT. 0.0) ****
C
C INFO(18) - option to get extra printing in initial condition
C calculation.
C **** Do you wish to have extra printing...
C no - set INFO(18) = 0
C yes - set INFO(18) = 1 for minimal printing, or
C set INFO(18) = 2 for full printing.
C If you have specified INFO(18) .ge. 1, data
C will be printed with the error handler routines.
C To print to a non-default unit number L, include
C the line CALL XSETUN(L) in your program. ****
C
C RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL)
C error tolerances to tell the code how accurately you
C want the solution to be computed. They must be defined
C as variables because the code may change them.
C you have two choices --
C Both RTOL and ATOL are scalars (INFO(2) = 0), or
C both RTOL and ATOL are vectors (INFO(2) = 1).
C In either case all components must be non-negative.
C
C The tolerances are used by the code in a local error
C test at each step which requires roughly that
C abs(local error in Y(i)) .le. EWT(i) ,
C where EWT(i) = RTOL*abs(Y(i)) + ATOL is an error weight
C quantity, for each vector component.
C (More specifically, a root-mean-square norm is used to
C measure the size of vectors, and the error test uses the
C magnitude of the solution at the beginning of the step.)
C
C The true (global) error is the difference between the
C true solution of the initial value problem and the
C computed approximation. Practically all present day
C codes, including this one, control the local error at
C each step and do not even attempt to control the global
C error directly.
C
C Usually, but not always, the true accuracy of
C the computed Y is comparable to the error tolerances.
C This code will usually, but not always, deliver a more
C accurate solution if you reduce the tolerances and
C integrate again. By comparing two such solutions you
C can get a fairly reliable idea of the true error in the
C solution at the larger tolerances.
C
C Setting ATOL = 0. results in a pure relative error test
C on that component. Setting RTOL = 0. results in a pure
C absolute error test on that component. A mixed test
C with non-zero RTOL and ATOL corresponds roughly to a
C relative error test when the solution component is
C much bigger than ATOL and to an absolute error test
C when the solution component is smaller than the
C threshold ATOL.
C
C The code will not attempt to compute a solution at an
C accuracy unreasonable for the machine being used. It
C will advise you if you ask for too much accuracy and
C inform you as to the maximum accuracy it believes
C possible.
C
C RWORK(*) -- a real work array, which should be dimensioned in your
C calling program with a length equal to the value of
C LRW (or greater).
C
C LRW -- Set it to the declared length of the RWORK array. The
C minimum length depends on the options you have selected,
C given by a base value plus additional storage as described
C below.
C
C If INFO(12) = 0 (standard direct method), the base value is
C base = 50 + max(MAXORD+4,7)*NEQ.
C The default value is MAXORD = 5 (see INFO(9)). With the
C default MAXORD, base = 50 + 9*NEQ.
C Additional storage must be added to the base value for
C any or all of the following options:
C if INFO(6) = 0 (dense matrix), add NEQ**2
C if INFO(6) = 1 (banded matrix), then
C if INFO(5) = 0, add (2*ML+MU+1)*NEQ + 2*(NEQ/(ML+MU+1)+1),
C if INFO(5) = 1, add (2*ML+MU+1)*NEQ,
C if INFO(16) = 1, add NEQ.
C
C If INFO(12) = 1 (Krylov method), the base value is
C base = 50 + (MAXORD+5)*NEQ + (MAXL+3+MIN0(1,MAXL-KMP))*NEQ +
C + (MAXL+3)*MAXL + 1 + LENWP.
C See PSOL for description of LENWP. The default values are:
C MAXORD = 5 (see INFO(9)), MAXL = min(5,NEQ) and KMP = MAXL
C (see INFO(13)).
C With the default values for MAXORD, MAXL and KMP,
C base = 91 + 18*NEQ + LENWP.
C Additional storage must be added to the base value for
C any or all of the following options:
C if INFO(16) = 1, add NEQ.
C
C
C IWORK(*) -- an integer work array, which should be dimensioned in
C your calling program with a length equal to the value
C of LIW (or greater).
C
C LIW -- Set it to the declared length of the IWORK array. The
C minimum length depends on the options you have selected,
C given by a base value plus additional storage as described
C below.
C
C If INFO(12) = 0 (standard direct method), the base value is
C base = 40 + NEQ.
C IF INFO(10) = 1 or 3, add NEQ to the base value.
C If INFO(11) = 1 or INFO(16) =1, add NEQ to the base value.
C
C If INFO(12) = 1 (Krylov method), the base value is
C base = 40 + LENIWP.
C See PSOL for description of LENIWP.
C IF INFO(10) = 1 or 3, add NEQ to the base value.
C If INFO(11) = 1 or INFO(16) = 1, add NEQ to the base value.
C
C
C RPAR, IPAR -- These are arrays of double precision and integer type,
C respectively, which are available for you to use
C for communication between your program that calls
C DDASPK and the RES subroutine (and the JAC and PSOL
C subroutines). They are not altered by DDASPK.
C If you do not need RPAR or IPAR, ignore these
C parameters by treating them as dummy arguments.
C If you do choose to use them, dimension them in
C your calling program and in RES (and in JAC and PSOL)
C as arrays of appropriate length.
C
C JAC -- This is the name of a routine that you may supply
C (optionally) that relates to the Jacobian matrix of the
C nonlinear system that the code must solve at each T step.
C The role of JAC (and its call sequence) depends on whether
C a direct (INFO(12) = 0) or Krylov (INFO(12) = 1) method
C is selected.
C
C **** INFO(12) = 0 (direct methods):
C If you are letting the code generate partial derivatives
C numerically (INFO(5) = 0), then JAC can be absent
C (or perhaps a dummy routine to satisfy the loader).
C Otherwise you must supply a JAC routine to compute
C the matrix A = dG/dY + CJ*dG/dYPRIME. It must have
C the form
C
C SUBROUTINE JAC (T, Y, YPRIME, PD, CJ, RPAR, IPAR)
C
C The JAC routine must dimension Y, YPRIME, and PD (and RPAR
C and IPAR if used). CJ is a scalar which is input to JAC.
C For the given values of T, Y, and YPRIME, the JAC routine
C must evaluate the nonzero elements of the matrix A, and
C store these values in the array PD. The elements of PD are
C set to zero before each call to JAC, so that only nonzero
C elements need to be defined.
C The way you store the elements into the PD array depends
C on the structure of the matrix indicated by INFO(6).
C *** INFO(6) = 0 (full or dense matrix) ***
C Give PD a first dimension of NEQ. When you evaluate the
C nonzero partial derivatives of equation i (i.e. of G(i))
C with respect to component j (of Y and YPRIME), you must
C store the element in PD according to
C PD(i,j) = dG(i)/dY(j) + CJ*dG(i)/dYPRIME(j).
C *** INFO(6) = 1 (banded matrix with half-bandwidths ML, MU
C as described under INFO(6)) ***
C Give PD a first dimension of 2*ML+MU+1. When you
C evaluate the nonzero partial derivatives of equation i
C (i.e. of G(i)) with respect to component j (of Y and
C YPRIME), you must store the element in PD according to
C IROW = i - j + ML + MU + 1
C PD(IROW,j) = dG(i)/dY(j) + CJ*dG(i)/dYPRIME(j).
C
C **** INFO(12) = 1 (Krylov method):
C If you are not calculating Jacobian data in advance for use
C in PSOL (INFO(15) = 0), JAC can be absent (or perhaps a
C dummy routine to satisfy the loader). Otherwise, you may
C supply a JAC routine to compute and preprocess any parts of
C of the Jacobian matrix A = dG/dY + CJ*dG/dYPRIME that are
C involved in the preconditioner matrix P.
C It is to have the form
C
C SUBROUTINE JAC (RES, IRES, NEQ, T, Y, YPRIME, REWT, SAVR,
C WK, H, CJ, WP, IWP, IER, RPAR, IPAR)
C
C The JAC routine must dimension Y, YPRIME, REWT, SAVR, WK,
C and (if used) WP, IWP, RPAR, and IPAR.
C The Y, YPRIME, and SAVR arrays contain the current values
C of Y, YPRIME, and the residual G, respectively.
C The array WK is work space of length NEQ.
C H is the step size. CJ is a scalar, input to JAC, that is
C normally proportional to 1/H. REWT is an array of
C reciprocal error weights, 1/EWT(i), where EWT(i) is
C RTOL*abs(Y(i)) + ATOL (unless you supplied routine DDAWTS
C instead), for use in JAC if needed. For example, if JAC
C computes difference quotient approximations to partial
C derivatives, the REWT array may be useful in setting the
C increments used. The JAC routine should do any
C factorization operations called for, in preparation for
C solving linear systems in PSOL. The matrix P should
C be an approximation to the Jacobian,
C A = dG/dY + CJ*dG/dYPRIME.
C
C WP and IWP are real and integer work arrays which you may
C use for communication between your JAC routine and your
C PSOL routine. These may be used to store elements of the
C preconditioner P, or related matrix data (such as factored
C forms). They are not altered by DDASPK.
C If you do not need WP or IWP, ignore these parameters by
C treating them as dummy arguments. If you do use them,
C dimension them appropriately in your JAC and PSOL routines.
C See the PSOL description for instructions on setting
C the lengths of WP and IWP.
C
C On return, JAC should set the error flag IER as follows..
C IER = 0 if JAC was successful,
C IER .ne. 0 if JAC was unsuccessful (e.g. if Y or YPRIME
C was illegal, or a singular matrix is found).
C (If IER .ne. 0, a smaller stepsize will be tried.)
C IER = 0 on entry to JAC, so need be reset only on a failure.
C If RES is used within JAC, then a nonzero value of IRES will
C override any nonzero value of IER (see the RES description).
C
C Regardless of the method type, subroutine JAC must not
C alter T, Y(*), YPRIME(*), H, CJ, or REWT(*).
C You must declare the name JAC in an EXTERNAL statement in
C your program that calls DDASPK.
C
C PSOL -- This is the name of a routine you must supply if you have
C selected a Krylov method (INFO(12) = 1) with preconditioning.
C In the direct case (INFO(12) = 0), PSOL can be absent
C (a dummy routine may have to be supplied to satisfy the
C loader). Otherwise, you must provide a PSOL routine to
C solve linear systems arising from preconditioning.
C When supplied with INFO(12) = 1, the PSOL routine is to
C have the form
C
C SUBROUTINE PSOL (NEQ, T, Y, YPRIME, SAVR, WK, CJ, WGHT,
C WP, IWP, B, EPLIN, IER, RPAR, IPAR)
C
C The PSOL routine must solve linear systems of the form
C P*x = b where P is the left preconditioner matrix.
C
C The right-hand side vector b is in the B array on input, and
C PSOL must return the solution vector x in B.
C The Y, YPRIME, and SAVR arrays contain the current values
C of Y, YPRIME, and the residual G, respectively.
C
C Work space required by JAC and/or PSOL, and space for data to
C be communicated from JAC to PSOL is made available in the form
C of arrays WP and IWP, which are parts of the RWORK and IWORK
C arrays, respectively. The lengths of these real and integer
C work spaces WP and IWP must be supplied in LENWP and LENIWP,
C respectively, as follows..
C IWORK(27) = LENWP = length of real work space WP
C IWORK(28) = LENIWP = length of integer work space IWP.
C
C WK is a work array of length NEQ for use by PSOL.
C CJ is a scalar, input to PSOL, that is normally proportional
C to 1/H (H = stepsize). If the old value of CJ
C (at the time of the last JAC call) is needed, it must have
C been saved by JAC in WP.
C
C WGHT is an array of weights, to be used if PSOL uses an
C iterative method and performs a convergence test. (In terms
C of the argument REWT to JAC, WGHT is REWT/sqrt(NEQ).)
C If PSOL uses an iterative method, it should use EPLIN
C (a heuristic parameter) as the bound on the weighted norm of
C the residual for the computed solution. Specifically, the
C residual vector R should satisfy
C SQRT (SUM ( (R(i)*WGHT(i))**2 ) ) .le. EPLIN
C
C PSOL must not alter NEQ, T, Y, YPRIME, SAVR, CJ, WGHT, EPLIN.
C
C On return, PSOL should set the error flag IER as follows..
C IER = 0 if PSOL was successful,
C IER .lt. 0 if an unrecoverable error occurred, meaning
C control will be passed to the calling routine,
C IER .gt. 0 if a recoverable error occurred, meaning that
C the step will be retried with the same step size
C but with a call to JAC to update necessary data,
C unless the Jacobian data is current, in which case
C the step will be retried with a smaller step size.
C IER = 0 on entry to PSOL so need be reset only on a failure.
C
C You must declare the name PSOL in an EXTERNAL statement in
C your program that calls DDASPK.
C
C
C OPTIONALLY REPLACEABLE SUBROUTINE:
C
C DDASPK uses a weighted root-mean-square norm to measure the
C size of various error vectors. The weights used in this norm
C are set in the following subroutine:
C
C SUBROUTINE DDAWTS (NEQ, IWT, RTOL, ATOL, Y, EWT, RPAR, IPAR)
C DIMENSION RTOL(*), ATOL(*), Y(*), EWT(*), RPAR(*), IPAR(*)
C
C A DDAWTS routine has been included with DDASPK which sets the
C weights according to
C EWT(I) = RTOL*ABS(Y(I)) + ATOL
C in the case of scalar tolerances (IWT = 0) or
C EWT(I) = RTOL(I)*ABS(Y(I)) + ATOL(I)
C in the case of array tolerances (IWT = 1). (IWT is INFO(2).)
C In some special cases, it may be appropriate for you to define
C your own error weights by writing a subroutine DDAWTS to be
C called instead of the version supplied. However, this should
C be attempted only after careful thought and consideration.
C If you supply this routine, you may use the tolerances and Y
C as appropriate, but do not overwrite these variables. You
C may also use RPAR and IPAR to communicate data as appropriate.
C ***Note: Aside from the values of the weights, the choice of
C norm used in DDASPK (weighted root-mean-square) is not subject
C to replacement by the user. In this respect, DDASPK is not
C downward-compatible with the original DDASSL solver (in which
C the norm routine was optionally user-replaceable).
C
C
C------OUTPUT - AFTER ANY RETURN FROM DDASPK----------------------------
C
C The principal aim of the code is to return a computed solution at
C T = TOUT, although it is also possible to obtain intermediate
C results along the way. To find out whether the code achieved its
C goal or if the integration process was interrupted before the task
C was completed, you must check the IDID parameter.
C
C
C T -- The output value of T is the point to which the solution
C was successfully advanced.
C
C Y(*) -- contains the computed solution approximation at T.
C
C YPRIME(*) -- contains the computed derivative approximation at T.
C
C IDID -- reports what the code did, described as follows:
C
C *** TASK COMPLETED ***
C Reported by positive values of IDID
C
C IDID = 1 -- a step was successfully taken in the
C intermediate-output mode. The code has not
C yet reached TOUT.
C
C IDID = 2 -- the integration to TSTOP was successfully
C completed (T = TSTOP) by stepping exactly to TSTOP.
C
C IDID = 3 -- the integration to TOUT was successfully
C completed (T = TOUT) by stepping past TOUT.
C Y(*) and YPRIME(*) are obtained by interpolation.
C
C IDID = 4 -- the initial condition calculation, with
C INFO(11) > 0, was successful, and INFO(14) = 1.
C No integration steps were taken, and the solution
C is not considered to have been started.
C
C *** TASK INTERRUPTED ***
C Reported by negative values of IDID
C
C IDID = -1 -- a large amount of work has been expended
C (about 500 steps).
C
C IDID = -2 -- the error tolerances are too stringent.
C
C IDID = -3 -- the local error test cannot be satisfied
C because you specified a zero component in ATOL
C and the corresponding computed solution component
C is zero. Thus, a pure relative error test is
C impossible for this component.
C
C IDID = -5 -- there were repeated failures in the evaluation
C or processing of the preconditioner (in JAC).
C
C IDID = -6 -- DDASPK had repeated error test failures on the
C last attempted step.
C
C IDID = -7 -- the nonlinear system solver in the time integration
C could not converge.
C
C IDID = -8 -- the matrix of partial derivatives appears
C to be singular (direct method).
C
C IDID = -9 -- the nonlinear system solver in the time integration
C failed to achieve convergence, and there were repeated
C error test failures in this step.
C
C IDID =-10 -- the nonlinear system solver in the time integration
C failed to achieve convergence because IRES was equal
C to -1.
C
C IDID =-11 -- IRES = -2 was encountered and control is
C being returned to the calling program.
C
C IDID =-12 -- DDASPK failed to compute the initial Y, YPRIME.
C
C IDID =-13 -- unrecoverable error encountered inside user's
C PSOL routine, and control is being returned to
C the calling program.
C
C IDID =-14 -- the Krylov linear system solver could not
C achieve convergence.
C
C IDID =-15,..,-32 -- Not applicable for this code.
C
C *** TASK TERMINATED ***
C reported by the value of IDID=-33
C
C IDID = -33 -- the code has encountered trouble from which
C it cannot recover. A message is printed
C explaining the trouble and control is returned
C to the calling program. For example, this occurs
C when invalid input is detected.
C
C RTOL, ATOL -- these quantities remain unchanged except when
C IDID = -2. In this case, the error tolerances have been
C increased by the code to values which are estimated to
C be appropriate for continuing the integration. However,
C the reported solution at T was obtained using the input
C values of RTOL and ATOL.
C
C RWORK, IWORK -- contain information which is usually of no interest
C to the user but necessary for subsequent calls.
C However, you may be interested in the performance data
C listed below. These quantities are accessed in RWORK
C and IWORK but have internal mnemonic names, as follows..
C
C RWORK(3)--contains H, the step size h to be attempted
C on the next step.
C
C RWORK(4)--contains TN, the current value of the
C independent variable, i.e. the farthest point
C integration has reached. This will differ
C from T if interpolation has been performed
C (IDID = 3).
C
C RWORK(7)--contains HOLD, the stepsize used on the last
C successful step. If INFO(11) = INFO(14) = 1,
C this contains the value of H used in the
C initial condition calculation.
C
C IWORK(7)--contains K, the order of the method to be
C attempted on the next step.
C
C IWORK(8)--contains KOLD, the order of the method used
C on the last step.
C
C IWORK(11)--contains NST, the number of steps (in T)
C taken so far.
C
C IWORK(12)--contains NRE, the number of calls to RES
C so far.
C
C IWORK(13)--contains NJE, the number of calls to JAC so
C far (Jacobian or preconditioner evaluations).
C
C IWORK(14)--contains NETF, the total number of error test
C failures so far.
C
C IWORK(15)--contains NCFN, the total number of nonlinear
C convergence failures so far (includes counts
C of singular iteration matrix or singular
C preconditioners).
C
C IWORK(16)--contains NCFL, the number of convergence
C failures of the linear iteration so far.
C
C IWORK(17)--contains LENIW, the length of IWORK actually
C required. This is defined on normal returns
C and on an illegal input return for
C insufficient storage.
C
C IWORK(18)--contains LENRW, the length of RWORK actually
C required. This is defined on normal returns
C and on an illegal input return for
C insufficient storage.
C
C IWORK(19)--contains NNI, the total number of nonlinear
C iterations so far (each of which calls a
C linear solver).
C
C IWORK(20)--contains NLI, the total number of linear
C (Krylov) iterations so far.
C
C IWORK(21)--contains NPS, the number of PSOL calls so
C far, for preconditioning solve operations or
C for solutions with the user-supplied method.
C
C Note: The various counters in IWORK do not include
C counts during a call made with INFO(11) > 0 and
C INFO(14) = 1.
C
C
C------INPUT - WHAT TO DO TO CONTINUE THE INTEGRATION -----------------
C (CALLS AFTER THE FIRST)
C
C This code is organized so that subsequent calls to continue the
C integration involve little (if any) additional effort on your
C part. You must monitor the IDID parameter in order to determine
C what to do next.
C
C Recalling that the principal task of the code is to integrate
C from T to TOUT (the interval mode), usually all you will need
C to do is specify a new TOUT upon reaching the current TOUT.
C
C Do not alter any quantity not specifically permitted below. In
C particular do not alter NEQ, T, Y(*), YPRIME(*), RWORK(*),
C IWORK(*), or the differential equation in subroutine RES. Any
C such alteration constitutes a new problem and must be treated
C as such, i.e. you must start afresh.
C
C You cannot change from array to scalar error control or vice
C versa (INFO(2)), but you can change the size of the entries of
C RTOL or ATOL. Increasing a tolerance makes the equation easier
C to integrate. Decreasing a tolerance will make the equation
C harder to integrate and should generally be avoided.
C
C You can switch from the intermediate-output mode to the
C interval mode (INFO(3)) or vice versa at any time.
C
C If it has been necessary to prevent the integration from going
C past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
C code will not integrate to any TOUT beyond the currently
C specified TSTOP. Once TSTOP has been reached, you must change
C 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 in
C RWORK(1) whenever you set INFO(4) = 1.
C
C Do not change INFO(5), INFO(6), INFO(12-17) or their associated
C IWORK/RWORK locations unless you are going to restart the code.
C
C *** FOLLOWING A COMPLETED TASK ***
C
C If..
C IDID = 1, call the code again to continue the integration
C another step in the direction of TOUT.
C
C IDID = 2 or 3, define a new TOUT and call the code again.
C TOUT must be different from T. You cannot change
C the direction of integration without restarting.
C
C IDID = 4, reset INFO(11) = 0 and call the code again to begin
C the integration. (If you leave INFO(11) > 0 and
C INFO(14) = 1, you may generate an infinite loop.)
C In this situation, the next call to DASPK is
C considered to be the first call for the problem,
C in that all initializations are done.
C
C *** FOLLOWING AN INTERRUPTED TASK ***
C
C To show the code that you realize the task was interrupted and
C that you want to continue, you must take appropriate action and
C set INFO(1) = 1.
C
C If..
C IDID = -1, the code has taken about 500 steps. If you want to
C continue, set INFO(1) = 1 and call the code again.
C An additional 500 steps will be allowed.
C
C
C IDID = -2, the error tolerances RTOL, ATOL have been increased
C to values the code estimates appropriate for
C continuing. You may want to change them yourself.
C If you are sure you want to continue with relaxed
C error tolerances, set INFO(1) = 1 and call the code
C again.
C
C IDID = -3, a solution component is zero and you set the
C corresponding component of ATOL to zero. If you
C are sure you want to continue, you must first alter
C the error criterion to use positive values of ATOL
C for those components corresponding to zero solution
C components, then set INFO(1) = 1 and call the code
C again.
C
C IDID = -4 --- cannot occur with this code.
C
C IDID = -5, your JAC routine failed with the Krylov method. Check
C for errors in JAC and restart the integration.
C
C IDID = -6, repeated error test failures occurred on the last
C attempted step in DDASPK. A singularity in the
C solution may be present. If you are absolutely
C certain you want to continue, you should restart
C the integration. (Provide initial values of Y and
C YPRIME which are consistent.)
C
C IDID = -7, repeated convergence test failures occurred on the last
C attempted step in DDASPK. An inaccurate or ill-
C conditioned Jacobian or preconditioner may be the
C problem. If you are absolutely certain you want
C to continue, you should restart the integration.
C
C
C IDID = -8, the matrix of partial derivatives is singular, with
C the use of direct methods. Some of your equations
C may be redundant. DDASPK cannot solve the problem
C as stated. It is possible that the redundant
C equations could be removed, and then DDASPK could
C solve the problem. It is also possible that a
C solution to your problem either does not exist
C or is not unique.
C
C IDID = -9, DDASPK had multiple convergence test failures, preceded
C by multiple error test failures, on the last
C attempted step. It is possible that your problem is
C ill-posed and cannot be solved using this code. Or,
C there may be a discontinuity or a singularity in the
C solution. If you are absolutely certain you want to
C continue, you should restart the integration.
C
C IDID = -10, DDASPK had multiple convergence test failures
C because IRES was equal to -1. If you are
C absolutely certain you want to continue, you
C should restart the integration.
C
C IDID = -11, there was an unrecoverable error (IRES = -2) from RES
C inside the nonlinear system solver. Determine the
C cause before trying again.
C
C IDID = -12, DDASPK failed to compute the initial Y and YPRIME
C vectors. This could happen because the initial
C approximation to Y or YPRIME was not very good, or
C because no consistent values of these vectors exist.
C The problem could also be caused by an inaccurate or
C singular iteration matrix, or a poor preconditioner.
C
C IDID = -13, there was an unrecoverable error encountered inside
C your PSOL routine. Determine the cause before
C trying again.
C
C IDID = -14, the Krylov linear system solver failed to achieve
C convergence. This may be due to ill-conditioning
C in the iteration matrix, or a singularity in the
C preconditioner (if one is being used).
C Another possibility is that there is a better
C choice of Krylov parameters (see INFO(13)).
C Possibly the failure is caused by redundant equations
C in the system, or by inconsistent equations.
C In that case, reformulate the system to make it
C consistent and non-redundant.
C
C IDID = -15,..,-32 --- Cannot occur with this code.
C
C *** FOLLOWING A TERMINATED TASK ***
C
C If IDID = -33, you cannot continue the solution of this problem.
C An attempt to do so will result in your run being
C terminated.
C
C ---------------------------------------------------------------------
C
C***REFERENCES
C 1. L. R. Petzold, A Description of DASSL: A Differential/Algebraic
C System Solver, in Scientific Computing, R. S. Stepleman et al.
C (Eds.), North-Holland, Amsterdam, 1983, pp. 65-68.
C 2. K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical
C Solution of Initial-Value Problems in Differential-Algebraic
C Equations, Elsevier, New York, 1989.
C 3. P. N. Brown and A. C. Hindmarsh, Reduced Storage Matrix Methods
C in Stiff ODE Systems, J. Applied Mathematics and Computation,
C 31 (1989), pp. 40-91.
C 4. P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Using Krylov
C Methods in the Solution of Large-Scale Differential-Algebraic
C Systems, SIAM J. Sci. Comp., 15 (1994), pp. 1467-1488.
C 5. P. N. Brown, A. C. Hindmarsh, and L. R. Petzold, Consistent
C Initial Condition Calculation for Differential-Algebraic
C Systems, SIAM J. Sci. Comp. 19 (1998), pp. 1495-1512.
C
C***ROUTINES CALLED
C
C The following are all the subordinate routines used by DDASPK.
C
C DDASIC computes consistent initial conditions.
C DYYPNW updates Y and YPRIME in linesearch for initial condition
C calculation.
C DDSTP carries out one step of the integration.
C DCNSTR/DCNST0 check the current solution for constraint violations.
C DDAWTS sets error weight quantities.
C DINVWT tests and inverts the error weights.
C DDATRP performs interpolation to get an output solution.
C DDWNRM computes the weighted root-mean-square norm of a vector.
C D1MACH provides the unit roundoff of the computer.
C XERRWD/XSETF/XSETUN/IXSAV is a package to handle error messages.
C DDASID nonlinear equation driver to initialize Y and YPRIME using
C direct linear system solver methods. Interfaces to Newton
C solver (direct case).
C DNSID solves the nonlinear system for unknown initial values by
C modified Newton iteration and direct linear system methods.
C DLINSD carries out linesearch algorithm for initial condition
C calculation (direct case).
C DFNRMD calculates weighted norm of preconditioned residual in
C initial condition calculation (direct case).
C DNEDD nonlinear equation driver for direct linear system solver
C methods. Interfaces to Newton solver (direct case).
C DMATD assembles the iteration matrix (direct case).
C DNSD solves the associated nonlinear system by modified
C Newton iteration and direct linear system methods.
C DSLVD interfaces to linear system solver (direct case).
C DDASIK nonlinear equation driver to initialize Y and YPRIME using
C Krylov iterative linear system methods. Interfaces to
C Newton solver (Krylov case).
C DNSIK solves the nonlinear system for unknown initial values by
C Newton iteration and Krylov iterative linear system methods.
C DLINSK carries out linesearch algorithm for initial condition
C calculation (Krylov case).
C DFNRMK calculates weighted norm of preconditioned residual in
C initial condition calculation (Krylov case).
C DNEDK nonlinear equation driver for iterative linear system solver
C methods. Interfaces to Newton solver (Krylov case).
C DNSK solves the associated nonlinear system by Inexact Newton
C iteration and (linear) Krylov iteration.
C DSLVK interfaces to linear system solver (Krylov case).
C DSPIGM solves a linear system by SPIGMR algorithm.
C DATV computes matrix-vector product in Krylov algorithm.
C DORTH performs orthogonalization of Krylov basis vectors.
C DHEQR performs QR factorization of Hessenberg matrix.
C DHELS finds least-squares solution of Hessenberg linear system.
C DGEFA, DGESL, DGBFA, DGBSL are LINPACK routines for solving
C linear systems (dense or band direct methods).
C DAXPY, DCOPY, DDOT, DNRM2, DSCAL are Basic Linear Algebra (BLAS)
C routines.
C
C The routines called directly by DDASPK are:
C DCNST0, DDAWTS, DINVWT, D1MACH, DDWNRM, DDASIC, DDATRP, DDSTP,
C XERRWD
C
C***END PROLOGUE DDASPK
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
LOGICAL DONE, LAVL, LCFN, LCFL, LWARN
DIMENSION Y(*),YPRIME(*)
DIMENSION INFO(20)
DIMENSION RWORK(LRW),IWORK(LIW)
DIMENSION RTOL(*),ATOL(*)
DIMENSION RPAR(*),IPAR(*)
CHARACTER MSG*80
EXTERNAL RES, JAC, PSOL, DDASID, DDASIK, DNEDD, DNEDK
C
C Set pointers into IWORK.
C
PARAMETER (LML=1, LMU=2, LMTYPE=4,
* LIWM=1, LMXORD=3, LJCALC=5, LPHASE=6, LK=7, LKOLD=8,
* LNS=9, LNSTL=10, LNST=11, LNRE=12, LNJE=13, LETF=14, LNCFN=15,
* LNCFL=16, LNIW=17, LNRW=18, LNNI=19, LNLI=20, LNPS=21,
* LNPD=22, LMITER=23, LMAXL=24, LKMP=25, LNRMAX=26, LLNWP=27,
* LLNIWP=28, LLOCWP=29, LLCIWP=30, LKPRIN=31,
* LMXNIT=32, LMXNJ=33, LMXNH=34, LLSOFF=35, LICNS=41)
C
C Set pointers into RWORK.
C
PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4, LCJ=5, LCJOLD=6,
* LHOLD=7, LS=8, LROUND=9, LEPLI=10, LSQRN=11, LRSQRN=12,
* LEPCON=13, LSTOL=14, LEPIN=15,
* LALPHA=21, LBETA=27, LGAMMA=33, LPSI=39, LSIGMA=45, LDELTA=51)
C
SAVE LID, LENID, NONNEG, NCPHI
C
C
C***FIRST EXECUTABLE STATEMENT DDASPK
C
C
IF(INFO(1).NE.0) GO TO 100
C
C-----------------------------------------------------------------------
C This block is executed for the initial call only.
C It contains checking of inputs and initializations.
C-----------------------------------------------------------------------
C
C First check INFO array to make sure all elements of INFO
C Are within the proper range. (INFO(1) is checked later, because
C it must be tested on every call.) ITEMP holds the location
C within INFO which may be out of range.
C
DO 10 I=2,9
ITEMP = I
IF (INFO(I) .NE. 0 .AND. INFO(I) .NE. 1) GO TO 701
10 CONTINUE
ITEMP = 10
IF(INFO(10).LT.0 .OR. INFO(10).GT.3) GO TO 701
ITEMP = 11
IF(INFO(11).LT.0 .OR. INFO(11).GT.2) GO TO 701
DO 15 I=12,17
ITEMP = I
IF (INFO(I) .NE. 0 .AND. INFO(I) .NE. 1) GO TO 701
15 CONTINUE
ITEMP = 18
IF(INFO(18).LT.0 .OR. INFO(18).GT.2) GO TO 701
C
C Check NEQ to see if it is positive.
C
IF (NEQ .LE. 0) GO TO 702
C
C Check and compute maximum order.
C
MXORD=5
IF (INFO(9) .NE. 0) THEN
MXORD=IWORK(LMXORD)
IF (MXORD .LT. 1 .OR. MXORD .GT. 5) GO TO 703
ENDIF
IWORK(LMXORD)=MXORD
C
C Set and/or check inputs for constraint checking (INFO(10) .NE. 0).
C Set values for ICNFLG, NONNEG, and pointer LID.
C
ICNFLG = 0
NONNEG = 0
LID = LICNS
IF (INFO(10) .EQ. 0) GO TO 20
IF (INFO(10) .EQ. 1) THEN
ICNFLG = 1
NONNEG = 0
LID = LICNS + NEQ
ELSEIF (INFO(10) .EQ. 2) THEN
ICNFLG = 0
NONNEG = 1
ELSE
ICNFLG = 1
NONNEG = 1
LID = LICNS + NEQ
ENDIF
C
20 CONTINUE
C
C Set and/or check inputs for Krylov solver (INFO(12) .NE. 0).
C If indicated, set default values for MAXL, KMP, NRMAX, and EPLI.
C Otherwise, verify inputs required for iterative solver.
C
IF (INFO(12) .EQ. 0) GO TO 25
C
IWORK(LMITER) = INFO(12)
IF (INFO(13) .EQ. 0) THEN
IWORK(LMAXL) = MIN(5,NEQ)
IWORK(LKMP) = IWORK(LMAXL)
IWORK(LNRMAX) = 5
RWORK(LEPLI) = 0.05D0
ELSE
IF(IWORK(LMAXL) .LT. 1 .OR. IWORK(LMAXL) .GT. NEQ) GO TO 720
IF(IWORK(LKMP) .LT. 1 .OR. IWORK(LKMP) .GT. IWORK(LMAXL))
1 GO TO 721
IF(IWORK(LNRMAX) .LT. 0) GO TO 722
IF(RWORK(LEPLI).LE.0.0D0 .OR. RWORK(LEPLI).GE.1.0D0)GO TO 723
ENDIF
C
25 CONTINUE
C
C Set and/or check controls for the initial condition calculation
C (INFO(11) .GT. 0). If indicated, set default values.
C Otherwise, verify inputs required for iterative solver.
C
IF (INFO(11) .EQ. 0) GO TO 30
IF (INFO(17) .EQ. 0) THEN
IWORK(LMXNIT) = 5
IF (INFO(12) .GT. 0) IWORK(LMXNIT) = 15
IWORK(LMXNJ) = 6
IF (INFO(12) .GT. 0) IWORK(LMXNJ) = 2
IWORK(LMXNH) = 5
IWORK(LLSOFF) = 0
RWORK(LEPIN) = 0.01D0
ELSE
IF (IWORK(LMXNIT) .LE. 0) GO TO 725
IF (IWORK(LMXNJ) .LE. 0) GO TO 725
IF (IWORK(LMXNH) .LE. 0) GO TO 725
LSOFF = IWORK(LLSOFF)
IF (LSOFF .LT. 0 .OR. LSOFF .GT. 1) GO TO 725
IF (RWORK(LEPIN) .LE. 0.0D0) GO TO 725
ENDIF
C
30 CONTINUE
C
C Below is the computation and checking of the work array lengths
C LENIW and LENRW, using direct methods (INFO(12) = 0) or
C the Krylov methods (INFO(12) = 1).
C
LENIC = 0
IF (INFO(10) .EQ. 1 .OR. INFO(10) .EQ. 3) LENIC = NEQ
LENID = 0
IF (INFO(11) .EQ. 1 .OR. INFO(16) .EQ. 1) LENID = NEQ
IF (INFO(12) .EQ. 0) THEN
C
C Compute MTYPE, etc. Check ML and MU.
C
NCPHI = MAX(MXORD + 1, 4)
IF(INFO(6).EQ.0) THEN
LENPD = NEQ**2
LENRW = 50 + (NCPHI+3)*NEQ + LENPD
IF(INFO(5).EQ.0) THEN
IWORK(LMTYPE)=2
ELSE
IWORK(LMTYPE)=1
ENDIF
ELSE
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).EQ.0) THEN
IWORK(LMTYPE)=5
MBAND=IWORK(LML)+IWORK(LMU)+1
MSAVE=(NEQ/MBAND)+1
LENRW = 50 + (NCPHI+3)*NEQ + LENPD + 2*MSAVE
ELSE
IWORK(LMTYPE)=4
LENRW = 50 + (NCPHI+3)*NEQ + LENPD
ENDIF
ENDIF
C
C Compute LENIW, LENWP, LENIWP.
C
LENIW = 40 + LENIC + LENID + NEQ
LENWP = 0
LENIWP = 0
C
ELSE IF (INFO(12) .EQ. 1) THEN
NCPHI = MXORD + 1
MAXL = IWORK(LMAXL)
LENWP = IWORK(LLNWP)
LENIWP = IWORK(LLNIWP)
LENPD = (MAXL+3+MIN0(1,MAXL-IWORK(LKMP)))*NEQ
1 + (MAXL+3)*MAXL + 1 + LENWP
LENRW = 50 + (MXORD+5)*NEQ + LENPD
LENIW = 40 + LENIC + LENID + LENIWP
C
ENDIF
IF(INFO(16) .NE. 0) LENRW = LENRW + NEQ
C
C Check lengths of RWORK and IWORK.
C
IWORK(LNIW)=LENIW
IWORK(LNRW)=LENRW
IWORK(LNPD)=LENPD
IWORK(LLOCWP) = LENPD-LENWP+1
IF(LRW.LT.LENRW)GO TO 704
IF(LIW.LT.LENIW)GO TO 705
C
C Check ICNSTR for legality.
C
IF (LENIC .GT. 0) THEN
DO 40 I = 1,NEQ
ICI = IWORK(LICNS-1+I)
IF (ICI .LT. -2 .OR. ICI .GT. 2) GO TO 726
40 CONTINUE
ENDIF
C
C Check Y for consistency with constraints.
C
IF (LENIC .GT. 0) THEN
CALL DCNST0(NEQ,Y,IWORK(LICNS),IRET)
IF (IRET .NE. 0) GO TO 727
ENDIF
C
C Check ID for legality and set INDEX = 0 or 1.
C
INDEX = 1
IF (LENID .GT. 0) THEN
INDEX = 0
DO 50 I = 1,NEQ
IDI = IWORK(LID-1+I)
IF (IDI .NE. 1 .AND. IDI .NE. -1) GO TO 724
IF (IDI .EQ. -1) INDEX = 1
50 CONTINUE
ENDIF
C
C Check to see that TOUT is different from T.
C
IF(TOUT .EQ. T)GO TO 719
C
C Check HMAX.
C
IF(INFO(7) .NE. 0) THEN
HMAX = RWORK(LHMAX)
IF (HMAX .LE. 0.0D0) GO TO 710
ENDIF
C
C Initialize counters and other flags.
C
IWORK(LNST)=0
IWORK(LNRE)=0
IWORK(LNJE)=0
IWORK(LETF)=0
IWORK(LNCFN)=0
IWORK(LNNI)=0
IWORK(LNLI)=0
IWORK(LNPS)=0
IWORK(LNCFL)=0
IWORK(LKPRIN)=INFO(18)
IDID=1
GO TO 200
C
C-----------------------------------------------------------------------
C This block is for continuation calls only.
C Here we check INFO(1), and if the last step was interrupted,
C we check whether appropriate action was taken.
C-----------------------------------------------------------------------
C
100 CONTINUE
IF(INFO(1).EQ.1)GO TO 110
ITEMP = 1
IF(INFO(1).NE.-1)GO TO 701
C
C If we are here, the last step was interrupted by an error
C condition from DDSTP, and appropriate action was not taken.
C This is a fatal error.
C
MSG = 'DASPK-- THE LAST STEP TERMINATED WITH A NEGATIVE'
CALL XERRWD(MSG,49,201,0,0,0,0,0,0.0D0,0.0D0)
MSG = 'DASPK-- VALUE (=I1) OF IDID AND NO APPROPRIATE'
CALL XERRWD(MSG,47,202,0,1,IDID,0,0,0.0D0,0.0D0)
MSG = 'DASPK-- ACTION WAS TAKEN. RUN TERMINATED'
CALL XERRWD(MSG,41,203,1,0,0,0,0,0.0D0,0.0D0)
RETURN
110 CONTINUE
C
C-----------------------------------------------------------------------
C This block is executed on all calls.
C
C Counters are saved for later checks of performance.
C Then the error tolerance parameters are checked, and the
C work array pointers are set.
C-----------------------------------------------------------------------
C
200 CONTINUE
C
C Save counters for use later.
C
IWORK(LNSTL)=IWORK(LNST)
NLI0 = IWORK(LNLI)
NNI0 = IWORK(LNNI)
NCFN0 = IWORK(LNCFN)
NCFL0 = IWORK(LNCFL)
NWARN = 0
C
C Check RTOL and ATOL.
C
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 707
210 CONTINUE
IF (NZFLG .EQ. 0) GO TO 708
C
C Set pointers to RWORK and IWORK segments.
C For direct methods, SAVR is not used.
C
IWORK(LLCIWP) = LID + LENID
LSAVR = LDELTA
IF (INFO(12) .NE. 0) LSAVR = LDELTA + NEQ
LE = LSAVR + NEQ
LWT = LE + NEQ
LVT = LWT
IF (INFO(16) .NE. 0) LVT = LWT + NEQ
LPHI = LVT + NEQ
LWM = LPHI + NCPHI*NEQ
IF (INFO(1) .EQ. 1) GO TO 400
C
C-----------------------------------------------------------------------
C This block is executed on the initial call only.
C Set the initial step size, the error weight vector, and PHI.
C Compute unknown initial components of Y and YPRIME, if requested.
C-----------------------------------------------------------------------
C
300 CONTINUE
TN=T
IDID=1
C
C Set error weight array WT and altered weight array VT.
C
CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
CALL DINVWT(NEQ,RWORK(LWT),IER)
IF (IER .NE. 0) GO TO 713
IF (INFO(16) .NE. 0) THEN
DO 305 I = 1, NEQ
305 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1)
ENDIF
C
C Compute unit roundoff and HMIN.
C
UROUND = D1MACH(4)
RWORK(LROUND) = UROUND
HMIN = 4.0D0*UROUND*MAX(ABS(T),ABS(TOUT))
C
C Set/check STPTOL control for initial condition calculation.
C
IF (INFO(11) .NE. 0) THEN
IF( INFO(17) .EQ. 0) THEN
RWORK(LSTOL) = UROUND**.6667D0
ELSE
IF (RWORK(LSTOL) .LE. 0.0D0) GO TO 725
ENDIF
ENDIF
C
C Compute EPCON and square root of NEQ and its reciprocal, used
C inside iterative solver.
C
RWORK(LEPCON) = 0.33D0
FLOATN = NEQ
RWORK(LSQRN) = SQRT(FLOATN)
RWORK(LRSQRN) = 1.D0/RWORK(LSQRN)
C
C Check initial interval to see that it is long enough.
C
TDIST = ABS(TOUT - T)
IF(TDIST .LT. HMIN) GO TO 714
C
C Check H0, if this was input.
C
IF (INFO(8) .EQ. 0) GO TO 310
H0 = RWORK(LH)
IF ((TOUT - T)*H0 .LT. 0.0D0) GO TO 711
IF (H0 .EQ. 0.0D0) GO TO 712
GO TO 320
310 CONTINUE
C
C Compute initial stepsize, to be used by either
C DDSTP or DDASIC, depending on INFO(11).
C
H0 = 0.001D0*TDIST
YPNORM = DDWNRM(NEQ,YPRIME,RWORK(LVT),RPAR,IPAR)
IF (YPNORM .GT. 0.5D0/H0) H0 = 0.5D0/YPNORM
H0 = SIGN(H0,TOUT-T)
C
C Adjust H0 if necessary to meet HMAX bound.
C
320 IF (INFO(7) .EQ. 0) GO TO 330
RH = ABS(H0)/RWORK(LHMAX)
IF (RH .GT. 1.0D0) H0 = H0/RH
C
C Check against TSTOP, if applicable.
C
330 IF (INFO(4) .EQ. 0) GO TO 340
TSTOP = RWORK(LTSTOP)
write(*,*) 'tstop = ',tstop
IF ((TSTOP - T)*H0 .LT. 0.0D0) GO TO 715
IF ((T + H0 - TSTOP)*H0 .GT. 0.0D0) H0 = TSTOP - T
IF ((TSTOP - TOUT)*H0 .LT. 0.0D0) GO TO 709
C
340 IF (INFO(11) .EQ. 0) GO TO 370
C
C Compute unknown components of initial Y and YPRIME, depending
C on INFO(11) and INFO(12). INFO(12) represents the nonlinear
C solver type (direct/Krylov). Pass the name of the specific
C nonlinear solver, depending on INFO(12). The location of the work
C arrays SAVR, YIC, YPIC, PWK also differ in the two cases.
C For use in stopping tests, pass TSCALE = TDIST if INDEX = 0.
C
NWT = 1
EPCONI = RWORK(LEPIN)*RWORK(LEPCON)
TSCALE = 0.0D0
IF (INDEX .EQ. 0) TSCALE = TDIST
350 IF (INFO(12) .EQ. 0) THEN
LYIC = LPHI + 2*NEQ
LYPIC = LYIC + NEQ
LPWK = LYPIC
CALL DDASIC(TN,Y,YPRIME,NEQ,INFO(11),IWORK(LID),
* RES,JAC,PSOL,H0,TSCALE,RWORK(LWT),NWT,IDID,RPAR,IPAR,
* RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE),
* RWORK(LYIC),RWORK(LYPIC),RWORK(LPWK),RWORK(LWM),IWORK(LIWM),
* RWORK(LROUND),RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
* EPCONI,RWORK(LSTOL),INFO(15),ICNFLG,IWORK(LICNS),DDASID)
ELSE IF (INFO(12) .EQ. 1) THEN
LYIC = LWM
LYPIC = LYIC + NEQ
LPWK = LYPIC + NEQ
CALL DDASIC(TN,Y,YPRIME,NEQ,INFO(11),IWORK(LID),
* RES,JAC,PSOL,H0,TSCALE,RWORK(LWT),NWT,IDID,RPAR,IPAR,
* RWORK(LPHI),RWORK(LSAVR),RWORK(LDELTA),RWORK(LE),
* RWORK(LYIC),RWORK(LYPIC),RWORK(LPWK),RWORK(LWM),IWORK(LIWM),
* RWORK(LROUND),RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
* EPCONI,RWORK(LSTOL),INFO(15),ICNFLG,IWORK(LICNS),DDASIK)
ENDIF
C
IF (IDID .LT. 0) GO TO 600
C
C DDASIC was successful. If this was the first call to DDASIC,
C update the WT array (with the current Y) and call it again.
C
IF (NWT .EQ. 2) GO TO 355
NWT = 2
CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
CALL DINVWT(NEQ,RWORK(LWT),IER)
IF (IER .NE. 0) GO TO 713
GO TO 350
C
C If INFO(14) = 1, return now with IDID = 4.
C
355 IF (INFO(14) .EQ. 1) THEN
IDID = 4
H = H0
IF (INFO(11) .EQ. 1) RWORK(LHOLD) = H0
GO TO 590
ENDIF
C
C Update the WT and VT arrays one more time, with the new Y.
C
CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
CALL DINVWT(NEQ,RWORK(LWT),IER)
IF (IER .NE. 0) GO TO 713
IF (INFO(16) .NE. 0) THEN
DO 357 I = 1, NEQ
357 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1)
ENDIF
C
C Reset the initial stepsize to be used by DDSTP.
C Use H0, if this was input. Otherwise, recompute H0,
C and adjust it if necessary to meet HMAX bound.
C
IF (INFO(8) .NE. 0) THEN
H0 = RWORK(LH)
GO TO 360
ENDIF
C
H0 = 0.001D0*TDIST
YPNORM = DDWNRM(NEQ,YPRIME,RWORK(LVT),RPAR,IPAR)
IF (YPNORM .GT. 0.5D0/H0) H0 = 0.5D0/YPNORM
H0 = SIGN(H0,TOUT-T)
C
360 IF (INFO(7) .NE. 0) THEN
RH = ABS(H0)/RWORK(LHMAX)
IF (RH .GT. 1.0D0) H0 = H0/RH
ENDIF
C
C Check against TSTOP, if applicable.
C
IF (INFO(4) .NE. 0) THEN
TSTOP = RWORK(LTSTOP)
write(*,*) 'tstop = ',tstop
IF ((T + H0 - TSTOP)*H0 .GT. 0.0D0) H0 = TSTOP - T
ENDIF
C
C Load H and RWORK(LH) with H0.
C
370 H = H0
RWORK(LH) = H
C
C Load Y and H*YPRIME into PHI(*,1) and PHI(*,2).
C
ITEMP = LPHI + NEQ
DO 380 I = 1,NEQ
RWORK(LPHI + I - 1) = Y(I)
380 RWORK(ITEMP + I - 1) = H*YPRIME(I)
C
GO TO 500
C
C-----------------------------------------------------------------------
C This block is for continuation calls only.
C Its purpose is to check stop conditions before taking a step.
C Adjust H if necessary to meet HMAX bound.
C-----------------------------------------------------------------------
C
400 CONTINUE
UROUND=RWORK(LROUND)
DONE = .FALSE.
TN=RWORK(LTN)
H=RWORK(LH)
IF(INFO(7) .EQ. 0) GO TO 410
RH = ABS(H)/RWORK(LHMAX)
IF(RH .GT. 1.0D0) H = H/RH
410 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 490
420 IF((TN-T)*H .LE. 0.0D0) GO TO 490
IF((TN - TOUT)*H .GE. 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 490
425 CONTINUE
CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
* RWORK(LPHI),RWORK(LPSI))
T = TOUT
IDID = 3
DONE = .TRUE.
GO TO 490
430 IF(INFO(3) .EQ. 1) GO TO 440
TSTOP=RWORK(LTSTOP)
write(*,*) 'tstop = ',tstop
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 490
440 TSTOP = RWORK(LTSTOP)
write(*,*) 'tstop = ',tstop
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 .GE. 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 490
445 CONTINUE
CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
* RWORK(LPHI),RWORK(LPSI))
T = TOUT
IDID = 3
DONE = .TRUE.
GO TO 490
450 CONTINUE
C
C Check whether we are within roundoff of TSTOP.
C
IF(ABS(TN-TSTOP).GT.100.0D0*UROUND*
* (ABS(TN)+ABS(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 490
460 TNEXT=TN+H
IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490
H=TSTOP-TN
RWORK(LH)=H
C
490 IF (DONE) GO TO 590
C
C-----------------------------------------------------------------------
C The next block contains the call to the one-step integrator DDSTP.
C This is a looping point for the integration steps.
C Check for too many steps.
C Check for poor Newton/Krylov performance.
C Update WT. Check for too much accuracy requested.
C Compute minimum stepsize.
C-----------------------------------------------------------------------
C
500 CONTINUE
C
C Check for too many steps.
C
IF((IWORK(LNST)-IWORK(LNSTL)).LT.500) GO TO 505
IDID=-1
GO TO 527
C
C Check for poor Newton/Krylov performance.
C
505 IF (INFO(12) .EQ. 0) GO TO 510
NSTD = IWORK(LNST) - IWORK(LNSTL)
NNID = IWORK(LNNI) - NNI0
IF (NSTD .LT. 10 .OR. NNID .EQ. 0) GO TO 510
AVLIN = REAL(IWORK(LNLI) - NLI0)/REAL(NNID)
RCFN = REAL(IWORK(LNCFN) - NCFN0)/REAL(NSTD)
RCFL = REAL(IWORK(LNCFL) - NCFL0)/REAL(NNID)
FMAXL = IWORK(LMAXL)
LAVL = AVLIN .GT. FMAXL
LCFN = RCFN .GT. 0.9D0
LCFL = RCFL .GT. 0.9D0
LWARN = LAVL .OR. LCFN .OR. LCFL
IF (.NOT.LWARN) GO TO 510
NWARN = NWARN + 1
IF (NWARN .GT. 10) GO TO 510
IF (LAVL) THEN
MSG = 'DASPK-- Warning. Poor iterative algorithm performance '
CALL XERRWD (MSG, 56, 501, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
MSG = ' at T = R1. Average no. of linear iterations = R2 '
CALL XERRWD (MSG, 56, 501, 0, 0, 0, 0, 2, TN, AVLIN)
ENDIF
IF (LCFN) THEN
MSG = 'DASPK-- Warning. Poor iterative algorithm performance '
CALL XERRWD (MSG, 56, 502, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
MSG = ' at T = R1. Nonlinear convergence failure rate = R2'
CALL XERRWD (MSG, 56, 502, 0, 0, 0, 0, 2, TN, RCFN)
ENDIF
IF (LCFL) THEN
MSG = 'DASPK-- Warning. Poor iterative algorithm performance '
CALL XERRWD (MSG, 56, 503, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
MSG = ' at T = R1. Linear convergence failure rate = R2 '
CALL XERRWD (MSG, 56, 503, 0, 0, 0, 0, 2, TN, RCFL)
ENDIF
C
C Update WT and VT, if this is not the first call.
C
510 CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),RWORK(LWT),
* RPAR,IPAR)
CALL DINVWT(NEQ,RWORK(LWT),IER)
IF (IER .NE. 0) THEN
IDID = -3
GO TO 527
ENDIF
IF (INFO(16) .NE. 0) THEN
DO 515 I = 1, NEQ
515 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1)
ENDIF
C
C Test for too much accuracy requested.
C
R = DDWNRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*100.0D0*UROUND
IF (R .LE. 1.0D0) GO TO 525
C
C Multiply RTOL and ATOL by R and return.
C
IF(INFO(2).EQ.1)GO TO 523
RTOL(1)=R*RTOL(1)
ATOL(1)=R*ATOL(1)
IDID=-2
GO TO 527
523 DO 524 I=1,NEQ
RTOL(I)=R*RTOL(I)
524 ATOL(I)=R*ATOL(I)
IDID=-2
GO TO 527
525 CONTINUE
C
C Compute minimum stepsize.
C
HMIN=4.0D0*UROUND*MAX(ABS(TN),ABS(TOUT))
C
C Test H vs. HMAX
IF (INFO(7) .NE. 0) THEN
RH = ABS(H)/RWORK(LHMAX)
IF (RH .GT. 1.0D0) H = H/RH
ENDIF
C
C Call the one-step integrator.
C Note that INFO(12) represents the nonlinear solver type.
C Pass the required nonlinear solver, depending upon INFO(12).
C
IF (INFO(12) .EQ. 0) THEN
CALL DDSTP(TN,Y,YPRIME,NEQ,
* RES,JAC,PSOL,H,RWORK(LWT),RWORK(LVT),INFO(1),IDID,RPAR,IPAR,
* RWORK(LPHI),RWORK(LSAVR),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), RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
* RWORK(LEPCON), IWORK(LPHASE),IWORK(LJCALC),INFO(15),
* IWORK(LK), IWORK(LKOLD),IWORK(LNS),NONNEG,INFO(12),
* DNEDD)
ELSE IF (INFO(12) .EQ. 1) THEN
CALL DDSTP(TN,Y,YPRIME,NEQ,
* RES,JAC,PSOL,H,RWORK(LWT),RWORK(LVT),INFO(1),IDID,RPAR,IPAR,
* RWORK(LPHI),RWORK(LSAVR),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), RWORK(LEPLI),RWORK(LSQRN),RWORK(LRSQRN),
* RWORK(LEPCON), IWORK(LPHASE),IWORK(LJCALC),INFO(15),
* IWORK(LK), IWORK(LKOLD),IWORK(LNS),NONNEG,INFO(12),
* DNEDK)
ENDIF
C
527 IF(IDID.LT.0)GO TO 600
C
C-----------------------------------------------------------------------
C This block handles the case of a successful return from DDSTP
C (IDID=1). Test for stop conditions.
C-----------------------------------------------------------------------
C
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 580
530 IF((TN-TOUT)*H.GE.0.0D0)GO TO 535
T=TN
IDID=1
GO TO 580
535 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
IDID=3
T=TOUT
GO TO 580
540 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 580
542 IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*
* (ABS(TN)+ABS(H)))GO TO 545
TNEXT=TN+H
IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500
H=TSTOP-TN
GO TO 500
545 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
IDID=2
T=TSTOP
GO TO 580
550 IF((TN-TOUT)*H.GE.0.0D0)GO TO 555
IF(ABS(TN-TSTOP).LE.100.0D0*UROUND*(ABS(TN)+ABS(H)))GO TO 552
T=TN
IDID=1
GO TO 580
552 CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
IDID=2
T=TSTOP
GO TO 580
555 CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
* IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
T=TOUT
IDID=3
580 CONTINUE
C
C-----------------------------------------------------------------------
C All successful returns from DDASPK are made from this block.
C-----------------------------------------------------------------------
C
590 CONTINUE
RWORK(LTN)=TN
RWORK(LH)=H
RETURN
C
C-----------------------------------------------------------------------
C This block handles all unsuccessful returns other than for
C illegal input.
C-----------------------------------------------------------------------
C
600 CONTINUE
ITEMP = -IDID
GO TO (610,620,630,700,655,640,650,660,670,675,
* 680,685,690,695), ITEMP
C
C The maximum number of steps was taken before
C reaching tout.
C
610 MSG = 'DASPK-- AT CURRENT T (=R1) 500 STEPS'
CALL XERRWD(MSG,38,610,0,0,0,0,1,TN,0.0D0)
MSG = 'DASPK-- TAKEN ON THIS CALL BEFORE REACHING TOUT'
CALL XERRWD(MSG,48,611,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C Too much accuracy for machine precision.
C
620 MSG = 'DASPK-- AT T (=R1) TOO MUCH ACCURACY REQUESTED'
CALL XERRWD(MSG,47,620,0,0,0,0,1,TN,0.0D0)
MSG = 'DASPK-- FOR PRECISION OF MACHINE. RTOL AND ATOL'
CALL XERRWD(MSG,48,621,0,0,0,0,0,0.0D0,0.0D0)
MSG = 'DASPK-- WERE INCREASED TO APPROPRIATE VALUES'
CALL XERRWD(MSG,45,622,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C WT(I) .LE. 0.0D0 for some I (not at start of problem).
C
630 MSG = 'DASPK-- AT T (=R1) SOME ELEMENT OF WT'
CALL XERRWD(MSG,38,630,0,0,0,0,1,TN,0.0D0)
MSG = 'DASPK-- HAS BECOME .LE. 0.0'
CALL XERRWD(MSG,28,631,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C Error test failed repeatedly or with H=HMIN.
C
640 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE'
CALL XERRWD(MSG,44,640,0,0,0,0,2,TN,H)
MSG='DASPK-- ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN'
CALL XERRWD(MSG,57,641,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C Nonlinear solver failed to converge repeatedly or with H=HMIN.
C
650 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE'
CALL XERRWD(MSG,44,650,0,0,0,0,2,TN,H)
MSG = 'DASPK-- NONLINEAR SOLVER FAILED TO CONVERGE'
CALL XERRWD(MSG,44,651,0,0,0,0,0,0.0D0,0.0D0)
MSG = 'DASPK-- REPEATEDLY OR WITH ABS(H)=HMIN'
CALL XERRWD(MSG,40,652,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C The preconditioner had repeated failures.
C
655 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE'
CALL XERRWD(MSG,44,655,0,0,0,0,2,TN,H)
MSG = 'DASPK-- PRECONDITIONER HAD REPEATED FAILURES.'
CALL XERRWD(MSG,46,656,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C The iteration matrix is singular.
C
660 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE'
CALL XERRWD(MSG,44,660,0,0,0,0,2,TN,H)
MSG = 'DASPK-- ITERATION MATRIX IS SINGULAR.'
CALL XERRWD(MSG,38,661,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C Nonlinear system failure preceded by error test failures.
C
670 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE'
CALL XERRWD(MSG,44,670,0,0,0,0,2,TN,H)
MSG = 'DASPK-- NONLINEAR SOLVER COULD NOT CONVERGE.'
CALL XERRWD(MSG,45,671,0,0,0,0,0,0.0D0,0.0D0)
MSG = 'DASPK-- ALSO, THE ERROR TEST FAILED REPEATEDLY.'
CALL XERRWD(MSG,49,672,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C Nonlinear system failure because IRES = -1.
C
675 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE'
CALL XERRWD(MSG,44,675,0,0,0,0,2,TN,H)
MSG = 'DASPK-- NONLINEAR SYSTEM SOLVER COULD NOT CONVERGE'
CALL XERRWD(MSG,51,676,0,0,0,0,0,0.0D0,0.0D0)
MSG = 'DASPK-- BECAUSE IRES WAS EQUAL TO MINUS ONE'
CALL XERRWD(MSG,44,677,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C Failure because IRES = -2.
C
680 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2)'
CALL XERRWD(MSG,40,680,0,0,0,0,2,TN,H)
MSG = 'DASPK-- IRES WAS EQUAL TO MINUS TWO'
CALL XERRWD(MSG,36,681,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C Failed to compute initial YPRIME.
C
685 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE'
CALL XERRWD(MSG,44,685,0,0,0,0,0,0.0D0,0.0D0)
MSG = 'DASPK-- INITIAL (Y,YPRIME) COULD NOT BE COMPUTED'
CALL XERRWD(MSG,49,686,0,0,0,0,2,TN,H0)
GO TO 700
C
C Failure because IER was negative from PSOL.
C
690 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2)'
CALL XERRWD(MSG,40,690,0,0,0,0,2,TN,H)
MSG = 'DASPK-- IER WAS NEGATIVE FROM PSOL'
CALL XERRWD(MSG,35,691,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C Failure because the linear system solver could not converge.
C
695 MSG = 'DASPK-- AT T (=R1) AND STEPSIZE H (=R2) THE'
CALL XERRWD(MSG,44,695,0,0,0,0,2,TN,H)
MSG = 'DASPK-- LINEAR SYSTEM SOLVER COULD NOT CONVERGE.'
CALL XERRWD(MSG,50,696,0,0,0,0,0,0.0D0,0.0D0)
GO TO 700
C
C
700 CONTINUE
INFO(1)=-1
T=TN
RWORK(LTN)=TN
RWORK(LH)=H
RETURN
C
C-----------------------------------------------------------------------
C This block handles all error returns due to illegal input,
C as detected before calling DDSTP.
C First the error message routine is called. If this happens
C twice in succession, execution is terminated.
C-----------------------------------------------------------------------
C
701 MSG = 'DASPK-- ELEMENT (=I1) OF INFO VECTOR IS NOT VALID'
CALL XERRWD(MSG,50,1,0,1,ITEMP,0,0,0.0D0,0.0D0)
GO TO 750
702 MSG = 'DASPK-- NEQ (=I1) .LE. 0'
CALL XERRWD(MSG,25,2,0,1,NEQ,0,0,0.0D0,0.0D0)
GO TO 750
703 MSG = 'DASPK-- MAXORD (=I1) NOT IN RANGE'
CALL XERRWD(MSG,34,3,0,1,MXORD,0,0,0.0D0,0.0D0)
GO TO 750
704 MSG='DASPK-- RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)'
CALL XERRWD(MSG,60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0)
GO TO 750
705 MSG='DASPK-- IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)'
CALL XERRWD(MSG,60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0)
GO TO 750
706 MSG = 'DASPK-- SOME ELEMENT OF RTOL IS .LT. 0'
CALL XERRWD(MSG,39,6,0,0,0,0,0,0.0D0,0.0D0)
GO TO 750
707 MSG = 'DASPK-- SOME ELEMENT OF ATOL IS .LT. 0'
CALL XERRWD(MSG,39,7,0,0,0,0,0,0.0D0,0.0D0)
GO TO 750
708 MSG = 'DASPK-- ALL ELEMENTS OF RTOL AND ATOL ARE ZERO'
CALL XERRWD(MSG,47,8,0,0,0,0,0,0.0D0,0.0D0)
GO TO 750
709 MSG='DASPK-- INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)'
CALL XERRWD(MSG,54,9,0,0,0,0,2,TSTOP,TOUT)
GO TO 750
710 MSG = 'DASPK-- HMAX (=R1) .LT. 0.0'
CALL XERRWD(MSG,28,10,0,0,0,0,1,HMAX,0.0D0)
GO TO 750
711 MSG = 'DASPK-- TOUT (=R1) BEHIND T (=R2)'
CALL XERRWD(MSG,34,11,0,0,0,0,2,TOUT,T)
GO TO 750
712 MSG = 'DASPK-- INFO(8)=1 AND H0=0.0'
CALL XERRWD(MSG,29,12,0,0,0,0,0,0.0D0,0.0D0)
GO TO 750
713 MSG = 'DASPK-- SOME ELEMENT OF WT IS .LE. 0.0'
CALL XERRWD(MSG,39,13,0,0,0,0,0,0.0D0,0.0D0)
GO TO 750
714 MSG='DASPK-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION'
CALL XERRWD(MSG,60,14,0,0,0,0,2,TOUT,T)
GO TO 750
715 MSG = 'DASPK-- INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)'
CALL XERRWD(MSG,49,15,0,0,0,0,2,TSTOP,T)
GO TO 750
717 MSG = 'DASPK-- 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 750
718 MSG = 'DASPK-- 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 750
719 MSG = 'DASPK-- TOUT (=R1) IS EQUAL TO T (=R2)'
CALL XERRWD(MSG,39,19,0,0,0,0,2,TOUT,T)
GO TO 750
720 MSG = 'DASPK-- MAXL (=I1) ILLEGAL. EITHER .LT. 1 OR .GT. NEQ'
CALL XERRWD(MSG,54,20,0,1,IWORK(LMAXL),0,0,0.0D0,0.0D0)
GO TO 750
721 MSG = 'DASPK-- KMP (=I1) ILLEGAL. EITHER .LT. 1 OR .GT. MAXL'
CALL XERRWD(MSG,54,21,0,1,IWORK(LKMP),0,0,0.0D0,0.0D0)
GO TO 750
722 MSG = 'DASPK-- NRMAX (=I1) ILLEGAL. .LT. 0'
CALL XERRWD(MSG,36,22,0,1,IWORK(LNRMAX),0,0,0.0D0,0.0D0)
GO TO 750
723 MSG = 'DASPK-- EPLI (=R1) ILLEGAL. EITHER .LE. 0.D0 OR .GE. 1.D0'
CALL XERRWD(MSG,58,23,0,0,0,0,1,RWORK(LEPLI),0.0D0)
GO TO 750
724 MSG = 'DASPK-- ILLEGAL IWORK VALUE FOR INFO(11) .NE. 0'
CALL XERRWD(MSG,48,24,0,0,0,0,0,0.0D0,0.0D0)
GO TO 750
725 MSG = 'DASPK-- ONE OF THE INPUTS FOR INFO(17) = 1 IS ILLEGAL'
CALL XERRWD(MSG,54,25,0,0,0,0,0,0.0D0,0.0D0)
GO TO 750
726 MSG = 'DASPK-- ILLEGAL IWORK VALUE FOR INFO(10) .NE. 0'
CALL XERRWD(MSG,48,26,0,0,0,0,0,0.0D0,0.0D0)
GO TO 750
727 MSG = 'DASPK-- Y(I) AND IWORK(40+I) (I=I1) INCONSISTENT'
CALL XERRWD(MSG,49,27,0,1,IRET,0,0,0.0D0,0.0D0)
GO TO 750
750 IF(INFO(1).EQ.-1) GO TO 760
INFO(1)=-1
IDID=-33
RETURN
760 MSG = 'DASPK-- REPEATED OCCURRENCES OF ILLEGAL INPUT'
CALL XERRWD(MSG,46,701,0,0,0,0,0,0.0D0,0.0D0)
770 MSG = 'DASPK-- RUN TERMINATED. APPARENT INFINITE LOOP'
CALL XERRWD(MSG,47,702,1,0,0,0,0,0.0D0,0.0D0)
RETURN
C
C------END OF SUBROUTINE DDASPK-----------------------------------------
END
SUBROUTINE DDASIC (X, Y, YPRIME, NEQ, ICOPT, ID, RES, JAC, PSOL,
* H, TSCALE, WT, NIC, IDID, RPAR, IPAR, PHI, SAVR, DELTA, E,
* YIC, YPIC, PWK, WM, IWM, UROUND, EPLI, SQRTN, RSQRTN,
* EPCONI, STPTOL, JFLG, ICNFLG, ICNSTR, NLSIC)
C
C***BEGIN PROLOGUE DDASIC
C***REFER TO DDASPK
C***DATE WRITTEN 940628 (YYMMDD)
C***REVISION DATE 941206 (YYMMDD)
C***REVISION DATE 950714 (YYMMDD)
C***REVISION DATE 000628 TSCALE argument added.
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DDASIC is a driver routine to compute consistent initial values
C for Y and YPRIME. There are two different options:
C Denoting the differential variables in Y by Y_d, and
C the algebraic variables by Y_a, the problem solved is either:
C 1. Given Y_d, calculate Y_a and Y_d', or
C 2. Given Y', calculate Y.
C In either case, initial values for the given components
C are input, and initial guesses for the unknown components
C must also be provided as input.
C
C The external routine NLSIC solves the resulting nonlinear system.
C
C The parameters represent
C
C X -- Independent variable.
C Y -- Solution vector at X.
C YPRIME -- Derivative of solution vector.
C NEQ -- Number of equations to be integrated.
C ICOPT -- Flag indicating initial condition option chosen.
C ICOPT = 1 for option 1 above.
C ICOPT = 2 for option 2.
C ID -- Array of dimension NEQ, which must be initialized
C if option 1 is chosen.
C ID(i) = +1 if Y_i is a differential variable,
C ID(i) = -1 if Y_i is an algebraic variable.
C RES -- External user-supplied subroutine to evaluate the
C residual. See RES description in DDASPK prologue.
C JAC -- External user-supplied routine to update Jacobian
C or preconditioner information in the nonlinear solver
C (optional). See JAC description in DDASPK prologue.
C PSOL -- External user-supplied routine to solve
C a linear system using preconditioning.
C See PSOL in DDASPK prologue.
C H -- Scaling factor in iteration matrix. DDASIC may
C reduce H to achieve convergence.
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
C WT -- Vector of weights for error criterion.
C NIC -- Input number of initial condition calculation call
C (= 1 or 2).
C IDID -- Completion code. See IDID in DDASPK prologue.
C RPAR,IPAR -- Real and integer parameter arrays that
C are used for communication between the
C calling program and external user routines.
C They are not altered by DNSK
C PHI -- Work space for DDASIC of length at least 2*NEQ.
C SAVR -- Work vector for DDASIC of length NEQ.
C DELTA -- Work vector for DDASIC of length NEQ.
C E -- Work vector for DDASIC of length NEQ.
C YIC,YPIC -- Work vectors for DDASIC, each of length NEQ.
C PWK -- Work vector for DDASIC of length NEQ.
C WM,IWM -- Real and integer arrays storing
C information required by the linear solver.
C EPCONI -- Test constant for Newton iteration convergence.
C ICNFLG -- Flag showing whether constraints on Y are to apply.
C ICNSTR -- Integer array of length NEQ with constraint types.
C
C The other parameters are for use internally by DDASIC.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C DCOPY, NLSIC
C
C***END PROLOGUE DDASIC
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),ID(*),WT(*),PHI(NEQ,*)
DIMENSION SAVR(*),DELTA(*),E(*),YIC(*),YPIC(*),PWK(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*), ICNSTR(*)
EXTERNAL RES, JAC, PSOL, NLSIC
C
PARAMETER (LCFN=15)
PARAMETER (LMXNH=34)
C
C The following parameters are data-loaded here:
C RHCUT = factor by which H is reduced on retry of Newton solve.
C RATEMX = maximum convergence rate for which Newton iteration
C is considered converging.
C
SAVE RHCUT, RATEMX
DATA RHCUT/0.1D0/, RATEMX/0.8D0/
C
C
C-----------------------------------------------------------------------
C BLOCK 1.
C Initializations.
C JSKIP is a flag set to 1 when NIC = 2 and NH = 1, to signal that
C the initial call to the JAC routine is to be skipped then.
C Save Y and YPRIME in PHI. Initialize IDID, NH, and CJ.
C-----------------------------------------------------------------------
C
MXNH = IWM(LMXNH)
IDID = 1
NH = 1
JSKIP = 0
IF (NIC .EQ. 2) JSKIP = 1
CALL DCOPY (NEQ, Y, 1, PHI(1,1), 1)
CALL DCOPY (NEQ, YPRIME, 1, PHI(1,2), 1)
C
IF (ICOPT .EQ. 2) THEN
CJ = 0.0D0
ELSE
CJ = 1.0D0/H
ENDIF
C
C-----------------------------------------------------------------------
C BLOCK 2
C Call the nonlinear system solver to obtain
C consistent initial values for Y and YPRIME.
C-----------------------------------------------------------------------
C
200 CONTINUE
CALL NLSIC(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JAC,PSOL,H,TSCALE,WT,
* JSKIP,RPAR,IPAR,SAVR,DELTA,E,YIC,YPIC,PWK,WM,IWM,CJ,UROUND,
* EPLI,SQRTN,RSQRTN,EPCONI,RATEMX,STPTOL,JFLG,ICNFLG,ICNSTR,
* IERNLS)
C
IF (IERNLS .EQ. 0) RETURN
C
C-----------------------------------------------------------------------
C BLOCK 3
C The nonlinear solver was unsuccessful. Increment NCFN.
C Return with IDID = -12 if either
C IERNLS = -1: error is considered unrecoverable,
C ICOPT = 2: we are doing initialization problem type 2, or
C NH = MXNH: the maximum number of H values has been tried.
C Otherwise (problem 1 with IERNLS .GE. 1), reduce H and try again.
C If IERNLS > 1, restore Y and YPRIME to their original values.
C-----------------------------------------------------------------------
C
IWM(LCFN) = IWM(LCFN) + 1
JSKIP = 0
C
IF (IERNLS .EQ. -1) GO TO 350
IF (ICOPT .EQ. 2) GO TO 350
IF (NH .EQ. MXNH) GO TO 350
C
NH = NH + 1
H = H*RHCUT
CJ = 1.0D0/H
C
IF (IERNLS .EQ. 1) GO TO 200
C
CALL DCOPY (NEQ, PHI(1,1), 1, Y, 1)
CALL DCOPY (NEQ, PHI(1,2), 1, YPRIME, 1)
GO TO 200
C
350 IDID = -12
RETURN
C
C------END OF SUBROUTINE DDASIC-----------------------------------------
END
SUBROUTINE DYYPNW (NEQ, Y, YPRIME, CJ, RL, P, ICOPT, ID,
* YNEW, YPNEW)
C
C***BEGIN PROLOGUE DYYPNW
C***REFER TO DLINSK
C***DATE WRITTEN 940830 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DYYPNW calculates the new (Y,YPRIME) pair needed in the
C linesearch algorithm based on the current lambda value. It is
C called by DLINSK and DLINSD. Based on the ICOPT and ID values,
C the corresponding entry in Y or YPRIME is updated.
C
C In addition to the parameters described in the calling programs,
C the parameters represent
C
C P -- Array of length NEQ that contains the current
C approximate Newton step.
C RL -- Scalar containing the current lambda value.
C YNEW -- Array of length NEQ containing the updated Y vector.
C YPNEW -- Array of length NEQ containing the updated YPRIME
C vector.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED (NONE)
C
C***END PROLOGUE DYYPNW
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION Y(*), YPRIME(*), YNEW(*), YPNEW(*), ID(*), P(*)
C
IF (ICOPT .EQ. 1) THEN
DO 10 I=1,NEQ
IF(ID(I) .LT. 0) THEN
YNEW(I) = Y(I) - RL*P(I)
YPNEW(I) = YPRIME(I)
ELSE
YNEW(I) = Y(I)
YPNEW(I) = YPRIME(I) - RL*CJ*P(I)
ENDIF
10 CONTINUE
ELSE
DO 20 I = 1,NEQ
YNEW(I) = Y(I) - RL*P(I)
YPNEW(I) = YPRIME(I)
20 CONTINUE
ENDIF
RETURN
C----------------------- END OF SUBROUTINE DYYPNW ----------------------
END
SUBROUTINE DDSTP(X,Y,YPRIME,NEQ,RES,JAC,PSOL,H,WT,VT,
* JSTART,IDID,RPAR,IPAR,PHI,SAVR,DELTA,E,WM,IWM,
* ALPHA,BETA,GAMMA,PSI,SIGMA,CJ,CJOLD,HOLD,S,HMIN,UROUND,
* EPLI,SQRTN,RSQRTN,EPCON,IPHASE,JCALC,JFLG,K,KOLD,NS,NONNEG,
* NTYPE,NLS)
C
C***BEGIN PROLOGUE DDSTP
C***REFER TO DDASPK
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***REVISION DATE 940909 (YYMMDD) (Reset PSI(1), PHI(*,2) at 690)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DDSTP solves a system of differential/algebraic equations of
C the form G(X,Y,YPRIME) = 0, for one step (normally from X to X+H).
C
C The methods used are modified divided difference, fixed leading
C coefficient forms of backward differentiation formulas.
C The code adjusts the stepsize and order to control the local error
C per step.
C
C
C The parameters represent
C X -- Independent variable.
C Y -- Solution vector at X.
C YPRIME -- Derivative of solution vector
C after successful step.
C NEQ -- Number of equations to be integrated.
C RES -- External user-supplied subroutine
C to evaluate the residual. See RES description
C in DDASPK prologue.
C JAC -- External user-supplied routine to update
C Jacobian or preconditioner information in the
C nonlinear solver. See JAC description in DDASPK
C prologue.
C PSOL -- External user-supplied routine to solve
C a linear system using preconditioning.
C (This is optional). See PSOL in DDASPK prologue.
C H -- Appropriate step size for next step.
C Normally determined by the code.
C WT -- Vector of weights for error criterion used in Newton test.
C VT -- Masked vector of weights used in error test.
C JSTART -- Integer variable set 0 for
C first step, 1 otherwise.
C IDID -- Completion code returned from the nonlinear solver.
C See IDID description in DDASPK prologue.
C RPAR,IPAR -- Real and integer parameter arrays that
C are used for communication between the
C calling program and external user routines.
C They are not altered by DNSK
C PHI -- Array of divided differences used by
C DDSTP. The length is NEQ*(K+1), where
C K is the maximum order.
C SAVR -- Work vector for DDSTP of length NEQ.
C DELTA,E -- Work vectors for DDSTP of length NEQ.
C WM,IWM -- Real and integer arrays storing
C information required by the linear solver.
C
C The other parameters are information
C which is needed internally by DDSTP to
C continue from step to step.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C NLS, DDWNRM, DDATRP
C
C***END PROLOGUE DDSTP
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),WT(*),VT(*)
DIMENSION PHI(NEQ,*),SAVR(*),DELTA(*),E(*)
DIMENSION WM(*),IWM(*)
DIMENSION PSI(*),ALPHA(*),BETA(*),GAMMA(*),SIGMA(*)
DIMENSION RPAR(*),IPAR(*)
EXTERNAL RES, JAC, PSOL, NLS
C
PARAMETER (LMXORD=3)
PARAMETER (LNST=11, LETF=14, LCFN=15)
C
C
C-----------------------------------------------------------------------
C BLOCK 1.
C Initialize. On the first call, set
C the order to 1 and initialize
C other variables.
C-----------------------------------------------------------------------
C
C Initializations for all calls
C
XOLD=X
NCF=0
NEF=0
IF(JSTART .NE. 0) GO TO 120
C
C If this is the first step, perform
C other initializations
C
K=1
KOLD=0
HOLD=0.0D0
PSI(1)=H
CJ = 1.D0/H
IPHASE = 0
NS=0
120 CONTINUE
C
C
C
C
C
C-----------------------------------------------------------------------
C BLOCK 2
C Compute coefficients of formulas for
C this step.
C-----------------------------------------------------------------------
200 CONTINUE
KP1=K+1
KP2=K+2
KM1=K-1
IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0
NS=MIN0(NS+1,KOLD+2)
NSP1=NS+1
IF(KP1 .LT. NS)GO TO 230
C
BETA(1)=1.0D0
ALPHA(1)=1.0D0
TEMP1=H
GAMMA(1)=0.0D0
SIGMA(1)=1.0D0
DO 210 I=2,KP1
TEMP2=PSI(I-1)
PSI(I-1)=TEMP1
BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2
TEMP1=TEMP2+H
ALPHA(I)=H/TEMP1
SIGMA(I)=(I-1)*SIGMA(I-1)*ALPHA(I)
GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H
210 CONTINUE
PSI(KP1)=TEMP1
230 CONTINUE
C
C Compute ALPHAS, ALPHA0
C
ALPHAS = 0.0D0
ALPHA0 = 0.0D0
DO 240 I = 1,K
ALPHAS = ALPHAS - 1.0D0/I
ALPHA0 = ALPHA0 - ALPHA(I)
240 CONTINUE
C
C Compute leading coefficient CJ
C
CJLAST = CJ
CJ = -ALPHAS/H
C
C Compute variable stepsize error coefficient CK
C
CK = ABS(ALPHA(KP1) + ALPHAS - ALPHA0)
CK = MAX(CK,ALPHA(KP1))
C
C Change PHI to PHI STAR
C
IF(KP1 .LT. NSP1) GO TO 280
DO 270 J=NSP1,KP1
DO 260 I=1,NEQ
260 PHI(I,J)=BETA(J)*PHI(I,J)
270 CONTINUE
280 CONTINUE
C
C Update time
C
X=X+H
C
C Initialize IDID to 1
C
IDID = 1
C
C
C
C
C
C-----------------------------------------------------------------------
C BLOCK 3
C Call the nonlinear system solver to obtain the solution and
C derivative.
C-----------------------------------------------------------------------
C
CALL NLS(X,Y,YPRIME,NEQ,
* RES,JAC,PSOL,H,WT,JSTART,IDID,RPAR,IPAR,PHI,GAMMA,
* SAVR,DELTA,E,WM,IWM,CJ,CJOLD,CJLAST,S,
* UROUND,EPLI,SQRTN,RSQRTN,EPCON,JCALC,JFLG,KP1,
* NONNEG,NTYPE,IERNLS)
C
IF(IERNLS .NE. 0)GO TO 600
C
C
C
C
C
C-----------------------------------------------------------------------
C BLOCK 4
C Estimate the errors at orders K,K-1,K-2
C as if constant stepsize was used. Estimate
C the local error at order K and test
C whether the current step is successful.
C-----------------------------------------------------------------------
C
C Estimate errors at orders K,K-1,K-2
C
ENORM = DDWNRM(NEQ,E,VT,RPAR,IPAR)
ERK = SIGMA(K+1)*ENORM
TERK = (K+1)*ERK
EST = ERK
KNEW=K
IF(K .EQ. 1)GO TO 430
DO 405 I = 1,NEQ
405 DELTA(I) = PHI(I,KP1) + E(I)
ERKM1=SIGMA(K)*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR)
TERKM1 = K*ERKM1
IF(K .GT. 2)GO TO 410
IF(TERKM1 .LE. 0.5*TERK)GO TO 420
GO TO 430
410 CONTINUE
DO 415 I = 1,NEQ
415 DELTA(I) = PHI(I,K) + DELTA(I)
ERKM2=SIGMA(K-1)*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR)
TERKM2 = (K-1)*ERKM2
IF(MAX(TERKM1,TERKM2).GT.TERK)GO TO 430
C
C Lower the order
C
420 CONTINUE
KNEW=K-1
EST = ERKM1
C
C
C Calculate the local error for the current step
C to see if the step was successful
C
430 CONTINUE
ERR = CK * ENORM
IF(ERR .GT. 1.0D0)GO TO 600
C
C
C
C
C
C-----------------------------------------------------------------------
C BLOCK 5
C The step is successful. Determine
C the best order and stepsize for
C the next step. Update the differences
C for the next step.
C-----------------------------------------------------------------------
IDID=1
IWM(LNST)=IWM(LNST)+1
KDIFF=K-KOLD
KOLD=K
HOLD=H
C
C
C Estimate the error at order K+1 unless
C already decided to lower order, or
C already using maximum order, or
C stepsize not constant, or
C order raised in previous step
C
IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1
IF(IPHASE .EQ. 0)GO TO 545
IF(KNEW.EQ.KM1)GO TO 540
IF(K.EQ.IWM(LMXORD)) GO TO 550
IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550
DO 510 I=1,NEQ
510 DELTA(I)=E(I)-PHI(I,KP2)
ERKP1 = (1.0D0/(K+2))*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR)
TERKP1 = (K+2)*ERKP1
IF(K.GT.1)GO TO 520
IF(TERKP1.GE.0.5D0*TERK)GO TO 550
GO TO 530
520 IF(TERKM1.LE.MIN(TERK,TERKP1))GO TO 540
IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550
C
C Raise order
C
530 K=KP1
EST = ERKP1
GO TO 550
C
C Lower order
C
540 K=KM1
EST = ERKM1
GO TO 550
C
C If IPHASE = 0, increase order by one and multiply stepsize by
C factor two
C
545 K = KP1
HNEW = H*2.0D0
H = HNEW
GO TO 575
C
C
C Determine the appropriate stepsize for
C the next step.
C
550 HNEW=H
TEMP2=K+1
R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
IF(R .LT. 2.0D0) GO TO 555
HNEW = 2.0D0*H
GO TO 560
555 IF(R .GT. 1.0D0) GO TO 560
R = MAX(0.5D0,MIN(0.9D0,R))
HNEW = H*R
560 H=HNEW
C
C
C Update differences for next step
C
575 CONTINUE
IF(KOLD.EQ.IWM(LMXORD))GO TO 585
DO 580 I=1,NEQ
580 PHI(I,KP2)=E(I)
585 CONTINUE
DO 590 I=1,NEQ
590 PHI(I,KP1)=PHI(I,KP1)+E(I)
DO 595 J1=2,KP1
J=KP1-J1+1
DO 595 I=1,NEQ
595 PHI(I,J)=PHI(I,J)+PHI(I,J+1)
JSTART = 1
RETURN
C
C
C
C
C
C-----------------------------------------------------------------------
C BLOCK 6
C The step is unsuccessful. Restore X,PSI,PHI
C Determine appropriate stepsize for
C continuing the integration, or exit with
C an error flag if there have been many
C failures.
C-----------------------------------------------------------------------
600 IPHASE = 1
C
C Restore X,PHI,PSI
C
X=XOLD
IF(KP1.LT.NSP1)GO TO 630
DO 620 J=NSP1,KP1
TEMP1=1.0D0/BETA(J)
DO 610 I=1,NEQ
610 PHI(I,J)=TEMP1*PHI(I,J)
620 CONTINUE
630 CONTINUE
DO 640 I=2,KP1
640 PSI(I-1)=PSI(I)-H
C
C
C Test whether failure is due to nonlinear solver
C or error test
C
IF(IERNLS .EQ. 0)GO TO 660
IWM(LCFN)=IWM(LCFN)+1
C
C
C The nonlinear solver failed to converge.
C Determine the cause of the failure and take appropriate action.
C If IERNLS .LT. 0, then return. Otherwise, reduce the stepsize
C and try again, unless too many failures have occurred.
C
IF (IERNLS .LT. 0) GO TO 675
NCF = NCF + 1
R = 0.25D0
H = H*R
IF (NCF .LT. 10 .AND. ABS(H) .GE. HMIN) GO TO 690
IF (IDID .EQ. 1) IDID = -7
IF (NEF .GE. 3) IDID = -9
GO TO 675
C
C
C The nonlinear solver converged, and the cause
C of the failure was the error estimate
C exceeding the tolerance.
C
660 NEF=NEF+1
IWM(LETF)=IWM(LETF)+1
IF (NEF .GT. 1) GO TO 665
C
C On first error test failure, keep current order or lower
C order by one. Compute new stepsize based on differences
C of the solution.
C
K = KNEW
TEMP2 = K + 1
R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
R = MAX(0.25D0,MIN(0.9D0,R))
H = H*R
IF (ABS(H) .GE. HMIN) GO TO 690
IDID = -6
GO TO 675
C
C On second error test failure, use the current order or
C decrease order by one. Reduce the stepsize by a factor of
C one quarter.
C
665 IF (NEF .GT. 2) GO TO 670
K = KNEW
R = 0.25D0
H = R*H
IF (ABS(H) .GE. HMIN) GO TO 690
IDID = -6
GO TO 675
C
C On third and subsequent error test failures, set the order to
C one, and reduce the stepsize by a factor of one quarter.
C
670 K = 1
R = 0.25D0
H = R*H
IF (ABS(H) .GE. HMIN) GO TO 690
IDID = -6
GO TO 675
C
C
C
C
C For all crashes, restore Y to its last value,
C interpolate to find YPRIME at last X, and return.
C
C Before returning, verify that the user has not set
C IDID to a nonnegative value. If the user has set IDID
C to a nonnegative value, then reset IDID to be -7, indicating
C a failure in the nonlinear system solver.
C
675 CONTINUE
CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI)
JSTART = 1
IF (IDID .GE. 0) IDID = -7
RETURN
C
C
C Go back and try this step again.
C If this is the first step, reset PSI(1) and rescale PHI(*,2).
C
690 IF (KOLD .EQ. 0) THEN
PSI(1) = H
DO 695 I = 1,NEQ
695 PHI(I,2) = R*PHI(I,2)
ENDIF
GO TO 200
C
C------END OF SUBROUTINE DDSTP------------------------------------------
END
SUBROUTINE DCNSTR (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
C
C***BEGIN PROLOGUE DCNSTR
C***DATE WRITTEN 950808 (YYMMDD)
C***REVISION DATE 950814 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This subroutine checks for constraint violations in the proposed
C new approximate solution YNEW.
C If a constraint violation occurs, then a new step length, TAU,
C is calculated, and this value is to be given to the linesearch routine
C to calculate a new approximate solution YNEW.
C
C On entry:
C
C NEQ -- size of the nonlinear system, and the length of arrays
C Y, YNEW and ICNSTR.
C
C Y -- real array containing the current approximate y.
C
C YNEW -- real array containing the new approximate y.
C
C ICNSTR -- INTEGER array of length NEQ containing flags indicating
C which entries in YNEW are to be constrained.
C if ICNSTR(I) = 2, then YNEW(I) must be .GT. 0,
C if ICNSTR(I) = 1, then YNEW(I) must be .GE. 0,
C if ICNSTR(I) = -1, then YNEW(I) must be .LE. 0, while
C if ICNSTR(I) = -2, then YNEW(I) must be .LT. 0, while
C if ICNSTR(I) = 0, then YNEW(I) is not constrained.
C
C RLX -- real scalar restricting update, if ICNSTR(I) = 2 or -2,
C to ABS( (YNEW-Y)/Y ) < FAC2*RLX in component I.
C
C TAU -- the current size of the step length for the linesearch.
C
C On return
C
C TAU -- the adjusted size of the step length if a constraint
C violation occurred (otherwise, it is unchanged). it is
C the step length to give to the linesearch routine.
C
C IRET -- output flag.
C IRET=0 means that YNEW satisfied all constraints.
C IRET=1 means that YNEW failed to satisfy all the
C constraints, and a new linesearch step
C must be computed.
C
C IVAR -- index of variable causing constraint to be violated.
C
C-----------------------------------------------------------------------
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(NEQ), YNEW(NEQ), ICNSTR(NEQ)
SAVE FAC, FAC2, ZERO
DATA FAC /0.6D0/, FAC2 /0.9D0/, ZERO/0.0D0/
C-----------------------------------------------------------------------
C Check constraints for proposed new step YNEW. If a constraint has
C been violated, then calculate a new step length, TAU, to be
C used in the linesearch routine.
C-----------------------------------------------------------------------
IRET = 0
RDYMX = ZERO
IVAR = 0
DO 100 I = 1,NEQ
C
IF (ICNSTR(I) .EQ. 2) THEN
RDY = ABS( (YNEW(I)-Y(I))/Y(I) )
IF (RDY .GT. RDYMX) THEN
RDYMX = RDY
IVAR = I
ENDIF
IF (YNEW(I) .LE. ZERO) THEN
TAU = FAC*TAU
IVAR = I
IRET = 1
RETURN
ENDIF
C
ELSEIF (ICNSTR(I) .EQ. 1) THEN
IF (YNEW(I) .LT. ZERO) THEN
TAU = FAC*TAU
IVAR = I
IRET = 1
RETURN
ENDIF
C
ELSEIF (ICNSTR(I) .EQ. -1) THEN
IF (YNEW(I) .GT. ZERO) THEN
TAU = FAC*TAU
IVAR = I
IRET = 1
RETURN
ENDIF
C
ELSEIF (ICNSTR(I) .EQ. -2) THEN
RDY = ABS( (YNEW(I)-Y(I))/Y(I) )
IF (RDY .GT. RDYMX) THEN
RDYMX = RDY
IVAR = I
ENDIF
IF (YNEW(I) .GE. ZERO) THEN
TAU = FAC*TAU
IVAR = I
IRET = 1
RETURN
ENDIF
C
ENDIF
100 CONTINUE
IF(RDYMX .GE. RLX) THEN
TAU = FAC2*TAU*RLX/RDYMX
IRET = 1
ENDIF
C
RETURN
C----------------------- END OF SUBROUTINE DCNSTR ----------------------
END
SUBROUTINE DCNST0 (NEQ, Y, ICNSTR, IRET)
C
C***BEGIN PROLOGUE DCNST0
C***DATE WRITTEN 950808 (YYMMDD)
C***REVISION DATE 950808 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This subroutine checks for constraint violations in the initial
C approximate solution u.
C
C On entry
C
C NEQ -- size of the nonlinear system, and the length of arrays
C Y and ICNSTR.
C
C Y -- real array containing the initial approximate root.
C
C ICNSTR -- INTEGER array of length NEQ containing flags indicating
C which entries in Y are to be constrained.
C if ICNSTR(I) = 2, then Y(I) must be .GT. 0,
C if ICNSTR(I) = 1, then Y(I) must be .GE. 0,
C if ICNSTR(I) = -1, then Y(I) must be .LE. 0, while
C if ICNSTR(I) = -2, then Y(I) must be .LT. 0, while
C if ICNSTR(I) = 0, then Y(I) is not constrained.
C
C On return
C
C IRET -- output flag.
C IRET=0 means that u satisfied all constraints.
C IRET.NE.0 means that Y(IRET) failed to satisfy its
C constraint.
C
C-----------------------------------------------------------------------
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(NEQ), ICNSTR(NEQ)
SAVE ZERO
DATA ZERO/0.D0/
C-----------------------------------------------------------------------
C Check constraints for initial Y. If a constraint has been violated,
C set IRET = I to signal an error return to calling routine.
C-----------------------------------------------------------------------
IRET = 0
DO 100 I = 1,NEQ
IF (ICNSTR(I) .EQ. 2) THEN
IF (Y(I) .LE. ZERO) THEN
IRET = I
RETURN
ENDIF
ELSEIF (ICNSTR(I) .EQ. 1) THEN
IF (Y(I) .LT. ZERO) THEN
IRET = I
RETURN
ENDIF
ELSEIF (ICNSTR(I) .EQ. -1) THEN
IF (Y(I) .GT. ZERO) THEN
IRET = I
RETURN
ENDIF
ELSEIF (ICNSTR(I) .EQ. -2) THEN
IF (Y(I) .GE. ZERO) THEN
IRET = I
RETURN
ENDIF
ENDIF
100 CONTINUE
RETURN
C----------------------- END OF SUBROUTINE DCNST0 ----------------------
END
SUBROUTINE DDAWTS(NEQ,IWT,RTOL,ATOL,Y,WT,RPAR,IPAR)
C
C***BEGIN PROLOGUE DDAWTS
C***REFER TO DDASPK
C***ROUTINES CALLED (NONE)
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***END PROLOGUE DDAWTS
C-----------------------------------------------------------------------
C This subroutine sets the error weight vector,
C WT, according to WT(I)=RTOL(I)*ABS(Y(I))+ATOL(I),
C I = 1 to NEQ.
C RTOL and ATOL are scalars if IWT = 0,
C and vectors if IWT = 1.
C-----------------------------------------------------------------------
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION RTOL(*),ATOL(*),Y(*),WT(*)
DIMENSION RPAR(*),IPAR(*)
RTOLI=RTOL(1)
ATOLI=ATOL(1)
DO 20 I=1,NEQ
IF (IWT .EQ.0) GO TO 10
RTOLI=RTOL(I)
ATOLI=ATOL(I)
10 WT(I)=RTOLI*ABS(Y(I))+ATOLI
20 CONTINUE
RETURN
C
C------END OF SUBROUTINE DDAWTS-----------------------------------------
END
SUBROUTINE DINVWT(NEQ,WT,IER)
C
C***BEGIN PROLOGUE DINVWT
C***REFER TO DDASPK
C***ROUTINES CALLED (NONE)
C***DATE WRITTEN 950125 (YYMMDD)
C***END PROLOGUE DINVWT
C-----------------------------------------------------------------------
C This subroutine checks the error weight vector WT, of length NEQ,
C for components that are .le. 0, and if none are found, it
C inverts the WT(I) in place. This replaces division operations
C with multiplications in all norm evaluations.
C IER is returned as 0 if all WT(I) were found positive,
C and the first I with WT(I) .le. 0.0 otherwise.
C-----------------------------------------------------------------------
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION WT(*)
C
DO 10 I = 1,NEQ
IF (WT(I) .LE. 0.0D0) GO TO 30
10 CONTINUE
DO 20 I = 1,NEQ
20 WT(I) = 1.0D0/WT(I)
IER = 0
RETURN
C
30 IER = I
RETURN
C
C------END OF SUBROUTINE DINVWT-----------------------------------------
END
SUBROUTINE DDATRP(X,XOUT,YOUT,YPOUT,NEQ,KOLD,PHI,PSI)
C
C***BEGIN PROLOGUE DDATRP
C***REFER TO DDASPK
C***ROUTINES CALLED (NONE)
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***END PROLOGUE DDATRP
C
C-----------------------------------------------------------------------
C The methods in subroutine DDSTP use polynomials
C to approximate the solution. DDATRP approximates the
C solution and its derivative at time XOUT by evaluating
C one of these polynomials, and its derivative, there.
C Information defining this polynomial is passed from
C DDSTP, so DDATRP cannot be used alone.
C
C The parameters are
C
C X The current time in the integration.
C XOUT The time at which the solution is desired.
C YOUT The interpolated approximation to Y at XOUT.
C (This is output.)
C YPOUT The interpolated approximation to YPRIME at XOUT.
C (This is output.)
C NEQ Number of equations.
C KOLD Order used on last successful step.
C PHI Array of scaled divided differences of Y.
C PSI Array of past stepsize history.
C-----------------------------------------------------------------------
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION YOUT(*),YPOUT(*)
DIMENSION PHI(NEQ,*),PSI(*)
KOLDP1=KOLD+1
TEMP1=XOUT-X
DO 10 I=1,NEQ
YOUT(I)=PHI(I,1)
10 YPOUT(I)=0.0D0
C=1.0D0
D=0.0D0
GAMMA=TEMP1/PSI(1)
DO 30 J=2,KOLDP1
D=D*GAMMA+C/PSI(J-1)
C=C*GAMMA
GAMMA=(TEMP1+PSI(J-1))/PSI(J)
DO 20 I=1,NEQ
YOUT(I)=YOUT(I)+C*PHI(I,J)
20 YPOUT(I)=YPOUT(I)+D*PHI(I,J)
30 CONTINUE
RETURN
C
C------END OF SUBROUTINE DDATRP-----------------------------------------
END
DOUBLE PRECISION FUNCTION DDWNRM(NEQ,V,RWT,RPAR,IPAR)
C
C***BEGIN PROLOGUE DDWNRM
C***ROUTINES CALLED (NONE)
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***END PROLOGUE DDWNRM
C-----------------------------------------------------------------------
C This function routine computes the weighted
C root-mean-square norm of the vector of length
C NEQ contained in the array V, with reciprocal weights
C contained in the array RWT of length NEQ.
C DDWNRM=SQRT((1/NEQ)*SUM(V(I)*RWT(I))**2)
C-----------------------------------------------------------------------
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION V(*),RWT(*)
DIMENSION RPAR(*),IPAR(*)
DDWNRM = 0.0D0
VMAX = 0.0D0
DO 10 I = 1,NEQ
IF(ABS(V(I)*RWT(I)) .GT. VMAX) VMAX = ABS(V(I)*RWT(I))
10 CONTINUE
IF(VMAX .LE. 0.0D0) GO TO 30
SUM = 0.0D0
DO 20 I = 1,NEQ
20 SUM = SUM + ((V(I)*RWT(I))/VMAX)**2
DDWNRM = VMAX*SQRT(SUM/NEQ)
30 CONTINUE
RETURN
C
C------END OF FUNCTION DDWNRM-------------------------------------------
END
SUBROUTINE DDASID(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JACD,PDUM,H,TSCALE,
* WT,JSDUM,RPAR,IPAR,DUMSVR,DELTA,R,YIC,YPIC,DUMPWK,WM,IWM,CJ,
* UROUND,DUME,DUMS,DUMR,EPCON,RATEMX,STPTOL,JFDUM,
* ICNFLG,ICNSTR,IERNLS)
C
C***BEGIN PROLOGUE DDASID
C***REFER TO DDASPK
C***DATE WRITTEN 940701 (YYMMDD)
C***REVISION DATE 950808 (YYMMDD)
C***REVISION DATE 951110 Removed unreachable block 390.
C***REVISION DATE 000628 TSCALE argument added.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C
C DDASID solves a nonlinear system of algebraic equations of the
C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in
C the initial conditions.
C
C The method used is a modified Newton scheme.
C
C The parameters represent
C
C X -- Independent variable.
C Y -- Solution vector.
C YPRIME -- Derivative of solution vector.
C NEQ -- Number of unknowns.
C ICOPT -- Initial condition option chosen (1 or 2).
C ID -- Array of dimension NEQ, which must be initialized
C if ICOPT = 1. See DDASIC.
C RES -- External user-supplied subroutine to evaluate the
C residual. See RES description in DDASPK prologue.
C JACD -- External user-supplied routine to evaluate the
C Jacobian. See JAC description for the case
C INFO(12) = 0 in the DDASPK prologue.
C PDUM -- Dummy argument.
C H -- Scaling factor for this initial condition calc.
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
C WT -- Vector of weights for error criterion.
C JSDUM -- Dummy argument.
C RPAR,IPAR -- Real and integer arrays used for communication
C between the calling program and external user
C routines. They are not altered within DASPK.
C DUMSVR -- Dummy argument.
C DELTA -- Work vector for NLS of length NEQ.
C R -- Work vector for NLS of length NEQ.
C YIC,YPIC -- Work vectors for NLS, each of length NEQ.
C DUMPWK -- Dummy argument.
C WM,IWM -- Real and integer arrays storing matrix information
C such as the matrix of partial derivatives,
C permutation vector, and various other information.
C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
C UROUND -- Unit roundoff.
C DUME -- Dummy argument.
C DUMS -- Dummy argument.
C DUMR -- Dummy argument.
C EPCON -- Tolerance to test for convergence of the Newton
C iteration.
C RATEMX -- Maximum convergence rate for which Newton iteration
C is considered converging.
C JFDUM -- Dummy argument.
C STPTOL -- Tolerance used in calculating the minimum lambda
C value allowed.
C ICNFLG -- Integer scalar. If nonzero, then constraint
C violations in the proposed new approximate solution
C will be checked for, and the maximum step length
C will be adjusted accordingly.
C ICNSTR -- Integer array of length NEQ containing flags for
C checking constraints.
C IERNLS -- Error flag for nonlinear solver.
C 0 ==> nonlinear solver converged.
C 1,2 ==> recoverable error inside nonlinear solver.
C 1 => retry with current Y, YPRIME
C 2 => retry with original Y, YPRIME
C -1 ==> unrecoverable error in nonlinear solver.
C
C All variables with "DUM" in their names are dummy variables
C which are not used in this routine.
C
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C RES, DMATD, DNSID
C
C***END PROLOGUE DDASID
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),ID(*),WT(*),ICNSTR(*)
DIMENSION DELTA(*),R(*),YIC(*),YPIC(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
EXTERNAL RES, JACD
C
PARAMETER (LNRE=12, LNJE=13, LMXNIT=32, LMXNJ=33)
C
C
C Perform initializations.
C
MXNIT = IWM(LMXNIT)
MXNJ = IWM(LMXNJ)
IERNLS = 0
NJ = 0
C
C Call RES to initialize DELTA.
C
IRES = 0
IWM(LNRE) = IWM(LNRE) + 1
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
IF (IRES .LT. 0) GO TO 370
C
C Looping point for updating the Jacobian.
C
300 CONTINUE
C
C Initialize all error flags to zero.
C
IERJ = 0
IRES = 0
IERNEW = 0
C
C Reevaluate the iteration matrix, J = dG/dY + CJ*dG/dYPRIME,
C where G(X,Y,YPRIME) = 0.
C
NJ = NJ + 1
IWM(LNJE)=IWM(LNJE)+1
CALL DMATD(NEQ,X,Y,YPRIME,DELTA,CJ,H,IERJ,WT,R,
* WM,IWM,RES,IRES,UROUND,JACD,RPAR,IPAR)
IF (IRES .LT. 0 .OR. IERJ .NE. 0) GO TO 370
C
C Call the nonlinear Newton solver for up to MXNIT iterations.
C
CALL DNSID(X,Y,YPRIME,NEQ,ICOPT,ID,RES,WT,RPAR,IPAR,DELTA,R,
* YIC,YPIC,WM,IWM,CJ,TSCALE,EPCON,RATEMX,MXNIT,STPTOL,
* ICNFLG,ICNSTR,IERNEW)
C
IF (IERNEW .EQ. 1 .AND. NJ .LT. MXNJ) THEN
C
C MXNIT iterations were done, the convergence rate is < 1,
C and the number of Jacobian evaluations is less than MXNJ.
C Call RES, reevaluate the Jacobian, and try again.
C
IWM(LNRE)=IWM(LNRE)+1
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
IF (IRES .LT. 0) GO TO 370
GO TO 300
ENDIF
C
IF (IERNEW .NE. 0) GO TO 380
RETURN
C
C
C Unsuccessful exits from nonlinear solver.
C Compute IERNLS accordingly.
C
370 IERNLS = 2
IF (IRES .LE. -2) IERNLS = -1
RETURN
C
380 IERNLS = MIN(IERNEW,2)
RETURN
C
C------END OF SUBROUTINE DDASID-----------------------------------------
END
SUBROUTINE DNSID(X,Y,YPRIME,NEQ,ICOPT,ID,RES,WT,RPAR,IPAR,
* DELTA,R,YIC,YPIC,WM,IWM,CJ,TSCALE,EPCON,RATEMX,MAXIT,STPTOL,
* ICNFLG,ICNSTR,IERNEW)
C
C***BEGIN PROLOGUE DNSID
C***REFER TO DDASPK
C***DATE WRITTEN 940701 (YYMMDD)
C***REVISION DATE 950713 (YYMMDD)
C***REVISION DATE 000628 TSCALE argument added.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DNSID solves a nonlinear system of algebraic equations of the
C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME
C in the initial conditions.
C
C The method used is a modified Newton scheme.
C
C The parameters represent
C
C X -- Independent variable.
C Y -- Solution vector.
C YPRIME -- Derivative of solution vector.
C NEQ -- Number of unknowns.
C ICOPT -- Initial condition option chosen (1 or 2).
C ID -- Array of dimension NEQ, which must be initialized
C if ICOPT = 1. See DDASIC.
C RES -- External user-supplied subroutine to evaluate the
C residual. See RES description in DDASPK prologue.
C WT -- Vector of weights for error criterion.
C RPAR,IPAR -- Real and integer arrays used for communication
C between the calling program and external user
C routines. They are not altered within DASPK.
C DELTA -- Residual vector on entry, and work vector of
C length NEQ for DNSID.
C WM,IWM -- Real and integer arrays storing matrix information
C such as the matrix of partial derivatives,
C permutation vector, and various other information.
C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
C R -- Array of length NEQ used as workspace by the
C linesearch routine DLINSD.
C YIC,YPIC -- Work vectors for DLINSD, each of length NEQ.
C EPCON -- Tolerance to test for convergence of the Newton
C iteration.
C RATEMX -- Maximum convergence rate for which Newton iteration
C is considered converging.
C MAXIT -- Maximum allowed number of Newton iterations.
C STPTOL -- Tolerance used in calculating the minimum lambda
C value allowed.
C ICNFLG -- Integer scalar. If nonzero, then constraint
C violations in the proposed new approximate solution
C will be checked for, and the maximum step length
C will be adjusted accordingly.
C ICNSTR -- Integer array of length NEQ containing flags for
C checking constraints.
C IERNEW -- Error flag for Newton iteration.
C 0 ==> Newton iteration converged.
C 1 ==> failed to converge, but RATE .le. RATEMX.
C 2 ==> failed to converge, RATE .gt. RATEMX.
C 3 ==> other recoverable error (IRES = -1, or
C linesearch failed).
C -1 ==> unrecoverable error (IRES = -2).
C
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C DSLVD, DDWNRM, DLINSD, DCOPY
C
C***END PROLOGUE DNSID
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),WT(*),R(*)
DIMENSION ID(*),DELTA(*), YIC(*), YPIC(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
DIMENSION ICNSTR(*)
EXTERNAL RES
C
PARAMETER (LNNI=19, LLSOFF=35)
C
C
C Initializations. M is the Newton iteration counter.
C
LSOFF = IWM(LLSOFF)
M = 0
RATE = 1.0D0
RLX = 0.4D0
C
C Compute a new step vector DELTA by back-substitution.
C
CALL DSLVD (NEQ, DELTA, WM, IWM)
C
C Get norm of DELTA. Return now if norm(DELTA) .le. EPCON.
C
DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
FNRM = DELNRM
IF (TSCALE .GT. 0.0D0) FNRM = FNRM*TSCALE*ABS(CJ)
IF (FNRM .LE. EPCON) RETURN
C
C Newton iteration loop.
C
300 CONTINUE
IWM(LNNI) = IWM(LNNI) + 1
C
C Call linesearch routine for global strategy and set RATE
C
OLDFNM = FNRM
C
CALL DLINSD (NEQ, Y, X, YPRIME, CJ, TSCALE, DELTA, DELNRM, WT,
* LSOFF, STPTOL, IRET, RES, IRES, WM, IWM, FNRM, ICOPT,
* ID, R, YIC, YPIC, ICNFLG, ICNSTR, RLX, RPAR, IPAR)
C
RATE = FNRM/OLDFNM
C
C Check for error condition from linesearch.
IF (IRET .NE. 0) GO TO 390
C
C Test for convergence of the iteration, and return or loop.
C
IF (FNRM .LE. EPCON) RETURN
C
C The iteration has not yet converged. Update M.
C Test whether the maximum number of iterations have been tried.
C
M = M + 1
IF (M .GE. MAXIT) GO TO 380
C
C Copy the residual to DELTA and its norm to DELNRM, and loop for
C another iteration.
C
CALL DCOPY (NEQ, R, 1, DELTA, 1)
DELNRM = FNRM
GO TO 300
C
C The maximum number of iterations was done. Set IERNEW and return.
C
380 IF (RATE .LE. RATEMX) THEN
IERNEW = 1
ELSE
IERNEW = 2
ENDIF
RETURN
C
390 IF (IRES .LE. -2) THEN
IERNEW = -1
ELSE
IERNEW = 3
ENDIF
RETURN
C
C
C------END OF SUBROUTINE DNSID------------------------------------------
END
SUBROUTINE DLINSD (NEQ, Y, T, YPRIME, CJ, TSCALE, P, PNRM, WT,
* LSOFF, STPTOL, IRET, RES, IRES, WM, IWM,
* FNRM, ICOPT, ID, R, YNEW, YPNEW, ICNFLG,
* ICNSTR, RLX, RPAR, IPAR)
C
C***BEGIN PROLOGUE DLINSD
C***REFER TO DNSID
C***DATE WRITTEN 941025 (YYMMDD)
C***REVISION DATE 941215 (YYMMDD)
C***REVISION DATE 960129 Moved line RL = ONE to top block.
C***REVISION DATE 000628 TSCALE argument added.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DLINSD uses a linesearch algorithm to calculate a new (Y,YPRIME)
C pair (YNEW,YPNEW) such that
C
C f(YNEW,YPNEW) .le. (1 - 2*ALPHA*RL)*f(Y,YPRIME) ,
C
C where 0 < RL <= 1. Here, f(y,y') is defined as
C
C f(y,y') = (1/2)*norm( (J-inverse)*G(t,y,y') )**2 ,
C
C where norm() is the weighted RMS vector norm, G is the DAE
C system residual function, and J is the system iteration matrix
C (Jacobian).
C
C In addition to the parameters defined elsewhere, we have
C
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
C P -- Approximate Newton step used in backtracking.
C PNRM -- Weighted RMS norm of P.
C LSOFF -- Flag showing whether the linesearch algorithm is
C to be invoked. 0 means do the linesearch, and
C 1 means turn off linesearch.
C STPTOL -- Tolerance used in calculating the minimum lambda
C value allowed.
C ICNFLG -- Integer scalar. If nonzero, then constraint violations
C in the proposed new approximate solution will be
C checked for, and the maximum step length will be
C adjusted accordingly.
C ICNSTR -- Integer array of length NEQ containing flags for
C checking constraints.
C RLX -- Real scalar restricting update size in DCNSTR.
C YNEW -- Array of length NEQ used to hold the new Y in
C performing the linesearch.
C YPNEW -- Array of length NEQ used to hold the new YPRIME in
C performing the linesearch.
C Y -- Array of length NEQ containing the new Y (i.e.,=YNEW).
C YPRIME -- Array of length NEQ containing the new YPRIME
C (i.e.,=YPNEW).
C FNRM -- Real scalar containing SQRT(2*f(Y,YPRIME)) for the
C current (Y,YPRIME) on input and output.
C R -- Work array of length NEQ, containing the scaled
C residual (J-inverse)*G(t,y,y') on return.
C IRET -- Return flag.
C IRET=0 means that a satisfactory (Y,YPRIME) was found.
C IRET=1 means that the routine failed to find a new
C (Y,YPRIME) that was sufficiently distinct from
C the current (Y,YPRIME) pair.
C IRET=2 means IRES .ne. 0 from RES.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C DFNRMD, DYYPNW, DCNSTR, DCOPY, XERRWD
C
C***END PROLOGUE DLINSD
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
EXTERNAL RES
DIMENSION Y(*), YPRIME(*), WT(*), R(*), ID(*)
DIMENSION WM(*), IWM(*)
DIMENSION YNEW(*), YPNEW(*), P(*), ICNSTR(*)
DIMENSION RPAR(*), IPAR(*)
CHARACTER MSG*80
C
PARAMETER (LNRE=12, LKPRIN=31)
C
SAVE ALPHA, ONE, TWO
DATA ALPHA/1.0D-4/, ONE/1.0D0/, TWO/2.0D0/
C
KPRIN=IWM(LKPRIN)
C
F1NRM = (FNRM*FNRM)/TWO
RATIO = ONE
IF (KPRIN .GE. 2) THEN
MSG = '------ IN ROUTINE DLINSD-- PNRM = (R1)'
CALL XERRWD(MSG, 38, 901, 0, 0, 0, 0, 1, PNRM, 0.0D0)
ENDIF
TAU = PNRM
RL = ONE
C-----------------------------------------------------------------------
C Check for violations of the constraints, if any are imposed.
C If any violations are found, the step vector P is rescaled, and the
C constraint check is repeated, until no violations are found.
C-----------------------------------------------------------------------
IF (ICNFLG .NE. 0) THEN
10 CONTINUE
CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW)
CALL DCNSTR (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
IF (IRET .EQ. 1) THEN
RATIO1 = TAU/PNRM
RATIO = RATIO*RATIO1
DO 20 I = 1,NEQ
20 P(I) = P(I)*RATIO1
PNRM = TAU
IF (KPRIN .GE. 2) THEN
MSG = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)'
CALL XERRWD(MSG, 50, 902, 0, 1, IVAR, 0, 1, PNRM, 0.0D0)
ENDIF
IF (PNRM .LE. STPTOL) THEN
IRET = 1
RETURN
ENDIF
GO TO 10
ENDIF
ENDIF
C
SLPI = (-TWO*F1NRM)*RATIO
RLMIN = STPTOL/PNRM
IF (LSOFF .EQ. 0 .AND. KPRIN .GE. 2) THEN
MSG = '------ MIN. LAMBDA = (R1)'
CALL XERRWD(MSG, 25, 903, 0, 0, 0, 0, 1, RLMIN, 0.0D0)
ENDIF
C-----------------------------------------------------------------------
C Begin iteration to find RL value satisfying alpha-condition.
C If RL becomes less than RLMIN, then terminate with IRET = 1.
C-----------------------------------------------------------------------
100 CONTINUE
CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW)
CALL DFNRMD (NEQ, YNEW, T, YPNEW, R, CJ, TSCALE, WT, RES, IRES,
* FNRMP, WM, IWM, RPAR, IPAR)
IWM(LNRE) = IWM(LNRE) + 1
IF (IRES .NE. 0) THEN
IRET = 2
RETURN
ENDIF
IF (LSOFF .EQ. 1) GO TO 150
C
F1NRMP = FNRMP*FNRMP/TWO
IF (KPRIN .GE. 2) THEN
MSG = '------ LAMBDA = (R1)'
CALL XERRWD(MSG, 20, 904, 0, 0, 0, 0, 1, RL, 0.0D0)
MSG = '------ NORM(F1) = (R1), NORM(F1NEW) = (R2)'
CALL XERRWD(MSG, 43, 905, 0, 0, 0, 0, 2, F1NRM, F1NRMP)
ENDIF
IF (F1NRMP .GT. F1NRM + ALPHA*SLPI*RL) GO TO 200
C-----------------------------------------------------------------------
C Alpha-condition is satisfied, or linesearch is turned off.
C Copy YNEW,YPNEW to Y,YPRIME and return.
C-----------------------------------------------------------------------
150 IRET = 0
CALL DCOPY (NEQ, YNEW, 1, Y, 1)
CALL DCOPY (NEQ, YPNEW, 1, YPRIME, 1)
FNRM = FNRMP
IF (KPRIN .GE. 1) THEN
MSG = '------ LEAVING ROUTINE DLINSD, FNRM = (R1)'
CALL XERRWD(MSG, 42, 906, 0, 0, 0, 0, 1, FNRM, 0.0D0)
ENDIF
RETURN
C-----------------------------------------------------------------------
C Alpha-condition not satisfied. Perform backtrack to compute new RL
C value. If no satisfactory YNEW,YPNEW can be found sufficiently
C distinct from Y,YPRIME, then return IRET = 1.
C-----------------------------------------------------------------------
200 CONTINUE
IF (RL .LT. RLMIN) THEN
IRET = 1
RETURN
ENDIF
C
RL = RL/TWO
GO TO 100
C
C----------------------- END OF SUBROUTINE DLINSD ----------------------
END
SUBROUTINE DFNRMD (NEQ, Y, T, YPRIME, R, CJ, TSCALE, WT,
* RES, IRES, FNORM, WM, IWM, RPAR, IPAR)
C
C***BEGIN PROLOGUE DFNRMD
C***REFER TO DLINSD
C***DATE WRITTEN 941025 (YYMMDD)
C***REVISION DATE 000628 TSCALE argument added.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DFNRMD calculates the scaled preconditioned norm of the nonlinear
C function used in the nonlinear iteration for obtaining consistent
C initial conditions. Specifically, DFNRMD calculates the weighted
C root-mean-square norm of the vector (J-inverse)*G(T,Y,YPRIME),
C where J is the Jacobian matrix.
C
C In addition to the parameters described in the calling program
C DLINSD, the parameters represent
C
C R -- Array of length NEQ that contains
C (J-inverse)*G(T,Y,YPRIME) on return.
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
C FNORM -- Scalar containing the weighted norm of R on return.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C RES, DSLVD, DDWNRM
C
C***END PROLOGUE DFNRMD
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
EXTERNAL RES
DIMENSION Y(*), YPRIME(*), WT(*), R(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
C-----------------------------------------------------------------------
C Call RES routine.
C-----------------------------------------------------------------------
IRES = 0
CALL RES(T,Y,YPRIME,CJ,R,IRES,RPAR,IPAR)
IF (IRES .LT. 0) RETURN
C-----------------------------------------------------------------------
C Apply inverse of Jacobian to vector R.
C-----------------------------------------------------------------------
CALL DSLVD(NEQ,R,WM,IWM)
C-----------------------------------------------------------------------
C Calculate norm of R.
C-----------------------------------------------------------------------
FNORM = DDWNRM(NEQ,R,WT,RPAR,IPAR)
IF (TSCALE .GT. 0.0D0) FNORM = FNORM*TSCALE*ABS(CJ)
C
RETURN
C----------------------- END OF SUBROUTINE DFNRMD ----------------------
END
SUBROUTINE DNEDD(X,Y,YPRIME,NEQ,RES,JACD,PDUM,H,WT,
* JSTART,IDID,RPAR,IPAR,PHI,GAMMA,DUMSVR,DELTA,E,
* WM,IWM,CJ,CJOLD,CJLAST,S,UROUND,DUME,DUMS,DUMR,
* EPCON,JCALC,JFDUM,KP1,NONNEG,NTYPE,IERNLS)
C
C***BEGIN PROLOGUE DNEDD
C***REFER TO DDASPK
C***DATE WRITTEN 891219 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DNEDD solves a nonlinear system of
C algebraic equations of the form
C G(X,Y,YPRIME) = 0 for the unknown Y.
C
C The method used is a modified Newton scheme.
C
C The parameters represent
C
C X -- Independent variable.
C Y -- Solution vector.
C YPRIME -- Derivative of solution vector.
C NEQ -- Number of unknowns.
C RES -- External user-supplied subroutine
C to evaluate the residual. See RES description
C in DDASPK prologue.
C JACD -- External user-supplied routine to evaluate the
C Jacobian. See JAC description for the case
C INFO(12) = 0 in the DDASPK prologue.
C PDUM -- Dummy argument.
C H -- Appropriate step size for next step.
C WT -- Vector of weights for error criterion.
C JSTART -- Indicates first call to this routine.
C If JSTART = 0, then this is the first call,
C otherwise it is not.
C IDID -- Completion flag, output by DNEDD.
C See IDID description in DDASPK prologue.
C RPAR,IPAR -- Real and integer arrays used for communication
C between the calling program and external user
C routines. They are not altered within DASPK.
C PHI -- Array of divided differences used by
C DNEDD. The length is NEQ*(K+1),where
C K is the maximum order.
C GAMMA -- Array used to predict Y and YPRIME. The length
C is MAXORD+1 where MAXORD is the maximum order.
C DUMSVR -- Dummy argument.
C DELTA -- Work vector for NLS of length NEQ.
C E -- Error accumulation vector for NLS of length NEQ.
C WM,IWM -- Real and integer arrays storing
C matrix information such as the matrix
C of partial derivatives, permutation
C vector, and various other information.
C CJ -- Parameter always proportional to 1/H.
C CJOLD -- Saves the value of CJ as of the last call to DMATD.
C Accounts for changes in CJ needed to
C decide whether to call DMATD.
C CJLAST -- Previous value of CJ.
C S -- A scalar determined by the approximate rate
C of convergence of the Newton iteration and used
C in the convergence test for the Newton iteration.
C
C If RATE is defined to be an estimate of the
C rate of convergence of the Newton iteration,
C then S = RATE/(1.D0-RATE).
C
C The closer RATE is to 0., the faster the Newton
C iteration is converging; the closer RATE is to 1.,
C the slower the Newton iteration is converging.
C
C On the first Newton iteration with an up-dated
C preconditioner S = 100.D0, Thus the initial
C RATE of convergence is approximately 1.
C
C S is preserved from call to call so that the rate
C estimate from a previous step can be applied to
C the current step.
C UROUND -- Unit roundoff.
C DUME -- Dummy argument.
C DUMS -- Dummy argument.
C DUMR -- Dummy argument.
C EPCON -- Tolerance to test for convergence of the Newton
C iteration.
C JCALC -- Flag used to determine when to update
C the Jacobian matrix. In general:
C
C JCALC = -1 ==> Call the DMATD routine to update
C the Jacobian matrix.
C JCALC = 0 ==> Jacobian matrix is up-to-date.
C JCALC = 1 ==> Jacobian matrix is out-dated,
C but DMATD will not be called unless
C JCALC is set to -1.
C JFDUM -- Dummy argument.
C KP1 -- The current order(K) + 1; updated across calls.
C NONNEG -- Flag to determine nonnegativity constraints.
C NTYPE -- Identification code for the NLS routine.
C 0 ==> modified Newton; direct solver.
C IERNLS -- Error flag for nonlinear solver.
C 0 ==> nonlinear solver converged.
C 1 ==> recoverable error inside nonlinear solver.
C -1 ==> unrecoverable error inside nonlinear solver.
C
C All variables with "DUM" in their names are dummy variables
C which are not used in this routine.
C
C Following is a list and description of local variables which
C may not have an obvious usage. They are listed in roughly the
C order they occur in this subroutine.
C
C The following group of variables are passed as arguments to
C the Newton iteration solver. They are explained in greater detail
C in DNSD:
C TOLNEW, MULDEL, MAXIT, IERNEW
C
C IERTYP -- Flag which tells whether this subroutine is correct.
C 0 ==> correct subroutine.
C 1 ==> incorrect subroutine.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C DDWNRM, RES, DMATD, DNSD
C
C***END PROLOGUE DNEDD
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),WT(*)
DIMENSION DELTA(*),E(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
DIMENSION PHI(NEQ,*),GAMMA(*)
EXTERNAL RES, JACD
C
PARAMETER (LNRE=12, LNJE=13)
C
SAVE MULDEL, MAXIT, XRATE
DATA MULDEL/1/, MAXIT/4/, XRATE/0.25D0/
C
C Verify that this is the correct subroutine.
C
IERTYP = 0
IF (NTYPE .NE. 0) THEN
IERTYP = 1
GO TO 380
ENDIF
C
C If this is the first step, perform initializations.
C
IF (JSTART .EQ. 0) THEN
CJOLD = CJ
JCALC = -1
ENDIF
C
C Perform all other initializations.
C
IERNLS = 0
C
C Decide whether new Jacobian is needed.
C
TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE)
TEMP2 = 1.0D0/TEMP1
IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1
IF (CJ .NE. CJLAST) S = 100.D0
C
C-----------------------------------------------------------------------
C Entry point for updating the Jacobian with current
C stepsize.
C-----------------------------------------------------------------------
300 CONTINUE
C
C Initialize all error flags to zero.
C
IERJ = 0
IRES = 0
IERNEW = 0
C
C Predict the solution and derivative and compute the tolerance
C for the Newton iteration.
C
DO 310 I=1,NEQ
Y(I)=PHI(I,1)
310 YPRIME(I)=0.0D0
DO 330 J=2,KP1
DO 320 I=1,NEQ
Y(I)=Y(I)+PHI(I,J)
320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J)
330 CONTINUE
PNORM = DDWNRM (NEQ,Y,WT,RPAR,IPAR)
TOLNEW = 100.D0*UROUND*PNORM
C
C Call RES to initialize DELTA.
C
IWM(LNRE)=IWM(LNRE)+1
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
IF (IRES .LT. 0) GO TO 380
C
C If indicated, reevaluate the iteration matrix
C J = dG/dY + CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0).
C Set JCALC to 0 as an indicator that this has been done.
C
IF(JCALC .EQ. -1) THEN
IWM(LNJE)=IWM(LNJE)+1
JCALC=0
CALL DMATD(NEQ,X,Y,YPRIME,DELTA,CJ,H,IERJ,WT,E,WM,IWM,
* RES,IRES,UROUND,JACD,RPAR,IPAR)
CJOLD=CJ
S = 100.D0
IF (IRES .LT. 0) GO TO 380
IF(IERJ .NE. 0)GO TO 380
ENDIF
C
C Call the nonlinear Newton solver.
C
TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD)
CALL DNSD(X,Y,YPRIME,NEQ,RES,PDUM,WT,RPAR,IPAR,DUMSVR,
* DELTA,E,WM,IWM,CJ,DUMS,DUMR,DUME,EPCON,S,TEMP1,
* TOLNEW,MULDEL,MAXIT,IRES,IDUM,IERNEW)
C
IF (IERNEW .GT. 0 .AND. JCALC .NE. 0) THEN
C
C The Newton iteration had a recoverable failure with an old
C iteration matrix. Retry the step with a new iteration matrix.
C
JCALC = -1
GO TO 300
ENDIF
C
IF (IERNEW .NE. 0) GO TO 380
C
C The Newton iteration has converged. If nonnegativity of
C solution is required, set the solution nonnegative, if the
C perturbation to do it is small enough. If the change is too
C large, then consider the corrector iteration to have failed.
C
375 IF(NONNEG .EQ. 0) GO TO 390
DO 377 I = 1,NEQ
377 DELTA(I) = MIN(Y(I),0.0D0)
DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
IF(DELNRM .GT. EPCON) GO TO 380
DO 378 I = 1,NEQ
378 E(I) = E(I) - DELTA(I)
GO TO 390
C
C
C Exits from nonlinear solver.
C No convergence with current iteration
C matrix, or singular iteration matrix.
C Compute IERNLS and IDID accordingly.
C
380 CONTINUE
IF (IRES .LE. -2 .OR. IERTYP .NE. 0) THEN
IERNLS = -1
IF (IRES .LE. -2) IDID = -11
IF (IERTYP .NE. 0) IDID = -15
ELSE
IERNLS = 1
IF (IRES .LT. 0) IDID = -10
IF (IERJ .NE. 0) IDID = -8
ENDIF
C
390 JCALC = 1
RETURN
C
C------END OF SUBROUTINE DNEDD------------------------------------------
END
SUBROUTINE DNSD(X,Y,YPRIME,NEQ,RES,PDUM,WT,RPAR,IPAR,
* DUMSVR,DELTA,E,WM,IWM,CJ,DUMS,DUMR,DUME,EPCON,
* S,CONFAC,TOLNEW,MULDEL,MAXIT,IRES,IDUM,IERNEW)
C
C***BEGIN PROLOGUE DNSD
C***REFER TO DDASPK
C***DATE WRITTEN 891219 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***REVISION DATE 950126 (YYMMDD)
C***REVISION DATE 000711 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DNSD solves a nonlinear system of
C algebraic equations of the form
C G(X,Y,YPRIME) = 0 for the unknown Y.
C
C The method used is a modified Newton scheme.
C
C The parameters represent
C
C X -- Independent variable.
C Y -- Solution vector.
C YPRIME -- Derivative of solution vector.
C NEQ -- Number of unknowns.
C RES -- External user-supplied subroutine
C to evaluate the residual. See RES description
C in DDASPK prologue.
C PDUM -- Dummy argument.
C WT -- Vector of weights for error criterion.
C RPAR,IPAR -- Real and integer arrays used for communication
C between the calling program and external user
C routines. They are not altered within DASPK.
C DUMSVR -- Dummy argument.
C DELTA -- Work vector for DNSD of length NEQ.
C E -- Error accumulation vector for DNSD of length NEQ.
C WM,IWM -- Real and integer arrays storing
C matrix information such as the matrix
C of partial derivatives, permutation
C vector, and various other information.
C CJ -- Parameter always proportional to 1/H (step size).
C DUMS -- Dummy argument.
C DUMR -- Dummy argument.
C DUME -- Dummy argument.
C EPCON -- Tolerance to test for convergence of the Newton
C iteration.
C S -- Used for error convergence tests.
C In the Newton iteration: S = RATE/(1 - RATE),
C where RATE is the estimated rate of convergence
C of the Newton iteration.
C The calling routine passes the initial value
C of S to the Newton iteration.
C CONFAC -- A residual scale factor to improve convergence.
C TOLNEW -- Tolerance on the norm of Newton correction in
C alternative Newton convergence test.
C MULDEL -- A flag indicating whether or not to multiply
C DELTA by CONFAC.
C 0 ==> do not scale DELTA by CONFAC.
C 1 ==> scale DELTA by CONFAC.
C MAXIT -- Maximum allowed number of Newton iterations.
C IRES -- Error flag returned from RES. See RES description
C in DDASPK prologue. If IRES = -1, then IERNEW
C will be set to 1.
C If IRES < -1, then IERNEW will be set to -1.
C IDUM -- Dummy argument.
C IERNEW -- Error flag for Newton iteration.
C 0 ==> Newton iteration converged.
C 1 ==> recoverable error inside Newton iteration.
C -1 ==> unrecoverable error inside Newton iteration.
C
C All arguments with "DUM" in their names are dummy arguments
C which are not used in this routine.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C DSLVD, DDWNRM, RES
C
C***END PROLOGUE DNSD
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),WT(*),DELTA(*),E(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
EXTERNAL RES
C
PARAMETER (LNRE=12, LNNI=19)
C
C Initialize Newton counter M and accumulation vector E.
C
M = 0
DO 100 I=1,NEQ
100 E(I)=0.0D0
C
C Corrector loop.
C
300 CONTINUE
IWM(LNNI) = IWM(LNNI) + 1
C
C If necessary, multiply residual by convergence factor.
C
IF (MULDEL .EQ. 1) THEN
DO 320 I = 1,NEQ
320 DELTA(I) = DELTA(I) * CONFAC
ENDIF
C
C Compute a new iterate (back-substitution).
C Store the correction in DELTA.
C
CALL DSLVD(NEQ,DELTA,WM,IWM)
C
C Update Y, E, and YPRIME.
C
DO 340 I=1,NEQ
Y(I)=Y(I)-DELTA(I)
E(I)=E(I)-DELTA(I)
340 YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
C
C Test for convergence of the iteration.
C
DELNRM=DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
IF (M .EQ. 0) THEN
OLDNRM = DELNRM
IF (DELNRM .LE. TOLNEW) GO TO 370
ELSE
RATE = (DELNRM/OLDNRM)**(1.0D0/M)
IF (RATE .GT. 0.9D0) GO TO 380
S = RATE/(1.0D0 - RATE)
ENDIF
IF (S*DELNRM .LE. EPCON) GO TO 370
C
C The corrector has not yet converged.
C Update M and test whether the
C maximum number of iterations have
C been tried.
C
M=M+1
IF(M.GE.MAXIT) GO TO 380
C
C Evaluate the residual,
C and go back to do another iteration.
C
IWM(LNRE)=IWM(LNRE)+1
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
IF (IRES .LT. 0) GO TO 380
GO TO 300
C
C The iteration has converged.
C
370 RETURN
C
C The iteration has not converged. Set IERNEW appropriately.
C
380 CONTINUE
IF (IRES .LE. -2 ) THEN
IERNEW = -1
ELSE
IERNEW = 1
ENDIF
RETURN
C
C
C------END OF SUBROUTINE DNSD-------------------------------------------
END
SUBROUTINE DMATD(NEQ,X,Y,YPRIME,DELTA,CJ,H,IER,EWT,E,
* WM,IWM,RES,IRES,UROUND,JACD,RPAR,IPAR)
C
C***BEGIN PROLOGUE DMATD
C***REFER TO DDASPK
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***REVISION DATE 940701 (YYMMDD) (new LIPVT)
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This routine computes the iteration matrix
C J = dG/dY+CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0).
C Here J is computed by:
C the user-supplied routine JACD if IWM(MTYPE) is 1 or 4, or
C by numerical difference quotients if IWM(MTYPE) is 2 or 5.
C
C The parameters have the following meanings.
C X = Independent variable.
C Y = Array containing predicted values.
C YPRIME = Array containing predicted derivatives.
C DELTA = Residual evaluated at (X,Y,YPRIME).
C (Used only if IWM(MTYPE)=2 or 5).
C CJ = Scalar parameter defining iteration matrix.
C H = Current stepsize in integration.
C IER = Variable which is .NE. 0 if iteration matrix
C is singular, and 0 otherwise.
C EWT = Vector of error weights for computing norms.
C E = Work space (temporary) of length NEQ.
C WM = Real work space for matrices. On output
C it contains the LU decomposition
C of the iteration matrix.
C IWM = Integer work space containing
C matrix information.
C RES = External user-supplied subroutine
C to evaluate the residual. See RES description
C in DDASPK prologue.
C IRES = Flag which is equal to zero if no illegal values
C in RES, and less than zero otherwise. (If IRES
C is less than zero, the matrix was not completed).
C In this case (if IRES .LT. 0), then IER = 0.
C UROUND = The unit roundoff error of the machine being used.
C JACD = Name of the external user-supplied routine
C to evaluate the iteration matrix. (This routine
C is only used if IWM(MTYPE) is 1 or 4)
C See JAC description for the case INFO(12) = 0
C in DDASPK prologue.
C RPAR,IPAR= Real and integer parameter arrays that
C are used for communication between the
C calling program and external user routines.
C They are not altered by DMATD.
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C JACD, RES, DGEFA, DGBFA
C
C***END PROLOGUE DMATD
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),DELTA(*),EWT(*),E(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
EXTERNAL RES, JACD
C
PARAMETER (LML=1, LMU=2, LMTYPE=4, LNRE=12, LNPD=22, LLCIWP=30)
C
LIPVT = IWM(LLCIWP)
IER = 0
MTYPE=IWM(LMTYPE)
GO TO (100,200,300,400,500),MTYPE
C
C
C Dense user-supplied matrix.
C
100 LENPD=IWM(LNPD)
DO 110 I=1,LENPD
110 WM(I)=0.0D0
CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR)
GO TO 230
C
C
C Dense finite-difference-generated matrix.
C
200 IRES=0
NROW=0
SQUR = SQRT(UROUND)
DO 210 I=1,NEQ
DEL=SQUR*MAX(ABS(Y(I)),ABS(H*YPRIME(I)),
* ABS(1.D0/EWT(I)))
DEL=SIGN(DEL,H*YPRIME(I))
DEL=(Y(I)+DEL)-Y(I)
YSAVE=Y(I)
YPSAVE=YPRIME(I)
Y(I)=Y(I)+DEL
YPRIME(I)=YPRIME(I)+CJ*DEL
IWM(LNRE)=IWM(LNRE)+1
CALL RES(X,Y,YPRIME,CJ,E,IRES,RPAR,IPAR)
IF (IRES .LT. 0) RETURN
DELINV=1.0D0/DEL
DO 220 L=1,NEQ
220 WM(NROW+L)=(E(L)-DELTA(L))*DELINV
NROW=NROW+NEQ
Y(I)=YSAVE
YPRIME(I)=YPSAVE
210 CONTINUE
C
C
C Do dense-matrix LU decomposition on J.
C
230 CALL DGEFA(WM,NEQ,NEQ,IWM(LIPVT),IER)
RETURN
C
C
C Dummy section for IWM(MTYPE)=3.
C
300 RETURN
C
C
C Banded user-supplied matrix.
C
400 LENPD=IWM(LNPD)
DO 410 I=1,LENPD
410 WM(I)=0.0D0
CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR)
MEBAND=2*IWM(LML)+IWM(LMU)+1
GO TO 550
C
C
C Banded finite-difference-generated matrix.
C
500 MBAND=IWM(LML)+IWM(LMU)+1
MBA=MIN0(MBAND,NEQ)
MEBAND=MBAND+IWM(LML)
MEB1=MEBAND-1
MSAVE=(NEQ/MBAND)+1
ISAVE=IWM(LNPD)
IPSAVE=ISAVE+MSAVE
IRES=0
SQUR=SQRT(UROUND)
DO 540 J=1,MBA
DO 510 N=J,NEQ,MBAND
K= (N-J)/MBAND + 1
WM(ISAVE+K)=Y(N)
WM(IPSAVE+K)=YPRIME(N)
DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),
* ABS(1.D0/EWT(N)))
DEL=SIGN(DEL,H*YPRIME(N))
DEL=(Y(N)+DEL)-Y(N)
Y(N)=Y(N)+DEL
510 YPRIME(N)=YPRIME(N)+CJ*DEL
IWM(LNRE)=IWM(LNRE)+1
CALL RES(X,Y,YPRIME,CJ,E,IRES,RPAR,IPAR)
IF (IRES .LT. 0) RETURN
DO 530 N=J,NEQ,MBAND
K= (N-J)/MBAND + 1
Y(N)=WM(ISAVE+K)
YPRIME(N)=WM(IPSAVE+K)
DEL=SQUR*MAX(ABS(Y(N)),ABS(H*YPRIME(N)),
* ABS(1.D0/EWT(N)))
DEL=SIGN(DEL,H*YPRIME(N))
DEL=(Y(N)+DEL)-Y(N)
DELINV=1.0D0/DEL
I1=MAX0(1,(N-IWM(LMU)))
I2=MIN0(NEQ,(N+IWM(LML)))
II=N*MEB1-IWM(LML)
DO 520 I=I1,I2
520 WM(II+I)=(E(I)-DELTA(I))*DELINV
530 CONTINUE
540 CONTINUE
C
C
C Do LU decomposition of banded J.
C
550 CALL DGBFA (WM,MEBAND,NEQ,IWM(LML),IWM(LMU),IWM(LIPVT),IER)
RETURN
C
C------END OF SUBROUTINE DMATD------------------------------------------
END
SUBROUTINE DSLVD(NEQ,DELTA,WM,IWM)
C
C***BEGIN PROLOGUE DSLVD
C***REFER TO DDASPK
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***REVISION DATE 940701 (YYMMDD) (new LIPVT)
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This routine manages the solution of the linear
C system arising in the Newton iteration.
C Real matrix information and real temporary storage
C is stored in the array WM.
C Integer matrix information is stored in the array IWM.
C For a dense matrix, the LINPACK routine DGESL is called.
C For a banded matrix, the LINPACK routine DGBSL is called.
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C DGESL, DGBSL
C
C***END PROLOGUE DSLVD
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION DELTA(*),WM(*),IWM(*)
C
PARAMETER (LML=1, LMU=2, LMTYPE=4, LLCIWP=30)
C
LIPVT = IWM(LLCIWP)
MTYPE=IWM(LMTYPE)
GO TO(100,100,300,400,400),MTYPE
C
C Dense matrix.
C
100 CALL DGESL(WM,NEQ,NEQ,IWM(LIPVT),DELTA,0)
RETURN
C
C Dummy section for MTYPE=3.
C
300 CONTINUE
RETURN
C
C Banded matrix.
C
400 MEBAND=2*IWM(LML)+IWM(LMU)+1
CALL DGBSL(WM,MEBAND,NEQ,IWM(LML),
* IWM(LMU),IWM(LIPVT),DELTA,0)
RETURN
C
C------END OF SUBROUTINE DSLVD------------------------------------------
END
SUBROUTINE DDASIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JACK,PSOL,H,TSCALE,
* WT,JSKIP,RPAR,IPAR,SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,UROUND,
* EPLI,SQRTN,RSQRTN,EPCON,RATEMX,STPTOL,JFLG,
* ICNFLG,ICNSTR,IERNLS)
C
C***BEGIN PROLOGUE DDASIK
C***REFER TO DDASPK
C***DATE WRITTEN 941026 (YYMMDD)
C***REVISION DATE 950808 (YYMMDD)
C***REVISION DATE 951110 Removed unreachable block 390.
C***REVISION DATE 000628 TSCALE argument added.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C
C DDASIK solves a nonlinear system of algebraic equations of the
C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in
C the initial conditions.
C
C An initial value for Y and initial guess for YPRIME are input.
C
C The method used is a Newton scheme with Krylov iteration and a
C linesearch algorithm.
C
C The parameters represent
C
C X -- Independent variable.
C Y -- Solution vector at x.
C YPRIME -- Derivative of solution vector.
C NEQ -- Number of equations to be integrated.
C ICOPT -- Initial condition option chosen (1 or 2).
C ID -- Array of dimension NEQ, which must be initialized
C if ICOPT = 1. See DDASIC.
C RES -- External user-supplied subroutine
C to evaluate the residual. See RES description
C in DDASPK prologue.
C JACK -- External user-supplied routine to update
C the preconditioner. (This is optional).
C See JAC description for the case
C INFO(12) = 1 in the DDASPK prologue.
C PSOL -- External user-supplied routine to solve
C a linear system using preconditioning.
C (This is optional). See explanation inside DDASPK.
C H -- Scaling factor for this initial condition calc.
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
C WT -- Vector of weights for error criterion.
C JSKIP -- input flag to signal if initial JAC call is to be
C skipped. 1 => skip the call, 0 => do not skip call.
C RPAR,IPAR -- Real and integer arrays used for communication
C between the calling program and external user
C routines. They are not altered within DASPK.
C SAVR -- Work vector for DDASIK of length NEQ.
C DELTA -- Work vector for DDASIK of length NEQ.
C R -- Work vector for DDASIK of length NEQ.
C YIC,YPIC -- Work vectors for DDASIK, each of length NEQ.
C PWK -- Work vector for DDASIK of length NEQ.
C WM,IWM -- Real and integer arrays storing
C matrix information for linear system
C solvers, and various other information.
C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
C UROUND -- Unit roundoff. Not used here.
C EPLI -- convergence test constant.
C See DDASPK prologue for more details.
C SQRTN -- Square root of NEQ.
C RSQRTN -- reciprical of square root of NEQ.
C EPCON -- Tolerance to test for convergence of the Newton
C iteration.
C RATEMX -- Maximum convergence rate for which Newton iteration
C is considered converging.
C JFLG -- Flag showing whether a Jacobian routine is supplied.
C ICNFLG -- Integer scalar. If nonzero, then constraint
C violations in the proposed new approximate solution
C will be checked for, and the maximum step length
C will be adjusted accordingly.
C ICNSTR -- Integer array of length NEQ containing flags for
C checking constraints.
C IERNLS -- Error flag for nonlinear solver.
C 0 ==> nonlinear solver converged.
C 1,2 ==> recoverable error inside nonlinear solver.
C 1 => retry with current Y, YPRIME
C 2 => retry with original Y, YPRIME
C -1 ==> unrecoverable error in nonlinear solver.
C
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C RES, JACK, DNSIK, DCOPY
C
C***END PROLOGUE DDASIK
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),ID(*),WT(*),ICNSTR(*)
DIMENSION SAVR(*),DELTA(*),R(*),YIC(*),YPIC(*),PWK(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
EXTERNAL RES, JACK, PSOL
C
PARAMETER (LNRE=12, LNJE=13, LLOCWP=29, LLCIWP=30)
PARAMETER (LMXNIT=32, LMXNJ=33)
C
C
C Perform initializations.
C
LWP = IWM(LLOCWP)
LIWP = IWM(LLCIWP)
MXNIT = IWM(LMXNIT)
MXNJ = IWM(LMXNJ)
IERNLS = 0
NJ = 0
EPLIN = EPLI*EPCON
C
C Call RES to initialize DELTA.
C
IRES = 0
IWM(LNRE) = IWM(LNRE) + 1
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
IF (IRES .LT. 0) GO TO 370
C
C Looping point for updating the preconditioner.
C
300 CONTINUE
C
C Initialize all error flags to zero.
C
IERPJ = 0
IRES = 0
IERNEW = 0
C
C If a Jacobian routine was supplied, call it.
C
IF (JFLG .EQ. 1 .AND. JSKIP .EQ. 0) THEN
NJ = NJ + 1
IWM(LNJE)=IWM(LNJE)+1
CALL JACK (RES, IRES, NEQ, X, Y, YPRIME, WT, DELTA, R, H, CJ,
* WM(LWP), IWM(LIWP), IERPJ, RPAR, IPAR)
IF (IRES .LT. 0 .OR. IERPJ .NE. 0) GO TO 370
ENDIF
JSKIP = 0
C
C Call the nonlinear Newton solver for up to MXNIT iterations.
C
CALL DNSIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,PSOL,WT,RPAR,IPAR,
* SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,TSCALE,SQRTN,RSQRTN,
* EPLIN,EPCON,RATEMX,MXNIT,STPTOL,ICNFLG,ICNSTR,IERNEW)
C
IF (IERNEW .EQ. 1 .AND. NJ .LT. MXNJ .AND. JFLG .EQ. 1) THEN
C
C Up to MXNIT iterations were done, the convergence rate is < 1,
C a Jacobian routine is supplied, and the number of JACK calls
C is less than MXNJ.
C Copy the residual SAVR to DELTA, call JACK, and try again.
C
CALL DCOPY (NEQ, SAVR, 1, DELTA, 1)
GO TO 300
ENDIF
C
IF (IERNEW .NE. 0) GO TO 380
RETURN
C
C
C Unsuccessful exits from nonlinear solver.
C Set IERNLS accordingly.
C
370 IERNLS = 2
IF (IRES .LE. -2) IERNLS = -1
RETURN
C
380 IERNLS = MIN(IERNEW,2)
RETURN
C
C----------------------- END OF SUBROUTINE DDASIK-----------------------
END
SUBROUTINE DNSIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,PSOL,WT,RPAR,IPAR,
* SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,TSCALE,SQRTN,RSQRTN,EPLIN,
* EPCON,RATEMX,MAXIT,STPTOL,ICNFLG,ICNSTR,IERNEW)
C
C***BEGIN PROLOGUE DNSIK
C***REFER TO DDASPK
C***DATE WRITTEN 940701 (YYMMDD)
C***REVISION DATE 950714 (YYMMDD)
C***REVISION DATE 000628 TSCALE argument added.
C***REVISION DATE 000628 Added criterion for IERNEW = 1 return.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DNSIK solves a nonlinear system of algebraic equations of the
C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in
C the initial conditions.
C
C The method used is a Newton scheme combined with a linesearch
C algorithm, using Krylov iterative linear system methods.
C
C The parameters represent
C
C X -- Independent variable.
C Y -- Solution vector.
C YPRIME -- Derivative of solution vector.
C NEQ -- Number of unknowns.
C ICOPT -- Initial condition option chosen (1 or 2).
C ID -- Array of dimension NEQ, which must be initialized
C if ICOPT = 1. See DDASIC.
C RES -- External user-supplied subroutine
C to evaluate the residual. See RES description
C in DDASPK prologue.
C PSOL -- External user-supplied routine to solve
C a linear system using preconditioning.
C See explanation inside DDASPK.
C WT -- Vector of weights for error criterion.
C RPAR,IPAR -- Real and integer arrays used for communication
C between the calling program and external user
C routines. They are not altered within DASPK.
C SAVR -- Work vector for DNSIK of length NEQ.
C DELTA -- Residual vector on entry, and work vector of
C length NEQ for DNSIK.
C R -- Work vector for DNSIK of length NEQ.
C YIC,YPIC -- Work vectors for DNSIK, each of length NEQ.
C PWK -- Work vector for DNSIK of length NEQ.
C WM,IWM -- Real and integer arrays storing
C matrix information such as the matrix
C of partial derivatives, permutation
C vector, and various other information.
C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
C SQRTN -- Square root of NEQ.
C RSQRTN -- reciprical of square root of NEQ.
C EPLIN -- Tolerance for linear system solver.
C EPCON -- Tolerance to test for convergence of the Newton
C iteration.
C RATEMX -- Maximum convergence rate for which Newton iteration
C is considered converging.
C MAXIT -- Maximum allowed number of Newton iterations.
C STPTOL -- Tolerance used in calculating the minimum lambda
C value allowed.
C ICNFLG -- Integer scalar. If nonzero, then constraint
C violations in the proposed new approximate solution
C will be checked for, and the maximum step length
C will be adjusted accordingly.
C ICNSTR -- Integer array of length NEQ containing flags for
C checking constraints.
C IERNEW -- Error flag for Newton iteration.
C 0 ==> Newton iteration converged.
C 1 ==> failed to converge, but RATE .lt. 1, or the
C residual norm was reduced by a factor of .1.
C 2 ==> failed to converge, RATE .gt. RATEMX.
C 3 ==> other recoverable error.
C -1 ==> unrecoverable error inside Newton iteration.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C DFNRMK, DSLVK, DDWNRM, DLINSK, DCOPY
C
C***END PROLOGUE DNSIK
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),WT(*),ID(*),DELTA(*),R(*),SAVR(*)
DIMENSION YIC(*),YPIC(*),PWK(*),WM(*),IWM(*), RPAR(*),IPAR(*)
DIMENSION ICNSTR(*)
EXTERNAL RES, PSOL
C
PARAMETER (LNNI=19, LNPS=21, LLOCWP=29, LLCIWP=30)
PARAMETER (LLSOFF=35, LSTOL=14)
C
C
C Initializations. M is the Newton iteration counter.
C
LSOFF = IWM(LLSOFF)
M = 0
RATE = 1.0D0
LWP = IWM(LLOCWP)
LIWP = IWM(LLCIWP)
RLX = 0.4D0
C
C Save residual in SAVR.
C
CALL DCOPY (NEQ, DELTA, 1, SAVR, 1)
C
C Compute norm of (P-inverse)*(residual).
C
CALL DFNRMK (NEQ, Y, X, YPRIME, SAVR, R, CJ, TSCALE, WT,
* SQRTN, RSQRTN, RES, IRES, PSOL, 1, IER, FNRM, EPLIN,
* WM(LWP), IWM(LIWP), PWK, RPAR, IPAR)
IWM(LNPS) = IWM(LNPS) + 1
IF (IER .NE. 0) THEN
IERNEW = 3
RETURN
ENDIF
C
C Return now if residual norm is .le. EPCON.
C
IF (FNRM .LE. EPCON) RETURN
C
C Newton iteration loop.
C
FNRM0 = FNRM
300 CONTINUE
IWM(LNNI) = IWM(LNNI) + 1
C
C Compute a new step vector DELTA.
C
CALL DSLVK (NEQ, Y, X, YPRIME, SAVR, DELTA, WT, WM, IWM,
* RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK,
* RPAR, IPAR)
IF (IRES .NE. 0 .OR. IERSL .NE. 0) GO TO 390
C
C Get norm of DELTA. Return now if DELTA is zero.
C
DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
IF (DELNRM .EQ. 0.0D0) RETURN
C
C Call linesearch routine for global strategy and set RATE.
C
OLDFNM = FNRM
C
CALL DLINSK (NEQ, Y, X, YPRIME, SAVR, CJ, TSCALE, DELTA, DELNRM,
* WT, SQRTN, RSQRTN, LSOFF, STPTOL, IRET, RES, IRES, PSOL,
* WM, IWM, RHOK, FNRM, ICOPT, ID, WM(LWP), IWM(LIWP), R, EPLIN,
* YIC, YPIC, PWK, ICNFLG, ICNSTR, RLX, RPAR, IPAR)
C
RATE = FNRM/OLDFNM
C
C Check for error condition from linesearch.
IF (IRET .NE. 0) GO TO 390
C
C Test for convergence of the iteration, and return or loop.
C
IF (FNRM .LE. EPCON) RETURN
C
C The iteration has not yet converged. Update M.
C Test whether the maximum number of iterations have been tried.
C
M = M + 1
IF(M .GE. MAXIT) GO TO 380
C
C Copy the residual SAVR to DELTA and loop for another iteration.
C
CALL DCOPY (NEQ, SAVR, 1, DELTA, 1)
GO TO 300
C
C The maximum number of iterations was done. Set IERNEW and return.
C
380 IF (RATE .LE. RATEMX .OR. FNRM .LE. 0.1D0*FNRM0) THEN
IERNEW = 1
ELSE
IERNEW = 2
ENDIF
RETURN
C
390 IF (IRES .LE. -2 .OR. IERSL .LT. 0) THEN
IERNEW = -1
ELSE
IERNEW = 3
IF (IRES .EQ. 0 .AND. IERSL .EQ. 1 .AND. M .GE. 2
1 .AND. RATE .LT. 1.0D0) IERNEW = 1
ENDIF
RETURN
C
C
C----------------------- END OF SUBROUTINE DNSIK------------------------
END
SUBROUTINE DLINSK (NEQ, Y, T, YPRIME, SAVR, CJ, TSCALE, P, PNRM,
* WT, SQRTN, RSQRTN, LSOFF, STPTOL, IRET, RES, IRES, PSOL,
* WM, IWM, RHOK, FNRM, ICOPT, ID, WP, IWP, R, EPLIN, YNEW, YPNEW,
* PWK, ICNFLG, ICNSTR, RLX, RPAR, IPAR)
C
C***BEGIN PROLOGUE DLINSK
C***REFER TO DNSIK
C***DATE WRITTEN 940830 (YYMMDD)
C***REVISION DATE 951006 (Arguments SQRTN, RSQRTN added.)
C***REVISION DATE 960129 Moved line RL = ONE to top block.
C***REVISION DATE 000628 TSCALE argument added.
C***REVISION DATE 000628 RHOK*RHOK term removed in alpha test.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DLINSK uses a linesearch algorithm to calculate a new (Y,YPRIME)
C pair (YNEW,YPNEW) such that
C
C f(YNEW,YPNEW) .le. (1 - 2*ALPHA*RL)*f(Y,YPRIME)
C
C where 0 < RL <= 1, and RHOK is the scaled preconditioned norm of
C the final residual vector in the Krylov iteration.
C Here, f(y,y') is defined as
C
C f(y,y') = (1/2)*norm( (P-inverse)*G(t,y,y') )**2 ,
C
C where norm() is the weighted RMS vector norm, G is the DAE
C system residual function, and P is the preconditioner used
C in the Krylov iteration.
C
C In addition to the parameters defined elsewhere, we have
C
C SAVR -- Work array of length NEQ, containing the residual
C vector G(t,y,y') on return.
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
C P -- Approximate Newton step used in backtracking.
C PNRM -- Weighted RMS norm of P.
C LSOFF -- Flag showing whether the linesearch algorithm is
C to be invoked. 0 means do the linesearch,
C 1 means turn off linesearch.
C STPTOL -- Tolerance used in calculating the minimum lambda
C value allowed.
C ICNFLG -- Integer scalar. If nonzero, then constraint violations
C in the proposed new approximate solution will be
C checked for, and the maximum step length will be
C adjusted accordingly.
C ICNSTR -- Integer array of length NEQ containing flags for
C checking constraints.
C RHOK -- Weighted norm of preconditioned Krylov residual.
C RLX -- Real scalar restricting update size in DCNSTR.
C YNEW -- Array of length NEQ used to hold the new Y in
C performing the linesearch.
C YPNEW -- Array of length NEQ used to hold the new YPRIME in
C performing the linesearch.
C PWK -- Work vector of length NEQ for use in PSOL.
C Y -- Array of length NEQ containing the new Y (i.e.,=YNEW).
C YPRIME -- Array of length NEQ containing the new YPRIME
C (i.e.,=YPNEW).
C FNRM -- Real scalar containing SQRT(2*f(Y,YPRIME)) for the
C current (Y,YPRIME) on input and output.
C R -- Work space length NEQ for residual vector.
C IRET -- Return flag.
C IRET=0 means that a satisfactory (Y,YPRIME) was found.
C IRET=1 means that the routine failed to find a new
C (Y,YPRIME) that was sufficiently distinct from
C the current (Y,YPRIME) pair.
C IRET=2 means a failure in RES or PSOL.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C DFNRMK, DYYPNW, DCNSTR, DCOPY, XERRWD
C
C***END PROLOGUE DLINSK
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
EXTERNAL RES, PSOL
DIMENSION Y(*), YPRIME(*), P(*), WT(*), SAVR(*), R(*), ID(*)
DIMENSION WM(*), IWM(*), YNEW(*), YPNEW(*), PWK(*), ICNSTR(*)
DIMENSION WP(*), IWP(*), RPAR(*), IPAR(*)
CHARACTER MSG*80
C
PARAMETER (LNRE=12, LNPS=21, LKPRIN=31)
C
SAVE ALPHA, ONE, TWO
DATA ALPHA/1.0D-4/, ONE/1.0D0/, TWO/2.0D0/
C
KPRIN=IWM(LKPRIN)
F1NRM = (FNRM*FNRM)/TWO
RATIO = ONE
C
IF (KPRIN .GE. 2) THEN
MSG = '------ IN ROUTINE DLINSK-- PNRM = (R1)'
CALL XERRWD(MSG, 38, 921, 0, 0, 0, 0, 1, PNRM, 0.0D0)
ENDIF
TAU = PNRM
RL = ONE
C-----------------------------------------------------------------------
C Check for violations of the constraints, if any are imposed.
C If any violations are found, the step vector P is rescaled, and the
C constraint check is repeated, until no violations are found.
C-----------------------------------------------------------------------
IF (ICNFLG .NE. 0) THEN
10 CONTINUE
CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW)
CALL DCNSTR (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR)
IF (IRET .EQ. 1) THEN
RATIO1 = TAU/PNRM
RATIO = RATIO*RATIO1
DO 20 I = 1,NEQ
20 P(I) = P(I)*RATIO1
PNRM = TAU
IF (KPRIN .GE. 2) THEN
MSG = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)'
CALL XERRWD(MSG, 50, 922, 0, 1, IVAR, 0, 1, PNRM, 0.0D0)
ENDIF
IF (PNRM .LE. STPTOL) THEN
IRET = 1
RETURN
ENDIF
GO TO 10
ENDIF
ENDIF
C
SLPI = -TWO*F1NRM*RATIO
RLMIN = STPTOL/PNRM
IF (LSOFF .EQ. 0 .AND. KPRIN .GE. 2) THEN
MSG = '------ MIN. LAMBDA = (R1)'
CALL XERRWD(MSG, 25, 923, 0, 0, 0, 0, 1, RLMIN, 0.0D0)
ENDIF
C-----------------------------------------------------------------------
C Begin iteration to find RL value satisfying alpha-condition.
C Update YNEW and YPNEW, then compute norm of new scaled residual and
C perform alpha condition test.
C-----------------------------------------------------------------------
100 CONTINUE
CALL DYYPNW (NEQ,Y,YPRIME,CJ,RL,P,ICOPT,ID,YNEW,YPNEW)
CALL DFNRMK (NEQ, YNEW, T, YPNEW, SAVR, R, CJ, TSCALE, WT,
* SQRTN, RSQRTN, RES, IRES, PSOL, 0, IER, FNRMP, EPLIN,
* WP, IWP, PWK, RPAR, IPAR)
IWM(LNRE) = IWM(LNRE) + 1
IF (IRES .GE. 0) IWM(LNPS) = IWM(LNPS) + 1
IF (IRES .NE. 0 .OR. IER .NE. 0) THEN
IRET = 2
RETURN
ENDIF
IF (LSOFF .EQ. 1) GO TO 150
C
F1NRMP = FNRMP*FNRMP/TWO
IF (KPRIN .GE. 2) THEN
MSG = '------ LAMBDA = (R1)'
CALL XERRWD(MSG, 20, 924, 0, 0, 0, 0, 1, RL, 0.0D0)
MSG = '------ NORM(F1) = (R1), NORM(F1NEW) = (R2)'
CALL XERRWD(MSG, 43, 925, 0, 0, 0, 0, 2, F1NRM, F1NRMP)
ENDIF
IF (F1NRMP .GT. F1NRM + ALPHA*SLPI*RL) GO TO 200
C-----------------------------------------------------------------------
C Alpha-condition is satisfied, or linesearch is turned off.
C Copy YNEW,YPNEW to Y,YPRIME and return.
C-----------------------------------------------------------------------
150 IRET = 0
CALL DCOPY(NEQ, YNEW, 1, Y, 1)
CALL DCOPY(NEQ, YPNEW, 1, YPRIME, 1)
FNRM = FNRMP
IF (KPRIN .GE. 1) THEN
MSG = '------ LEAVING ROUTINE DLINSK, FNRM = (R1)'
CALL XERRWD(MSG, 42, 926, 0, 0, 0, 0, 1, FNRM, 0.0D0)
ENDIF
RETURN
C-----------------------------------------------------------------------
C Alpha-condition not satisfied. Perform backtrack to compute new RL
C value. If RL is less than RLMIN, i.e. no satisfactory YNEW,YPNEW can
C be found sufficiently distinct from Y,YPRIME, then return IRET = 1.
C-----------------------------------------------------------------------
200 CONTINUE
IF (RL .LT. RLMIN) THEN
IRET = 1
RETURN
ENDIF
C
RL = RL/TWO
GO TO 100
C
C----------------------- END OF SUBROUTINE DLINSK ----------------------
END
SUBROUTINE DFNRMK (NEQ, Y, T, YPRIME, SAVR, R, CJ, TSCALE, WT,
* SQRTN, RSQRTN, RES, IRES, PSOL, IRIN, IER,
* FNORM, EPLIN, WP, IWP, PWK, RPAR, IPAR)
C
C***BEGIN PROLOGUE DFNRMK
C***REFER TO DLINSK
C***DATE WRITTEN 940830 (YYMMDD)
C***REVISION DATE 951006 (SQRTN, RSQRTN, and scaling of WT added.)
C***REVISION DATE 000628 TSCALE argument added.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DFNRMK calculates the scaled preconditioned norm of the nonlinear
C function used in the nonlinear iteration for obtaining consistent
C initial conditions. Specifically, DFNRMK calculates the weighted
C root-mean-square norm of the vector (P-inverse)*G(T,Y,YPRIME),
C where P is the preconditioner matrix.
C
C In addition to the parameters described in the calling program
C DLINSK, the parameters represent
C
C TSCALE -- Scale factor in T, used for stopping tests if nonzero.
C IRIN -- Flag showing whether the current residual vector is
C input in SAVR. 1 means it is, 0 means it is not.
C R -- Array of length NEQ that contains
C (P-inverse)*G(T,Y,YPRIME) on return.
C FNORM -- Scalar containing the weighted norm of R on return.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C RES, DCOPY, DSCAL, PSOL, DDWNRM
C
C***END PROLOGUE DFNRMK
C
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
EXTERNAL RES, PSOL
DIMENSION Y(*), YPRIME(*), WT(*), SAVR(*), R(*), PWK(*)
DIMENSION WP(*), IWP(*), RPAR(*), IPAR(*)
C-----------------------------------------------------------------------
C Call RES routine if IRIN = 0.
C-----------------------------------------------------------------------
IF (IRIN .EQ. 0) THEN
IRES = 0
CALL RES (T, Y, YPRIME, CJ, SAVR, IRES, RPAR, IPAR)
IF (IRES .LT. 0) RETURN
ENDIF
C-----------------------------------------------------------------------
C Apply inverse of left preconditioner to vector R.
C First scale WT array by 1/sqrt(N), and undo scaling afterward.
C-----------------------------------------------------------------------
CALL DCOPY(NEQ, SAVR, 1, R, 1)
CALL DSCAL (NEQ, RSQRTN, WT, 1)
IER = 0
CALL PSOL (NEQ, T, Y, YPRIME, SAVR, PWK, CJ, WT, WP, IWP,
* R, EPLIN, IER, RPAR, IPAR)
CALL DSCAL (NEQ, SQRTN, WT, 1)
IF (IER .NE. 0) RETURN
C-----------------------------------------------------------------------
C Calculate norm of R.
C-----------------------------------------------------------------------
FNORM = DDWNRM (NEQ, R, WT, RPAR, IPAR)
IF (TSCALE .GT. 0.0D0) FNORM = FNORM*TSCALE*ABS(CJ)
C
RETURN
C----------------------- END OF SUBROUTINE DFNRMK ----------------------
END
SUBROUTINE DNEDK(X,Y,YPRIME,NEQ,RES,JACK,PSOL,
* H,WT,JSTART,IDID,RPAR,IPAR,PHI,GAMMA,SAVR,DELTA,E,
* WM,IWM,CJ,CJOLD,CJLAST,S,UROUND,EPLI,SQRTN,RSQRTN,
* EPCON,JCALC,JFLG,KP1,NONNEG,NTYPE,IERNLS)
C
C***BEGIN PROLOGUE DNEDK
C***REFER TO DDASPK
C***DATE WRITTEN 891219 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***REVISION DATE 940701 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DNEDK solves a nonlinear system of
C algebraic equations of the form
C G(X,Y,YPRIME) = 0 for the unknown Y.
C
C The method used is a matrix-free Newton scheme.
C
C The parameters represent
C X -- Independent variable.
C Y -- Solution vector at x.
C YPRIME -- Derivative of solution vector
C after successful step.
C NEQ -- Number of equations to be integrated.
C RES -- External user-supplied subroutine
C to evaluate the residual. See RES description
C in DDASPK prologue.
C JACK -- External user-supplied routine to update
C the preconditioner. (This is optional).
C See JAC description for the case
C INFO(12) = 1 in the DDASPK prologue.
C PSOL -- External user-supplied routine to solve
C a linear system using preconditioning.
C (This is optional). See explanation inside DDASPK.
C H -- Appropriate step size for this step.
C WT -- Vector of weights for error criterion.
C JSTART -- Indicates first call to this routine.
C If JSTART = 0, then this is the first call,
C otherwise it is not.
C IDID -- Completion flag, output by DNEDK.
C See IDID description in DDASPK prologue.
C RPAR,IPAR -- Real and integer arrays used for communication
C between the calling program and external user
C routines. They are not altered within DASPK.
C PHI -- Array of divided differences used by
C DNEDK. The length is NEQ*(K+1), where
C K is the maximum order.
C GAMMA -- Array used to predict Y and YPRIME. The length
C is K+1, where K is the maximum order.
C SAVR -- Work vector for DNEDK of length NEQ.
C DELTA -- Work vector for DNEDK of length NEQ.
C E -- Error accumulation vector for DNEDK of length NEQ.
C WM,IWM -- Real and integer arrays storing
C matrix information for linear system
C solvers, and various other information.
C CJ -- Parameter always proportional to 1/H.
C CJOLD -- Saves the value of CJ as of the last call to DITMD.
C Accounts for changes in CJ needed to
C decide whether to call DITMD.
C CJLAST -- Previous value of CJ.
C S -- A scalar determined by the approximate rate
C of convergence of the Newton iteration and used
C in the convergence test for the Newton iteration.
C
C If RATE is defined to be an estimate of the
C rate of convergence of the Newton iteration,
C then S = RATE/(1.D0-RATE).
C
C The closer RATE is to 0., the faster the Newton
C iteration is converging; the closer RATE is to 1.,
C the slower the Newton iteration is converging.
C
C On the first Newton iteration with an up-dated
C preconditioner S = 100.D0, Thus the initial
C RATE of convergence is approximately 1.
C
C S is preserved from call to call so that the rate
C estimate from a previous step can be applied to
C the current step.
C UROUND -- Unit roundoff. Not used here.
C EPLI -- convergence test constant.
C See DDASPK prologue for more details.
C SQRTN -- Square root of NEQ.
C RSQRTN -- reciprical of square root of NEQ.
C EPCON -- Tolerance to test for convergence of the Newton
C iteration.
C JCALC -- Flag used to determine when to update
C the Jacobian matrix. In general:
C
C JCALC = -1 ==> Call the DITMD routine to update
C the Jacobian matrix.
C JCALC = 0 ==> Jacobian matrix is up-to-date.
C JCALC = 1 ==> Jacobian matrix is out-dated,
C but DITMD will not be called unless
C JCALC is set to -1.
C JFLG -- Flag showing whether a Jacobian routine is supplied.
C KP1 -- The current order + 1; updated across calls.
C NONNEG -- Flag to determine nonnegativity constraints.
C NTYPE -- Identification code for the DNEDK routine.
C 1 ==> modified Newton; iterative linear solver.
C 2 ==> modified Newton; user-supplied linear solver.
C IERNLS -- Error flag for nonlinear solver.
C 0 ==> nonlinear solver converged.
C 1 ==> recoverable error inside non-linear solver.
C -1 ==> unrecoverable error inside non-linear solver.
C
C The following group of variables are passed as arguments to
C the Newton iteration solver. They are explained in greater detail
C in DNSK:
C TOLNEW, MULDEL, MAXIT, IERNEW
C
C IERTYP -- Flag which tells whether this subroutine is correct.
C 0 ==> correct subroutine.
C 1 ==> incorrect subroutine.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C RES, JACK, DDWNRM, DNSK
C
C***END PROLOGUE DNEDK
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),WT(*)
DIMENSION PHI(NEQ,*),SAVR(*),DELTA(*),E(*)
DIMENSION WM(*),IWM(*)
DIMENSION GAMMA(*),RPAR(*),IPAR(*)
EXTERNAL RES, JACK, PSOL
C
PARAMETER (LNRE=12, LNJE=13, LLOCWP=29, LLCIWP=30)
C
SAVE MULDEL, MAXIT, XRATE
DATA MULDEL/0/, MAXIT/4/, XRATE/0.25D0/
C
C Verify that this is the correct subroutine.
C
IERTYP = 0
IF (NTYPE .NE. 1) THEN
IERTYP = 1
GO TO 380
ENDIF
C
C If this is the first step, perform initializations.
C
IF (JSTART .EQ. 0) THEN
CJOLD = CJ
JCALC = -1
S = 100.D0
ENDIF
C
C Perform all other initializations.
C
IERNLS = 0
LWP = IWM(LLOCWP)
LIWP = IWM(LLCIWP)
C
C Decide whether to update the preconditioner.
C
IF (JFLG .NE. 0) THEN
TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE)
TEMP2 = 1.0D0/TEMP1
IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1
IF (CJ .NE. CJLAST) S = 100.D0
ELSE
JCALC = 0
ENDIF
C
C Looping point for updating preconditioner with current stepsize.
C
300 CONTINUE
C
C Initialize all error flags to zero.
C
IERPJ = 0
IRES = 0
IERSL = 0
IERNEW = 0
C
C Predict the solution and derivative and compute the tolerance
C for the Newton iteration.
C
DO 310 I=1,NEQ
Y(I)=PHI(I,1)
310 YPRIME(I)=0.0D0
DO 330 J=2,KP1
DO 320 I=1,NEQ
Y(I)=Y(I)+PHI(I,J)
320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J)
330 CONTINUE
EPLIN = EPLI*EPCON
TOLNEW = EPLIN
C
C Call RES to initialize DELTA.
C
IWM(LNRE)=IWM(LNRE)+1
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
IF (IRES .LT. 0) GO TO 380
C
C
C If indicated, update the preconditioner.
C Set JCALC to 0 as an indicator that this has been done.
C
IF(JCALC .EQ. -1)THEN
IWM(LNJE) = IWM(LNJE) + 1
JCALC=0
CALL JACK (RES, IRES, NEQ, X, Y, YPRIME, WT, DELTA, E, H, CJ,
* WM(LWP), IWM(LIWP), IERPJ, RPAR, IPAR)
CJOLD=CJ
S = 100.D0
IF (IRES .LT. 0) GO TO 380
IF (IERPJ .NE. 0) GO TO 380
ENDIF
C
C Call the nonlinear Newton solver.
C
CALL DNSK(X,Y,YPRIME,NEQ,RES,PSOL,WT,RPAR,IPAR,SAVR,
* DELTA,E,WM,IWM,CJ,SQRTN,RSQRTN,EPLIN,EPCON,
* S,TEMP1,TOLNEW,MULDEL,MAXIT,IRES,IERSL,IERNEW)
C
IF (IERNEW .GT. 0 .AND. JCALC .NE. 0) THEN
C
C The Newton iteration had a recoverable failure with an old
C preconditioner. Retry the step with a new preconditioner.
C
JCALC = -1
GO TO 300
ENDIF
C
IF (IERNEW .NE. 0) GO TO 380
C
C The Newton iteration has converged. If nonnegativity of
C solution is required, set the solution nonnegative, if the
C perturbation to do it is small enough. If the change is too
C large, then consider the corrector iteration to have failed.
C
IF(NONNEG .EQ. 0) GO TO 390
DO 360 I = 1,NEQ
360 DELTA(I) = MIN(Y(I),0.0D0)
DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
IF(DELNRM .GT. EPCON) GO TO 380
DO 370 I = 1,NEQ
370 E(I) = E(I) - DELTA(I)
GO TO 390
C
C
C Exits from nonlinear solver.
C No convergence with current preconditioner.
C Compute IERNLS and IDID accordingly.
C
380 CONTINUE
IF (IRES .LE. -2 .OR. IERSL .LT. 0 .OR. IERTYP .NE. 0) THEN
IERNLS = -1
IF (IRES .LE. -2) IDID = -11
IF (IERSL .LT. 0) IDID = -13
IF (IERTYP .NE. 0) IDID = -15
ELSE
IERNLS = 1
IF (IRES .EQ. -1) IDID = -10
IF (IERPJ .NE. 0) IDID = -5
IF (IERSL .GT. 0) IDID = -14
ENDIF
C
C
390 JCALC = 1
RETURN
C
C------END OF SUBROUTINE DNEDK------------------------------------------
END
SUBROUTINE DNSK(X,Y,YPRIME,NEQ,RES,PSOL,WT,RPAR,IPAR,
* SAVR,DELTA,E,WM,IWM,CJ,SQRTN,RSQRTN,EPLIN,EPCON,
* S,CONFAC,TOLNEW,MULDEL,MAXIT,IRES,IERSL,IERNEW)
C
C***BEGIN PROLOGUE DNSK
C***REFER TO DDASPK
C***DATE WRITTEN 891219 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***REVISION DATE 950126 (YYMMDD)
C***REVISION DATE 000711 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DNSK solves a nonlinear system of
C algebraic equations of the form
C G(X,Y,YPRIME) = 0 for the unknown Y.
C
C The method used is a modified Newton scheme.
C
C The parameters represent
C
C X -- Independent variable.
C Y -- Solution vector.
C YPRIME -- Derivative of solution vector.
C NEQ -- Number of unknowns.
C RES -- External user-supplied subroutine
C to evaluate the residual. See RES description
C in DDASPK prologue.
C PSOL -- External user-supplied routine to solve
C a linear system using preconditioning.
C See explanation inside DDASPK.
C WT -- Vector of weights for error criterion.
C RPAR,IPAR -- Real and integer arrays used for communication
C between the calling program and external user
C routines. They are not altered within DASPK.
C SAVR -- Work vector for DNSK of length NEQ.
C DELTA -- Work vector for DNSK of length NEQ.
C E -- Error accumulation vector for DNSK of length NEQ.
C WM,IWM -- Real and integer arrays storing
C matrix information such as the matrix
C of partial derivatives, permutation
C vector, and various other information.
C CJ -- Parameter always proportional to 1/H (step size).
C SQRTN -- Square root of NEQ.
C RSQRTN -- reciprical of square root of NEQ.
C EPLIN -- Tolerance for linear system solver.
C EPCON -- Tolerance to test for convergence of the Newton
C iteration.
C S -- Used for error convergence tests.
C In the Newton iteration: S = RATE/(1.D0-RATE),
C where RATE is the estimated rate of convergence
C of the Newton iteration.
C
C The closer RATE is to 0., the faster the Newton
C iteration is converging; the closer RATE is to 1.,
C the slower the Newton iteration is converging.
C
C The calling routine sends the initial value
C of S to the Newton iteration.
C CONFAC -- A residual scale factor to improve convergence.
C TOLNEW -- Tolerance on the norm of Newton correction in
C alternative Newton convergence test.
C MULDEL -- A flag indicating whether or not to multiply
C DELTA by CONFAC.
C 0 ==> do not scale DELTA by CONFAC.
C 1 ==> scale DELTA by CONFAC.
C MAXIT -- Maximum allowed number of Newton iterations.
C IRES -- Error flag returned from RES. See RES description
C in DDASPK prologue. If IRES = -1, then IERNEW
C will be set to 1.
C If IRES < -1, then IERNEW will be set to -1.
C IERSL -- Error flag for linear system solver.
C See IERSL description in subroutine DSLVK.
C If IERSL = 1, then IERNEW will be set to 1.
C If IERSL < 0, then IERNEW will be set to -1.
C IERNEW -- Error flag for Newton iteration.
C 0 ==> Newton iteration converged.
C 1 ==> recoverable error inside Newton iteration.
C -1 ==> unrecoverable error inside Newton iteration.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C RES, DSLVK, DDWNRM
C
C***END PROLOGUE DNSK
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),WT(*),DELTA(*),E(*),SAVR(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
EXTERNAL RES, PSOL
C
PARAMETER (LNNI=19, LNRE=12)
C
C Initialize Newton counter M and accumulation vector E.
C
M = 0
DO 100 I=1,NEQ
100 E(I) = 0.0D0
C
C Corrector loop.
C
300 CONTINUE
IWM(LNNI) = IWM(LNNI) + 1
C
C If necessary, multiply residual by convergence factor.
C
IF (MULDEL .EQ. 1) THEN
DO 320 I = 1,NEQ
320 DELTA(I) = DELTA(I) * CONFAC
ENDIF
C
C Save residual in SAVR.
C
DO 340 I = 1,NEQ
340 SAVR(I) = DELTA(I)
C
C Compute a new iterate. Store the correction in DELTA.
C
CALL DSLVK (NEQ, Y, X, YPRIME, SAVR, DELTA, WT, WM, IWM,
* RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK,
* RPAR, IPAR)
IF (IRES .NE. 0 .OR. IERSL .NE. 0) GO TO 380
C
C Update Y, E, and YPRIME.
C
DO 360 I=1,NEQ
Y(I) = Y(I) - DELTA(I)
E(I) = E(I) - DELTA(I)
360 YPRIME(I) = YPRIME(I) - CJ*DELTA(I)
C
C Test for convergence of the iteration.
C
DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
IF (M .EQ. 0) THEN
OLDNRM = DELNRM
IF (DELNRM .LE. TOLNEW) GO TO 370
ELSE
RATE = (DELNRM/OLDNRM)**(1.0D0/M)
IF (RATE .GT. 0.9D0) GO TO 380
S = RATE/(1.0D0 - RATE)
ENDIF
IF (S*DELNRM .LE. EPCON) GO TO 370
C
C The corrector has not yet converged. Update M and test whether
C the maximum number of iterations have been tried.
C
M = M + 1
IF (M .GE. MAXIT) GO TO 380
C
C Evaluate the residual, and go back to do another iteration.
C
IWM(LNRE) = IWM(LNRE) + 1
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
IF (IRES .LT. 0) GO TO 380
GO TO 300
C
C The iteration has converged.
C
370 RETURN
C
C The iteration has not converged. Set IERNEW appropriately.
C
380 CONTINUE
IF (IRES .LE. -2 .OR. IERSL .LT. 0) THEN
IERNEW = -1
ELSE
IERNEW = 1
ENDIF
RETURN
C
C
C------END OF SUBROUTINE DNSK-------------------------------------------
END
SUBROUTINE DSLVK (NEQ, Y, TN, YPRIME, SAVR, X, EWT, WM, IWM,
* RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK,
* RPAR, IPAR)
C
C***BEGIN PROLOGUE DSLVK
C***REFER TO DDASPK
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***REVISION DATE 940928 Removed MNEWT and added RHOK in call list.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DSLVK uses a restart algorithm and interfaces to DSPIGM for
C the solution of the linear system arising from a Newton iteration.
C
C In addition to variables described elsewhere,
C communication with DSLVK uses the following variables..
C WM = Real work space containing data for the algorithm
C (Krylov basis vectors, Hessenberg matrix, etc.).
C IWM = Integer work space containing data for the algorithm.
C X = The right-hand side vector on input, and the solution vector
C on output, of length NEQ.
C IRES = Error flag from RES.
C IERSL = Output flag ..
C IERSL = 0 means no trouble occurred (or user RES routine
C returned IRES < 0)
C IERSL = 1 means the iterative method failed to converge
C (DSPIGM returned IFLAG > 0.)
C IERSL = -1 means there was a nonrecoverable error in the
C iterative solver, and an error exit will occur.
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C DSCAL, DCOPY, DSPIGM
C
C***END PROLOGUE DSLVK
C
INTEGER NEQ, IWM, IRES, IERSL, IPAR
DOUBLE PRECISION Y, TN, YPRIME, SAVR, X, EWT, WM, CJ, EPLIN,
1 SQRTN, RSQRTN, RHOK, RPAR
DIMENSION Y(*), YPRIME(*), SAVR(*), X(*), EWT(*),
1 WM(*), IWM(*), RPAR(*), IPAR(*)
C
INTEGER IFLAG, IRST, NRSTS, NRMAX, LR, LDL, LHES, LGMR, LQ, LV,
1 LWK, LZ, MAXLP1, NPSL
INTEGER NLI, NPS, NCFL, NRE, MAXL, KMP, MITER
EXTERNAL RES, PSOL
C
PARAMETER (LNRE=12, LNCFL=16, LNLI=20, LNPS=21)
PARAMETER (LLOCWP=29, LLCIWP=30)
PARAMETER (LMITER=23, LMAXL=24, LKMP=25, LNRMAX=26)
C
C-----------------------------------------------------------------------
C IRST is set to 1, to indicate restarting is in effect.
C NRMAX is the maximum number of restarts.
C-----------------------------------------------------------------------
DATA IRST/1/
C
LIWP = IWM(LLCIWP)
NLI = IWM(LNLI)
NPS = IWM(LNPS)
NCFL = IWM(LNCFL)
NRE = IWM(LNRE)
LWP = IWM(LLOCWP)
MAXL = IWM(LMAXL)
KMP = IWM(LKMP)
NRMAX = IWM(LNRMAX)
MITER = IWM(LMITER)
IERSL = 0
IRES = 0
C-----------------------------------------------------------------------
C Use a restarting strategy to solve the linear system
C P*X = -F. Parse the work vector, and perform initializations.
C Note that zero is the initial guess for X.
C-----------------------------------------------------------------------
MAXLP1 = MAXL + 1
LV = 1
LR = LV + NEQ*MAXL
LHES = LR + NEQ + 1
LQ = LHES + MAXL*MAXLP1
LWK = LQ + 2*MAXL
LDL = LWK + MIN0(1,MAXL-KMP)*NEQ
LZ = LDL + NEQ
CALL DSCAL (NEQ, RSQRTN, EWT, 1)
CALL DCOPY (NEQ, X, 1, WM(LR), 1)
DO 110 I = 1,NEQ
110 X(I) = 0.D0
C-----------------------------------------------------------------------
C Top of loop for the restart algorithm. Initial pass approximates
C X and sets up a transformed system to perform subsequent restarts
C to update X. NRSTS is initialized to -1, because restarting
C does not occur until after the first pass.
C Update NRSTS; conditionally copy DL to R; call the DSPIGM
C algorithm to solve A*Z = R; updated counters; update X with
C the residual solution.
C Note: if convergence is not achieved after NRMAX restarts,
C then the linear solver is considered to have failed.
C-----------------------------------------------------------------------
NRSTS = -1
115 CONTINUE
NRSTS = NRSTS + 1
IF (NRSTS .GT. 0) CALL DCOPY (NEQ, WM(LDL), 1, WM(LR),1)
CALL DSPIGM (NEQ, TN, Y, YPRIME, SAVR, WM(LR), EWT, MAXL, MAXLP1,
1 KMP, EPLIN, CJ, RES, IRES, NRES, PSOL, NPSL, WM(LZ), WM(LV),
2 WM(LHES), WM(LQ), LGMR, WM(LWP), IWM(LIWP), WM(LWK),
3 WM(LDL), RHOK, IFLAG, IRST, NRSTS, RPAR, IPAR)
NLI = NLI + LGMR
NPS = NPS + NPSL
NRE = NRE + NRES
DO 120 I = 1,NEQ
120 X(I) = X(I) + WM(LZ+I-1)
IF ((IFLAG .EQ. 1) .AND. (NRSTS .LT. NRMAX) .AND. (IRES .EQ. 0))
1 GO TO 115
C-----------------------------------------------------------------------
C The restart scheme is finished. Test IRES and IFLAG to see if
C convergence was not achieved, and set flags accordingly.
C-----------------------------------------------------------------------
IF (IRES .LT. 0) THEN
NCFL = NCFL + 1
ELSE IF (IFLAG .NE. 0) THEN
NCFL = NCFL + 1
IF (IFLAG .GT. 0) IERSL = 1
IF (IFLAG .LT. 0) IERSL = -1
ENDIF
C-----------------------------------------------------------------------
C Update IWM with counters, rescale EWT, and return.
C-----------------------------------------------------------------------
IWM(LNLI) = NLI
IWM(LNPS) = NPS
IWM(LNCFL) = NCFL
IWM(LNRE) = NRE
CALL DSCAL (NEQ, SQRTN, EWT, 1)
RETURN
C
C------END OF SUBROUTINE DSLVK------------------------------------------
END
SUBROUTINE DSPIGM (NEQ, TN, Y, YPRIME, SAVR, R, WGHT, MAXL,
* MAXLP1, KMP, EPLIN, CJ, RES, IRES, NRE, PSOL, NPSL, Z, V,
* HES, Q, LGMR, WP, IWP, WK, DL, RHOK, IFLAG, IRST, NRSTS,
* RPAR, IPAR)
C
C***BEGIN PROLOGUE DSPIGM
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***REVISION DATE 940927 Removed MNEWT and added RHOK in call list.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This routine solves the linear system A * Z = R using a scaled
C preconditioned version of the generalized minimum residual method.
C An initial guess of Z = 0 is assumed.
C
C On entry
C
C NEQ = Problem size, passed to PSOL.
C
C TN = Current Value of T.
C
C Y = Array Containing current dependent variable vector.
C
C YPRIME = Array Containing current first derivative of Y.
C
C SAVR = Array containing current value of G(T,Y,YPRIME).
C
C R = The right hand side of the system A*Z = R.
C R is also used as work space when computing
C the final approximation and will therefore be
C destroyed.
C (R is the same as V(*,MAXL+1) in the call to DSPIGM.)
C
C WGHT = The vector of length NEQ containing the nonzero
C elements of the diagonal scaling matrix.
C
C MAXL = The maximum allowable order of the matrix H.
C
C MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES.
C
C KMP = The number of previous vectors the new vector, VNEW,
C must be made orthogonal to. (KMP .LE. MAXL.)
C
C EPLIN = Tolerance on residuals R-A*Z in weighted rms norm.
C
C CJ = Scalar proportional to current value of
C 1/(step size H).
C
C WK = Real work array used by routine DATV and PSOL.
C
C DL = Real work array used for calculation of the residual
C norm RHO when the method is incomplete (KMP.LT.MAXL)
C and/or when using restarting.
C
C WP = Real work array used by preconditioner PSOL.
C
C IWP = Integer work array used by preconditioner PSOL.
C
C IRST = Method flag indicating if restarting is being
C performed. IRST .GT. 0 means restarting is active,
C while IRST = 0 means restarting is not being used.
C
C NRSTS = Counter for the number of restarts on the current
C call to DSPIGM. If NRSTS .GT. 0, then the residual
C R is already scaled, and so scaling of R is not
C necessary.
C
C
C On Return
C
C Z = The final computed approximation to the solution
C of the system A*Z = R.
C
C LGMR = The number of iterations performed and
C the current order of the upper Hessenberg
C matrix HES.
C
C NRE = The number of calls to RES (i.e. DATV)
C
C NPSL = The number of calls to PSOL.
C
C V = The neq by (LGMR+1) array containing the LGMR
C orthogonal vectors V(*,1) to V(*,LGMR).
C
C HES = The upper triangular factor of the QR decomposition
C of the (LGMR+1) by LGMR upper Hessenberg matrix whose
C entries are the scaled inner-products of A*V(*,I)
C and V(*,K).
C
C Q = Real array of length 2*MAXL containing the components
C of the givens rotations used in the QR decomposition
C of HES. It is loaded in DHEQR and used in DHELS.
C
C IRES = Error flag from RES.
C
C DL = Scaled preconditioned residual,
C (D-inverse)*(P-inverse)*(R-A*Z). Only loaded when
C performing restarts of the Krylov iteration.
C
C RHOK = Weighted norm of final preconditioned residual.
C
C IFLAG = Integer error flag..
C 0 Means convergence in LGMR iterations, LGMR.LE.MAXL.
C 1 Means the convergence test did not pass in MAXL
C iterations, but the new residual norm (RHO) is
C .LT. the old residual norm (RNRM), and so Z is
C computed.
C 2 Means the convergence test did not pass in MAXL
C iterations, new residual norm (RHO) .GE. old residual
C norm (RNRM), and the initial guess, Z = 0, is
C returned.
C 3 Means there was a recoverable error in PSOL
C caused by the preconditioner being out of date.
C -1 Means there was an unrecoverable error in PSOL.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C PSOL, DNRM2, DSCAL, DATV, DORTH, DHEQR, DCOPY, DHELS, DAXPY
C
C***END PROLOGUE DSPIGM
C
INTEGER NEQ,MAXL,MAXLP1,KMP,IRES,NRE,NPSL,LGMR,IWP,
1 IFLAG,IRST,NRSTS,IPAR
DOUBLE PRECISION TN,Y,YPRIME,SAVR,R,WGHT,EPLIN,CJ,Z,V,HES,Q,WP,WK,
1 DL,RHOK,RPAR
DIMENSION Y(*), YPRIME(*), SAVR(*), R(*), WGHT(*), Z(*),
1 V(NEQ,*), HES(MAXLP1,*), Q(*), WP(*), IWP(*), WK(*), DL(*),
2 RPAR(*), IPAR(*)
INTEGER I, IER, INFO, IP1, I2, J, K, LL, LLP1
DOUBLE PRECISION RNRM,C,DLNRM,PROD,RHO,S,SNORMW,DNRM2,TEM
EXTERNAL RES, PSOL
C
IER = 0
IFLAG = 0
LGMR = 0
NPSL = 0
NRE = 0
C-----------------------------------------------------------------------
C The initial guess for Z is 0. The initial residual is therefore
C the vector R. Initialize Z to 0.
C-----------------------------------------------------------------------
DO 10 I = 1,NEQ
10 Z(I) = 0.0D0
C-----------------------------------------------------------------------
C Apply inverse of left preconditioner to vector R if NRSTS .EQ. 0.
C Form V(*,1), the scaled preconditioned right hand side.
C-----------------------------------------------------------------------
IF (NRSTS .EQ. 0) THEN
CALL PSOL (NEQ, TN, Y, YPRIME, SAVR, WK, CJ, WGHT, WP, IWP,
1 R, EPLIN, IER, RPAR, IPAR)
NPSL = 1
IF (IER .NE. 0) GO TO 300
DO 30 I = 1,NEQ
30 V(I,1) = R(I)*WGHT(I)
ELSE
DO 35 I = 1,NEQ
35 V(I,1) = R(I)
ENDIF
C-----------------------------------------------------------------------
C Calculate norm of scaled vector V(*,1) and normalize it
C If, however, the norm of V(*,1) (i.e. the norm of the preconditioned
C residual) is .le. EPLIN, then return with Z=0.
C-----------------------------------------------------------------------
RNRM = DNRM2 (NEQ, V, 1)
IF (RNRM .LE. EPLIN) THEN
RHOK = RNRM
RETURN
ENDIF
TEM = 1.0D0/RNRM
CALL DSCAL (NEQ, TEM, V(1,1), 1)
C-----------------------------------------------------------------------
C Zero out the HES array.
C-----------------------------------------------------------------------
DO 65 J = 1,MAXL
DO 60 I = 1,MAXLP1
60 HES(I,J) = 0.0D0
65 CONTINUE
C-----------------------------------------------------------------------
C Main loop to compute the vectors V(*,2) to V(*,MAXL).
C The running product PROD is needed for the convergence test.
C-----------------------------------------------------------------------
PROD = 1.0D0
DO 90 LL = 1,MAXL
LGMR = LL
C-----------------------------------------------------------------------
C Call routine DATV to compute VNEW = ABAR*V(LL), where ABAR is
C the matrix A with scaling and inverse preconditioner factors applied.
C Call routine DORTH to orthogonalize the new vector VNEW = V(*,LL+1).
C call routine DHEQR to update the factors of HES.
C-----------------------------------------------------------------------
CALL DATV (NEQ, Y, TN, YPRIME, SAVR, V(1,LL), WGHT, Z,
1 RES, IRES, PSOL, V(1,LL+1), WK, WP, IWP, CJ, EPLIN,
1 IER, NRE, NPSL, RPAR, IPAR)
IF (IRES .LT. 0) RETURN
IF (IER .NE. 0) GO TO 300
CALL DORTH (V(1,LL+1), V, HES, NEQ, LL, MAXLP1, KMP, SNORMW)
HES(LL+1,LL) = SNORMW
CALL DHEQR (HES, MAXLP1, LL, Q, INFO, LL)
IF (INFO .EQ. LL) GO TO 120
C-----------------------------------------------------------------------
C Update RHO, the estimate of the norm of the residual R - A*ZL.
C If KMP .LT. MAXL, then the vectors V(*,1),...,V(*,LL+1) are not
C necessarily orthogonal for LL .GT. KMP. The vector DL must then
C be computed, and its norm used in the calculation of RHO.
C-----------------------------------------------------------------------
PROD = PROD*Q(2*LL)
RHO = ABS(PROD*RNRM)
IF ((LL.GT.KMP) .AND. (KMP.LT.MAXL)) THEN
IF (LL .EQ. KMP+1) THEN
CALL DCOPY (NEQ, V(1,1), 1, DL, 1)
DO 75 I = 1,KMP
IP1 = I + 1
I2 = I*2
S = Q(I2)
C = Q(I2-1)
DO 70 K = 1,NEQ
70 DL(K) = S*DL(K) + C*V(K,IP1)
75 CONTINUE
ENDIF
S = Q(2*LL)
C = Q(2*LL-1)/SNORMW
LLP1 = LL + 1
DO 80 K = 1,NEQ
80 DL(K) = S*DL(K) + C*V(K,LLP1)
DLNRM = DNRM2 (NEQ, DL, 1)
RHO = RHO*DLNRM
ENDIF
C-----------------------------------------------------------------------
C Test for convergence. If passed, compute approximation ZL.
C If failed and LL .LT. MAXL, then continue iterating.
C-----------------------------------------------------------------------
IF (RHO .LE. EPLIN) GO TO 200
IF (LL .EQ. MAXL) GO TO 100
C-----------------------------------------------------------------------
C Rescale so that the norm of V(1,LL+1) is one.
C-----------------------------------------------------------------------
TEM = 1.0D0/SNORMW
CALL DSCAL (NEQ, TEM, V(1,LL+1), 1)
90 CONTINUE
100 CONTINUE
IF (RHO .LT. RNRM) GO TO 150
120 CONTINUE
IFLAG = 2
DO 130 I = 1,NEQ
130 Z(I) = 0.D0
RETURN
150 IFLAG = 1
C-----------------------------------------------------------------------
C The tolerance was not met, but the residual norm was reduced.
C If performing restarting (IRST .gt. 0) calculate the residual vector
C RL and store it in the DL array. If the incomplete version is
C being used (KMP .lt. MAXL) then DL has already been calculated.
C-----------------------------------------------------------------------
IF (IRST .GT. 0) THEN
IF (KMP .EQ. MAXL) THEN
C
C Calculate DL from the V(I)'s.
C
CALL DCOPY (NEQ, V(1,1), 1, DL, 1)
MAXLM1 = MAXL - 1
DO 175 I = 1,MAXLM1
IP1 = I + 1
I2 = I*2
S = Q(I2)
C = Q(I2-1)
DO 170 K = 1,NEQ
170 DL(K) = S*DL(K) + C*V(K,IP1)
175 CONTINUE
S = Q(2*MAXL)
C = Q(2*MAXL-1)/SNORMW
DO 180 K = 1,NEQ
180 DL(K) = S*DL(K) + C*V(K,MAXLP1)
ENDIF
C
C Scale DL by RNRM*PROD to obtain the residual RL.
C
TEM = RNRM*PROD
CALL DSCAL(NEQ, TEM, DL, 1)
ENDIF
C-----------------------------------------------------------------------
C Compute the approximation ZL to the solution.
C Since the vector Z was used as work space, and the initial guess
C of the Newton correction is zero, Z must be reset to zero.
C-----------------------------------------------------------------------
200 CONTINUE
LL = LGMR
LLP1 = LL + 1
DO 210 K = 1,LLP1
210 R(K) = 0.0D0
R(1) = RNRM
CALL DHELS (HES, MAXLP1, LL, Q, R)
DO 220 K = 1,NEQ
220 Z(K) = 0.0D0
DO 230 I = 1,LL
CALL DAXPY (NEQ, R(I), V(1,I), 1, Z, 1)
230 CONTINUE
DO 240 I = 1,NEQ
240 Z(I) = Z(I)/WGHT(I)
C Load RHO into RHOK.
RHOK = RHO
RETURN
C-----------------------------------------------------------------------
C This block handles error returns forced by routine PSOL.
C-----------------------------------------------------------------------
300 CONTINUE
IF (IER .LT. 0) IFLAG = -1
IF (IER .GT. 0) IFLAG = 3
C
RETURN
C
C------END OF SUBROUTINE DSPIGM-----------------------------------------
END
SUBROUTINE DATV (NEQ, Y, TN, YPRIME, SAVR, V, WGHT, YPTEM, RES,
* IRES, PSOL, Z, VTEM, WP, IWP, CJ, EPLIN, IER, NRE, NPSL,
* RPAR,IPAR)
C
C***BEGIN PROLOGUE DATV
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This routine computes the product
C
C Z = (D-inverse)*(P-inverse)*(dF/dY)*(D*V),
C
C where F(Y) = G(T, Y, CJ*(Y-A)), CJ is a scalar proportional to 1/H,
C and A involves the past history of Y. The quantity CJ*(Y-A) is
C an approximation to the first derivative of Y and is stored
C in the array YPRIME. Note that dF/dY = dG/dY + CJ*dG/dYPRIME.
C
C D is a diagonal scaling matrix, and P is the left preconditioning
C matrix. V is assumed to have L2 norm equal to 1.
C The product is stored in Z and is computed by means of a
C difference quotient, a call to RES, and one call to PSOL.
C
C On entry
C
C NEQ = Problem size, passed to RES and PSOL.
C
C Y = Array containing current dependent variable vector.
C
C YPRIME = Array containing current first derivative of y.
C
C SAVR = Array containing current value of G(T,Y,YPRIME).
C
C V = Real array of length NEQ (can be the same array as Z).
C
C WGHT = Array of length NEQ containing scale factors.
C 1/WGHT(I) are the diagonal elements of the matrix D.
C
C YPTEM = Work array of length NEQ.
C
C VTEM = Work array of length NEQ used to store the
C unscaled version of V.
C
C WP = Real work array used by preconditioner PSOL.
C
C IWP = Integer work array used by preconditioner PSOL.
C
C CJ = Scalar proportional to current value of
C 1/(step size H).
C
C
C On return
C
C Z = Array of length NEQ containing desired scaled
C matrix-vector product.
C
C IRES = Error flag from RES.
C
C IER = Error flag from PSOL.
C
C NRE = The number of calls to RES.
C
C NPSL = The number of calls to PSOL.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C RES, PSOL
C
C***END PROLOGUE DATV
C
INTEGER NEQ, IRES, IWP, IER, NRE, NPSL, IPAR
DOUBLE PRECISION Y, TN, YPRIME, SAVR, V, WGHT, YPTEM, Z, VTEM,
1 WP, CJ, RPAR
DIMENSION Y(*), YPRIME(*), SAVR(*), V(*), WGHT(*), YPTEM(*),
1 Z(*), VTEM(*), WP(*), IWP(*), RPAR(*), IPAR(*)
INTEGER I
DOUBLE PRECISION EPLIN
EXTERNAL RES, PSOL
C
IRES = 0
C-----------------------------------------------------------------------
C Set VTEM = D * V.
C-----------------------------------------------------------------------
DO 10 I = 1,NEQ
10 VTEM(I) = V(I)/WGHT(I)
IER = 0
C-----------------------------------------------------------------------
C Store Y in Z and increment Z by VTEM.
C Store YPRIME in YPTEM and increment YPTEM by VTEM*CJ.
C-----------------------------------------------------------------------
DO 20 I = 1,NEQ
YPTEM(I) = YPRIME(I) + VTEM(I)*CJ
20 Z(I) = Y(I) + VTEM(I)
C-----------------------------------------------------------------------
C Call RES with incremented Y, YPRIME arguments
C stored in Z, YPTEM. VTEM is overwritten with new residual.
C-----------------------------------------------------------------------
CONTINUE
CALL RES(TN,Z,YPTEM,CJ,VTEM,IRES,RPAR,IPAR)
NRE = NRE + 1
IF (IRES .LT. 0) RETURN
C-----------------------------------------------------------------------
C Set Z = (dF/dY) * VBAR using difference quotient.
C (VBAR is old value of VTEM before calling RES)
C-----------------------------------------------------------------------
DO 70 I = 1,NEQ
70 Z(I) = VTEM(I) - SAVR(I)
C-----------------------------------------------------------------------
C Apply inverse of left preconditioner to Z.
C-----------------------------------------------------------------------
CALL PSOL (NEQ, TN, Y, YPRIME, SAVR, YPTEM, CJ, WGHT, WP, IWP,
1 Z, EPLIN, IER, RPAR, IPAR)
NPSL = NPSL + 1
IF (IER .NE. 0) RETURN
C-----------------------------------------------------------------------
C Apply D-inverse to Z and return.
C-----------------------------------------------------------------------
DO 90 I = 1,NEQ
90 Z(I) = Z(I)*WGHT(I)
RETURN
C
C------END OF SUBROUTINE DATV-------------------------------------------
END
SUBROUTINE DORTH (VNEW, V, HES, N, LL, LDHES, KMP, SNORMW)
C
C***BEGIN PROLOGUE DORTH
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This routine orthogonalizes the vector VNEW against the previous
C KMP vectors in the V array. It uses a modified Gram-Schmidt
C orthogonalization procedure with conditional reorthogonalization.
C
C On entry
C
C VNEW = The vector of length N containing a scaled product
C OF The Jacobian and the vector V(*,LL).
C
C V = The N x LL array containing the previous LL
C orthogonal vectors V(*,1) to V(*,LL).
C
C HES = An LL x LL upper Hessenberg matrix containing,
C in HES(I,K), K.LT.LL, scaled inner products of
C A*V(*,K) and V(*,I).
C
C LDHES = The leading dimension of the HES array.
C
C N = The order of the matrix A, and the length of VNEW.
C
C LL = The current order of the matrix HES.
C
C KMP = The number of previous vectors the new vector VNEW
C must be made orthogonal to (KMP .LE. MAXL).
C
C
C On return
C
C VNEW = The new vector orthogonal to V(*,I0),
C where I0 = MAX(1, LL-KMP+1).
C
C HES = Upper Hessenberg matrix with column LL filled in with
C scaled inner products of A*V(*,LL) and V(*,I).
C
C SNORMW = L-2 norm of VNEW.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C DDOT, DNRM2, DAXPY
C
C***END PROLOGUE DORTH
C
INTEGER N, LL, LDHES, KMP
DOUBLE PRECISION VNEW, V, HES, SNORMW
DIMENSION VNEW(*), V(N,*), HES(LDHES,*)
INTEGER I, I0
DOUBLE PRECISION ARG, DDOT, DNRM2, SUMDSQ, TEM, VNRM
C
C-----------------------------------------------------------------------
C Get norm of unaltered VNEW for later use.
C-----------------------------------------------------------------------
VNRM = DNRM2 (N, VNEW, 1)
C-----------------------------------------------------------------------
C Do Modified Gram-Schmidt on VNEW = A*V(LL).
C Scaled inner products give new column of HES.
C Projections of earlier vectors are subtracted from VNEW.
C-----------------------------------------------------------------------
I0 = MAX0(1,LL-KMP+1)
DO 10 I = I0,LL
HES(I,LL) = DDOT (N, V(1,I), 1, VNEW, 1)
TEM = -HES(I,LL)
CALL DAXPY (N, TEM, V(1,I), 1, VNEW, 1)
10 CONTINUE
C-----------------------------------------------------------------------
C Compute SNORMW = norm of VNEW.
C If VNEW is small compared to its input value (in norm), then
C Reorthogonalize VNEW to V(*,1) through V(*,LL).
C Correct if relative correction exceeds 1000*(unit roundoff).
C Finally, correct SNORMW using the dot products involved.
C-----------------------------------------------------------------------
SNORMW = DNRM2 (N, VNEW, 1)
IF (VNRM + 0.001D0*SNORMW .NE. VNRM) RETURN
SUMDSQ = 0.0D0
DO 30 I = I0,LL
TEM = -DDOT (N, V(1,I), 1, VNEW, 1)
IF (HES(I,LL) + 0.001D0*TEM .EQ. HES(I,LL)) GO TO 30
HES(I,LL) = HES(I,LL) - TEM
CALL DAXPY (N, TEM, V(1,I), 1, VNEW, 1)
SUMDSQ = SUMDSQ + TEM**2
30 CONTINUE
IF (SUMDSQ .EQ. 0.0D0) RETURN
ARG = MAX(0.0D0,SNORMW**2 - SUMDSQ)
SNORMW = SQRT(ARG)
RETURN
C
C------END OF SUBROUTINE DORTH------------------------------------------
END
SUBROUTINE DHEQR (A, LDA, N, Q, INFO, IJOB)
C
C***BEGIN PROLOGUE DHEQR
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This routine performs a QR decomposition of an upper
C Hessenberg matrix A. There are two options available:
C
C (1) performing a fresh decomposition
C (2) updating the QR factors by adding a row and A
C column to the matrix A.
C
C DHEQR decomposes an upper Hessenberg matrix by using Givens
C rotations.
C
C On entry
C
C A DOUBLE PRECISION(LDA, N)
C The matrix to be decomposed.
C
C LDA INTEGER
C The leading dimension of the array A.
C
C N INTEGER
C A is an (N+1) by N Hessenberg matrix.
C
C IJOB INTEGER
C = 1 Means that a fresh decomposition of the
C matrix A is desired.
C .GE. 2 Means that the current decomposition of A
C will be updated by the addition of a row
C and a column.
C On return
C
C A The upper triangular matrix R.
C The factorization can be written Q*A = R, where
C Q is a product of Givens rotations and R is upper
C triangular.
C
C Q DOUBLE PRECISION(2*N)
C The factors C and S of each Givens rotation used
C in decomposing A.
C
C INFO INTEGER
C = 0 normal value.
C = K If A(K,K) .EQ. 0.0. This is not an error
C condition for this subroutine, but it does
C indicate that DHELS will divide by zero
C if called.
C
C Modification of LINPACK.
C Peter Brown, Lawrence Livermore Natl. Lab.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED (NONE)
C
C***END PROLOGUE DHEQR
C
INTEGER LDA, N, INFO, IJOB
DOUBLE PRECISION A(LDA,*), Q(*)
INTEGER I, IQ, J, K, KM1, KP1, NM1
DOUBLE PRECISION C, S, T, T1, T2
C
IF (IJOB .GT. 1) GO TO 70
C-----------------------------------------------------------------------
C A new factorization is desired.
C-----------------------------------------------------------------------
C
C QR decomposition without pivoting.
C
INFO = 0
DO 60 K = 1, N
KM1 = K - 1
KP1 = K + 1
C
C Compute Kth column of R.
C First, multiply the Kth column of A by the previous
C K-1 Givens rotations.
C
IF (KM1 .LT. 1) GO TO 20
DO 10 J = 1, KM1
I = 2*(J-1) + 1
T1 = A(J,K)
T2 = A(J+1,K)
C = Q(I)
S = Q(I+1)
A(J,K) = C*T1 - S*T2
A(J+1,K) = S*T1 + C*T2
10 CONTINUE
C
C Compute Givens components C and S.
C
20 CONTINUE
IQ = 2*KM1 + 1
T1 = A(K,K)
T2 = A(KP1,K)
IF (T2 .NE. 0.0D0) GO TO 30
C = 1.0D0
S = 0.0D0
GO TO 50
30 CONTINUE
IF (ABS(T2) .LT. ABS(T1)) GO TO 40
T = T1/T2
S = -1.0D0/SQRT(1.0D0+T*T)
C = -S*T
GO TO 50
40 CONTINUE
T = T2/T1
C = 1.0D0/SQRT(1.0D0+T*T)
S = -C*T
50 CONTINUE
Q(IQ) = C
Q(IQ+1) = S
A(K,K) = C*T1 - S*T2
IF (A(K,K) .EQ. 0.0D0) INFO = K
60 CONTINUE
RETURN
C-----------------------------------------------------------------------
C The old factorization of A will be updated. A row and a column
C has been added to the matrix A.
C N by N-1 is now the old size of the matrix.
C-----------------------------------------------------------------------
70 CONTINUE
NM1 = N - 1
C-----------------------------------------------------------------------
C Multiply the new column by the N previous Givens rotations.
C-----------------------------------------------------------------------
DO 100 K = 1,NM1
I = 2*(K-1) + 1
T1 = A(K,N)
T2 = A(K+1,N)
C = Q(I)
S = Q(I+1)
A(K,N) = C*T1 - S*T2
A(K+1,N) = S*T1 + C*T2
100 CONTINUE
C-----------------------------------------------------------------------
C Complete update of decomposition by forming last Givens rotation,
C and multiplying it times the column vector (A(N,N),A(NP1,N)).
C-----------------------------------------------------------------------
INFO = 0
T1 = A(N,N)
T2 = A(N+1,N)
IF (T2 .NE. 0.0D0) GO TO 110
C = 1.0D0
S = 0.0D0
GO TO 130
110 CONTINUE
IF (ABS(T2) .LT. ABS(T1)) GO TO 120
T = T1/T2
S = -1.0D0/SQRT(1.0D0+T*T)
C = -S*T
GO TO 130
120 CONTINUE
T = T2/T1
C = 1.0D0/SQRT(1.0D0+T*T)
S = -C*T
130 CONTINUE
IQ = 2*N - 1
Q(IQ) = C
Q(IQ+1) = S
A(N,N) = C*T1 - S*T2
IF (A(N,N) .EQ. 0.0D0) INFO = N
RETURN
C
C------END OF SUBROUTINE DHEQR------------------------------------------
END
SUBROUTINE DHELS (A, LDA, N, Q, B)
C
C***BEGIN PROLOGUE DHELS
C***DATE WRITTEN 890101 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This is similar to the LINPACK routine DGESL except that
C A is an upper Hessenberg matrix.
C
C DHELS solves the least squares problem
C
C MIN (B-A*X,B-A*X)
C
C using the factors computed by DHEQR.
C
C On entry
C
C A DOUBLE PRECISION (LDA, N)
C The output from DHEQR which contains the upper
C triangular factor R in the QR decomposition of A.
C
C LDA INTEGER
C The leading dimension of the array A .
C
C N INTEGER
C A is originally an (N+1) by N matrix.
C
C Q DOUBLE PRECISION(2*N)
C The coefficients of the N givens rotations
C used in the QR factorization of A.
C
C B DOUBLE PRECISION(N+1)
C The right hand side vector.
C
C
C On return
C
C B The solution vector X.
C
C
C Modification of LINPACK.
C Peter Brown, Lawrence Livermore Natl. Lab.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C DAXPY
C
C***END PROLOGUE DHELS
C
INTEGER LDA, N
DOUBLE PRECISION A(LDA,*), B(*), Q(*)
INTEGER IQ, K, KB, KP1
DOUBLE PRECISION C, S, T, T1, T2
C
C Minimize (B-A*X,B-A*X).
C First form Q*B.
C
DO 20 K = 1, N
KP1 = K + 1
IQ = 2*(K-1) + 1
C = Q(IQ)
S = Q(IQ+1)
T1 = B(K)
T2 = B(KP1)
B(K) = C*T1 - S*T2
B(KP1) = S*T1 + C*T2
20 CONTINUE
C
C Now solve R*X = Q*B.
C
DO 40 KB = 1, N
K = N + 1 - KB
B(K) = B(K)/A(K,K)
T = -B(K)
CALL DAXPY (K-1, T, A(1,K), 1, B(1), 1)
40 CONTINUE
RETURN
C
C------END OF SUBROUTINE DHELS------------------------------------------
END
Reference in a new issue View git blame Copy permalink