/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:04 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "sdaslx.h"
#include <float.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/* PARAMETER translations */
#define ICNSTR 10
#define IDB 6
#define IFIRST 16
#define IH0 8
#define IMAT 5
#define IMAXH 7
#define INITYP 11
#define IORD 9
#define IOUT 3
#define ISTOP 4
#define ITOL 2
#define IXSTEP 12
#define LALPHA 11
#define LBETA 17
#define LCJ 1
#define LCNSTR 10
#define LDELT 4
#define LDELTA 46
#define LE 23
#define LERRAB 32
#define LERRAC 62
#define LERRAD 62
#define LERRAE 84
#define LERRAF 84
#define LERRAG 86
#define LERRAH 88
#define LERRAI 88
#define LERRAJ 88
#define LERRAK 88
#define LERRAL 88
#define LERRAM 88
#define LERRAN 88
#define LERRAO 88
#define LERRAP 88
#define LERRAQ 88
#define LERRAR 88
#define LERRAS 88
#define LERRAT 88
#define LERRAU 88
#define LERRAV 88
#define LERRAW 88
#define LERRAX 88
#define LERRAY 88
#define LERRAZ 88
#define LERRBA 88
#define LERRBB 88
#define LERRBC 88
#define LERRBD 88
#define LERRBE 88
#define LERRBF 99
#define LGAMMA 23
#define LH 4
#define LHMAX 3
#define LHMIN 10
#define LK 7
#define LKOLD 8
#define LMAT 9
#define LML 1
#define LMU 2
#define LMXORD 6
#define LNJE 15
#define LNPD 18
#define LNRE 14
#define LNST 13
#define LNSTL 12
#define LPHI 25
#define LPSI 29
#define LROUND 9
#define LSIGMA 35
#define LTN 5
#define LTSTOP 2
#define LTXTAA 1
#define LTXTAB 9
#define LTXTAC 100
#define LTXTAD 259
#define LTXTAE 341
#define LTXTAF 444
#define LTXTAG 515
#define LTXTAH 624
#define LTXTAI 744
#define LTXTAJ 813
#define LTXTAK 941
#define LTXTAL 1095
#define LTXTAM 1146
#define LTXTAN 1196
#define LTXTAO 1270
#define LTXTAP 1344
#define LTXTAQ 1417
#define LTXTAR 1437
#define LTXTAS 1457
#define LTXTAT 1499
#define LTXTAU 1568
#define LTXTAV 1618
#define LTXTAW 1674
#define LTXTAX 1712
#define LTXTAY 1767
#define LTXTAZ 1832
#define LTXTBA 1930
#define LTXTBB 1978
#define LTXTBC 2038
#define LTXTBD 2146
#define LTXTBE 2201
#define LTXTBF 2378
#define LWM 5
#define LWT 24
#define MECONT 50
#define MEEMES 52
#define MEFVEC 61
#define MERET 51
#define METDIG 22
#define METEXT 53
#define MXSTEP 22
#define NTEMP 26
#define REVLOC 21
/* end of PARAMETER translations */
void /*FUNCTION*/ sdaslx(
void (*sdasf)(),
long neq,
float *t,
float y[],
float yprime[],
float tout,
long info[],
float rtol[],
float atol[],
long *idid,
float rwork[],
long lrw,
long iwork[],
long liw)
{
static LOGICAL32 done, nomat;
static long int _l0, i, idat[14], ires, itemp, ldd, leniw, lenpd,
lenrw, locate, mband, msave, nzflg;
static float atoli, edat[5], h, hmax, ho, r, rh, rtoli, tdist,
tn, tnext, tstop, ypnorm;
static char mtxtaa[11][232]={"SDASLX$BEvaluation terminated integration by returning IRES=-2.$NCurrent T = $F, Stepsize H = $F.$EToo much accuracy requested for machine precision.$NRTOL(:) and ATOL(:) replaced by appropriate values.$NCurrent T = $F. Last and ne",
"xt step BDF order =$I,$I.$EDid not reach TOUT in MXSTEPS steps.$NCurrent T = $F, TOUT = $1$F, MXSTEPS=$2$I.$ECorrector could not converge because evaluation returned IRES = -1.$NCurrent T = $F, Stepsize H = $F.$EAfter starting, an ",
"entry of WT(:) has become <= 0.0.$NCurrent T = $F.$EThe error test failed repeatedly.$NCurrent T = $F, Stepsize H =$ $F. $NLast and next step BDF order =$I,$I.$EThe corrector failed$ to converge repeatedly.$NCurrent T = $F, Stepsiz",
"e H = $F. $NLast and next step BDF order =$I,$I.$EThe iteration matrix is singular.$NCurrent T = $F, Stepsize H = $F.$ERepeated corrector convergence failures with singular matrices$ flagged.$NCurrent T = $F, Last and next step BDF",
" order =$I,$I.$ECould not initially solve G(t,y,y')=0 for y'. Consider using larger tolerances.$NCurrent T = $F,$ Stepsize H = $F. Last and next step BDF order =$I,$I.$ESetting INFO($4$I) = $I is not an allowed option.$ENumber of ",
"equations, NEQ, is <= 0. NEQ = $12$I.$EMAXORD, maximum order of BDF used, must be >= 1 and <= 5. MAXORD = $3$I.$ESize of RWORK(:) is too small.$NNow have LRW = $6$I, but need LRW =$ $I.$ESize of IWORK(:) is too small.$NNow have LI",
"W = $8$I, but need LIW = $I.$ERTOL($4$I) is < 0.$EATOL($4$I) is < 0.$EAll entries of RTOL(:) and ATOL(:) = 0.0$EYour TOUT of $2$F is > than your TSTOP of $1$F whidh is not allowed$EValue of maximum stepsize HMAX = $3$F is <= 0.0.$E",
"The current T = $F with H=$F is preceded by TOUT = $F.$EValue of initial stepsize H0 is = 0.$EWith a starting T = $F, the TOUT = $1$F is too close.$EWith the current T=$F and H=$F, the TSTOP=$2$F is inconsistent.$ELower bandwith is",
" $10$I, upper bandwidth is $I and both should be .ge. 0, and < neq which is $I.$EThe current T = $F and TOUT = $1$F are equal.$ECode does not handle constraints when using a band matrix.$EFlag not set for computing constraint corr",
"ection as it must be with a user defined matrix and constraints.$EAt the initial point constraints were not consistent.$EAt T = $F, with H = $F, the constraints appear inconsistent. The norm of the residual with maximal perturbati",
"on of the solution is $F and we did not reduce the norm below $F.$ELast step terminated with IDID < 0. No appropriate action was taken. IDID = $I. Repeated occurrences of illegal input. Run terminated. Apparent infinite loop.$E "};
static char mtxtab[1][67]={"Adjustments to Y in order to get back on the constraints were:$N$B"};
static char mtxtac[1][39]={"Residuals in the constraints were:$N$B"};
static long mloc[32]={LTXTAA,LTXTAB,LTXTAC,LTXTAD,LTXTAE,LTXTAF,
LTXTAG,LTXTAH,LTXTAI,LTXTAJ,LTXTAK,LTXTAL,LTXTAM,LTXTAN,LTXTAO,
LTXTAP,LTXTAQ,LTXTAR,LTXTAS,LTXTAT,LTXTAU,LTXTAV,LTXTAW,LTXTAX,
LTXTAY,LTXTAZ,LTXTBA,LTXTBB,LTXTBC,LTXTBD,LTXTBE,LTXTBF};
static long mact[5]={MEEMES,0,0,0,MERET};
static long mactv[6]={METDIG,6,METEXT,MEFVEC,0,MECONT};
static long infobd[16-(2)+1]={1,1,1,14,3422444,1,1,5,1000,1,1000000,
0,0,0,0};
static long maperr[31]={LERRAB,LERRAC,LERRAD,LERRAE,LERRAF,LERRAG,
LERRAH,LERRAI,LERRAJ,LERRAK,LERRAL,LERRAM,LERRAN,LERRAO,LERRAP,
LERRAQ,LERRAR,LERRAS,LERRAT,LERRAU,LERRAV,LERRAW,LERRAX,LERRAY,
LERRAZ,LERRBA,LERRBB,LERRBC,LERRBD,LERRBE,LERRBF};
/* OFFSET Vectors w/subscript range: 1 to dimension */
float *const Atol = &atol[0] - 1;
float *const Edat = &edat[0] - 1;
long *const Idat = &idat[0] - 1;
long *const Info = &info[0] - 1;
long *const Iwork = &iwork[0] - 1;
long *const Mact = &mact[0] - 1;
long *const Mactv = &mactv[0] - 1;
long *const Maperr = &maperr[0] - 1;
long *const Mloc = &mloc[0] - 1;
float *const Rtol = &rtol[0] - 1;
float *const Rwork = &rwork[0] - 1;
float *const Y = &y[0] - 1;
float *const Yprime = &yprime[0] - 1;
/* end of OFFSET VECTORS */
/* Copyright (c) 2006, Math a la Carte, Inc.
*>> 2009-10-31 sdaslx Krogh Insured that lenpd was defined.
*>> 2008-11-13 sdaslx Hanson correct out of date comments about starting
*>> 2008-08-27 sdaslx Hanson correct leading dimension for banded proble
*>> 2008-08-26 sdaslx Hanson add argument of leading dimension to sdasf
*>> 2008-04-18 sdaslx Hanson allowed INFO(11) =0,1,2 for new startup
*>> 2006-05-19 sdaslx Hanson flag (-29 -> -31) & new const errs -29, -30
*>> 2006-04-24 sdaslx Krogh Added logical variable, nomat.
*>> 2003-07-15 sdaslx Hanson Install Soderlind stepsize code.
*>> 2002-06-26 sdaslx Krogh Insured iwork(lk) has current value.
*>> 2001-12-29 sdaslx Krogh Many changes in usage and documentation.
*>> 2001-12-12 sdaslx Krogh Changed code for reverse communication
*>> 2001-11-23 sdaslx Krogh Changed many names per library conventions.
*>> 2001-11-04 sdaslx Krogh Fixes for F77 and conversion to single
*>> 2001-11-01 sdaslx Hanson Provide code to Math a la Carte.
*--S replaces "?":?daslx, ?daswt, ?das1, ?dasnm, ?dastp,
*--& ?dasin, ?daswt, ?dasf, ?dasdb, ?mess
****BEGIN PROLOGUE SDASLX
****PURPOSE This code solves a system of differential/algebraic
* equations of the form G(T,Y,YPRIME) = 0.
****LIBRARY SLATEC (SDASLX)
****CATEGORY I1A2
****TYPE DOUBLE PRECISION (SDASSL-S, SDASLX-D)
****KEYWORDS BACKWARD DIFFERENTIATION FORMULAS, DASSL,
* DIFFERENTIAL/ALGEBRAIC, IMPLICIT DIFFERENTIAL SYSTEMS
****AUTHOR Petzold, Linda R., (LLNL)
* Computing and Mathematics Research Division
* Lawrence Livermore National Laboratory
* L - 316, P.O. Box 808,
* Livermore, CA. 94550
* Changes made by R. J. Hanson, Math a la Carte.
****DESCRIPTION
* Subroutine SDASLX uses the backward differentiation formulas of
* orders one through five to solve a system of the above form for Y and
* YPRIME. Values for T and Y at the initial time must be given as
* input. The subroutine solves the system from T to TOUT. It is easy
* easy to continue the solution to get results at additional TOUT.
* This is the interval mode of operation. Intermediate results can
* also be obtained easily by using the intermediate-output capability.
* The code automatically selects the order and step size.
*
* *Usage:
*
* EXTERNAL SDASF
* INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW),
* REAL T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL,
* * RWORK(LRW)
*
* CALL SDASLX (SDASF, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
* IDID, RWORK, LRW, IWORK, LIW)
*
* ************************* Arguments **********************************
*
* SDASF:EXT -- Depending on how INFO(5) is set, this routine may be
* used to compute g, compute the iteration matrix, and/or to factor
* and solve the linear systems as requested. It has the form
* SUBROUTINE SDASF(T,Y,YPRIME,DELTA,PD,LDD,CJ,IRES,RWORK,IWORK)
* The actions to be taken depend on the value of the IRES argument
* to SDASF. IRES also provides a way for the user to signal certain
* things to SDASLX. The values of IRES are
* = 0 Initialize subroutine parameters or any storage requirements.
* If this does not apply, then ignore this value and return
* without any further action. Initial conditions for t, y and y'
* can be assigned to T, Y(:) and YPRIME(:). Do not alter
* DELTA, D or C when IRES = 0. */
/* When reverse communication is used, IRES is in IWORK(3), DELTA
* starts in RWORK(IWORK(4)), the matrix PD starts in
* RWORK(IWORK(5)), and CJ is in RWORK(1).
* In your code you can ignore any values for IRES which can not
* occur for the value of INFO(5) which is set. */
/* = 1 Evaluate the residual of the differential/algebraic system.
* Place values in the array DELTA = g}(t, y, y'). */
/* = 2 This case occurs for |INFO(5)| = 2, 4, 5, 6, 8, and 10. */
/* Evaluate the partial derivatives (p used for partial
* symbol), D(i,j) = (pg_i / py_j) + c (pg_i) / py'_j.
* The value LDD is the row dimension of D(*,*). The
* scalar c is an input argument to SDASF. For the given values of
* T, Y, YPRIME, the routine must evaluate the non-zero partial
* derivatives for each equation, and store these values in the
* matrix D. For dense or banded matrices, stored in the work
* space RWORK(:), the elements of D are set to zero before each
* call to SDASF so only non-zero elements need to be defined. Do
* not alter T, Y, YPRIME, or C. You must dimension D with a first
* first dimension of LDD. The way you must store the elements
* into D depends on the structure of the matrix, indicated by
* INFO(5) as follows. */
/* INFO(5)=2, 8, and 12: Full (dense) matrix. When you evaluate
* the (non-zero) partial derivative of equation i with
* respect to variable j, you must store it in D according
* to D(i,j) = (p g_{i} / py_{j})+c (pg_{i}/py'_{j}).
* INFO(5)=4, 10, and 14: Banded Jacobian with ML lower and MU
* upper diagonal bands (refer to INFO(6) description of ML
* and MU). Give D a first dimension of 2*ML+MU+1. When you
* evaluate the (non-zero) partial derivative of equation i
* with respect to variable j, you must store it in D */
/* according to irow=i-j+ ML+MU+1, D(irow,j) = (pg_i /py_j) +
* c (pg_i / py'_j).
* INFO(5) = 5 and 6: The array D is an unused dummy argument.
* Save the partials (pg_i / py_j) + c (p g_{i} / py'_j) in
* the subroutine SDASF. This case requires that you factor
* the matrix (if you ever do) at this time as the value
* IRES=3 is not provided in this case. As for that case you
* should set IRES = 0 to flag the fact that there were no
* problems in obtaining the factorization. */
/* = 3 This case occurs for |INFO(5)| > 6. Factor the matrix of
* partials. Prepare to solve systems involving the matrix D =
* pg_i / py_j) + c (pg_i / py'_j). The solution method can be a
* convenient one of the user's choice. If the matrix is non-
* singular it is important to return IRES = 0 as a signal.
* Otherwise return IRES = 3, if the system is nearly singular. */
/* = 4 Solve a linear algebraic system using the matrix of partials D,
* i.e. solve Dw = r for w. The vector r is input in array
* DELTA(:). The solution w overwrites DELTA(:). If for any
* reason the system cannot be solved, return w = r as the
* approximate solution. This may cause the integrator to take
* corrective action such as reducing the step-size or the order of
* the formulas used. This situation may occur when iteratively
* solving a linear system, but requiring an excessive number of
* iterations. */
/* = 5 Compute the residual and partials for projecting the solution
* Y(:) onto the constraints after a step has been computed and the
* corrector equation has converged. If you are handling the
* linear algebra (|INFO(5)| .ge. 5) then you should also compute
* the projection and store it in DELTA(:). The SDASLX code
* applies the projection, Y(:) = Y(:) - DELTA(:). (Note that the
* flag is given to SDASF if INFO(5) .ge. 0, and else is provided
* using reverse commiunication. If for example you are using
* forward communication for derivatives and doing you own linear
* algebra using reverse communication, you will need to either
* deal with the linear algebra in SDASF, or call a routine from
* there that will do the job.) */
/* Ordinarily subroutine SDASF should not change the value of IRES.
* The following values can be set for special cases.
* 0 This must be set if you are factoring the iteration matrix, to
* let SDASLX know that your matrix is not singular.
* -1 Some kind of difficulty has been encountered. This causes
* SDASLX to reduce the stepsize or order which may cure the
* problem.
* -2 Return immediately to the main program, for one reason or the
* other it is time to quit.
* <-2 This has the same effect as setting INFO(6) to the negative of
* this number. It provides a way of turning on debug print at
* any time. */
/* NEQ:IN This is the number of equations to be solved. NEQ .ge. 1.
*
* T:INOUT This is the current value of the independent variable.
* Initialize to the value at the initial point.
*
* Y(.ge. NEQ):INOUT This array contains the solution components at T.
* Initialize to the value at the initial point.
*
* YPRIME(.ge. NEQ):INOUT This array contains the derivatives of the
* solution components at T. You can either provide intial values,
* or use INFO(11) = 1 to have them computed.
*
* If =1, use the method given in the origial DDASL code of Petzold.
* If =2, use the Newton method of Krogh and Hanson,
* Solving Constrained Differential-Algebraic Systems Using
* Projections, 2008. This is recommended but requires more storage
* than the option with vaue =1.
*
* TOUT:IN This is the first point at which a solution is desired. The
* value of TOUT relative to T defines the direction of the
* integration, so TOUT = T is not allowed at the intial point.
*
* INFO(.ge. 16):IN The basic task of the code is to solve the system
* from T to TOUT and return an answer at TOUT. INFO is an integer
* array which is used to communicate exactly how you want this task
* to be carried out. You must provide inputs for the first 16
* entries. The simplest use of the code is obtained by setting all
* entries to 0, although you probably want to check INFO(5) at
* least
*
* INFO(1) - This parameter enables the code to initialize itself. You
* must set it to indicate the start of every new problem.
*
* **** Is this the first call for this problem ...
* Yes - Set INFO(1) = 0
* No - Not applicable here.
* See below for continuation calls. ****
*
* INFO(2) - How much accuracy you want of your solution is specified
* by the error tolerances RTOL and ATOL. The simplest use is to
* take them both to be scalars. To obtain more flexibility, they
* can both be vectors. The code must be told your choice.
*
* **** Are both error tolerances RTOL, ATOL scalars ...
* Yes - Set INFO(2) = 0 */
/* and input scalars for both RTOL and ATOL
* No - Set INFO(2) = 1
* and input arrays for both RTOL and ATOL ****
*
* INFO(3) - The code integrates from T in the direction of TOUT by
* steps. If you wish, it will return the computed solution and
* derivative at the next intermediate step (the intermediate-output
* mode) or TOUT, whichever comes first. This is a good way to
* proceed if you want to see the behavior of the solution. If you
* must have solutions at a great many specific TOUT points, this
* code will compute them efficiently. You can change this mode
* at any time.
*
* **** Do you want the solution only at
* TOUT (and not at the next intermediate step) ...
* Yes - Set INFO(3) = 0
* No - Set INFO(3) = 1 ****
*
* INFO(4) - To handle solutions at a great many specific values TOUT
* efficiently, this code may integrate past TOUT and interpolate to
* obtain the result at TOUT. Sometimes it is not possible to
* integrate beyond some point TSTOP because the equation changes
* there or it is not defined past TSTOP. Then you must tell the code
* not to go past. The code will never give values beyond TSTOP, but
* you can change this mode at any time. However when you set
* INFO(4) to 1, you must always set RWORK(1) to a value past the
* current time.
*
* **** Can the integration be carried out without any
* restrictions on the independent variable T ...
* Yes - Set INFO(4)=0
* No - Set INFO(4)=1
* and define the stopping point TSTOP by
* setting RWORK(1)=TSTOP ****
*
* INFO(5) - To solve differential/algebraic problems it is necessary
* to use a matrix of partial derivatives for the system of
* differential equations. This flag must be set to define your
* mtarix.
* = 1 (or 0) Matrix is full and dense and numerical derivatives are
* generated internally. Although it is less trouble for you to have
* the code compute derivatives by numerical differencing, the
* solution will be more reliable if you provide the derivatives
* using formulas.
* = 2 Matrix is full and dense and the user is computing derivatives.
* = 3 As for 1, but the matrix, (p stands for the partial operator)
* D = pg/py+ c (pg}/py', is banded. Here c is a factor determined
* by SDASLX. The differential equation is said to have
* half-bandwidths ML (lower) and MU (upper) if D involves only some
* of the unknowns y_{J} with i - ML .le. J .le. i + MU for all
* i=1,...,NEQ. Thus, ML and MU are the widths of the lower and
* upper parts of the band, with the main diagonal being excluded.
* If the matrix is stored in the RWORK(:) array and you do not
* indicate that the equation has a banded matrix of partial
* derivatives, the code works with a full array of NEQ**2 elements,
* stored in the conventional Fortran array style. Computations with
* banded matrices typically require less time and storage than with
* full matrices, if 2 * ML + MU < NEQ. If you tell the code that
* the matrix of partial derivatives has a banded structure and you
* want to provide subroutine SDASF to compute the partial
* derivatives, then you must store the elements of the matrix in the
* Linpack band-matrix specification, indicated in the description of
* SDASF.
* Provide the lower ML and upper MU bandwidths by setting
* IWORK(1) = ML and IWORK(2) = MU.
* = 4 As for 3, but derivatives are to be computed by the user.
* = 5 The user is going to take care of all matters connected with the
* factorization and solving of the systems.
* = 6 As for 5, except the linear algebra is to be accomplished using
* reverse communication.
* = 7 - 10 As for 1 - 4, except user is doing linear algebra in SDASF.
* = 11 - 14 As for 7 - 10, except using reverse communication.
* = k<0 Same as -k, except the evaluation of g and all partials is done
* using reverse communciation rather than by calling SDASF.
*
* INFO(6) - Set the level of debugging print. A value of 0 gives no
* print. Otherwise this is a 7 digit number d_6d_5d_4d_3d_2d_1d_0
* defining what to print as follows.
* d_0 Print on entry to SDASLX
* = 0 No print
* = 1 IDID, INFO(1), NEQ, T, TOUT
* = 2 y, y'
* = 3 Tolerances
* = 4 INFO (all of it)
* d_1 Print on exit from SDASLX. Print is as for d_0
* d_2 After a call to SDASF (or return from reverse communication
* that would ordinarily call SDASF.
* = 0 No print
* = 1 Print whatever was just computed, except no matrices.
* = 2 Print whatever was just computed, including matrices.
* d_3 As for the case above, except print is for what is in the
* locations about to be stored to.
* d_4 Internal print inside subroutine SDASTP.
* = 0 No print
* = 1 y, y' and corrections.
* = 2 tolerances
* = 3 difference tables
* = 4 integration coefficients
* d_5 Determines how much of WORK and IWORK are printed, when
* there is other print. */
/* = 0 No print
* = 1 Always print IWORK(1:16)
* = 2 Always print WORK(1:9)
* = 3 Always print both of the above.
* d_6 For turning off, or limiting the amount of print.
* = 0 No effect
* = 1 No effect, but gives a way to specify a value of 0, 1 or 2
* when passing a negative value of IRES after starting.
* > 1 Print data for just this many of the first variables, and
* just this many of the first rows in matrices when
* variables or matrices are printed.
*
* INFO(7) - You can specify a maximum (absolute value of) stepsize, so
* that the code will avoid passing over very large regions.
*
* **** Do you want the code to decide on its own maximum stepsize?
* Yes - Set INFO(7)=0
* No - Set INFO(7)=1
* and define HMAX by setting
* RWORK(2)=HMAX ****
*
* INFO(8) - Differential/algebraic problems may occasionally suffer
* from severe scaling difficulties on the first step. If you know a
* great deal about the scaling of your problem, you can help to
* alleviate this problem by specifying an initial stepsize HO.
*
* **** Do you want the code to define its own initial stepsize?
* Yes - Set INFO(8)=0
* No - Set INFO(8)=1
* and define HO by setting
* RWORK(3)=HO ****
*
* INFO(9) - If storage is a severe problem, you can save some
* locations by restricting the maximum order MAXORD. the default
* value is 5. for each order decrease below 5, the code requires NEQ
* fewer locations, however it is likely to be slower. In any case,
* you must have 1 .LE. MAXORD .LE. 5
*
* **** Do you want the maximum order to default to 5?
* Yes - Set INFO(9)=0
* No - Set INFO(9)=MAXORD .le. 5
*
* INFO(10) - If you know that the solutions to your equations satisfy
* some constraints you may want to use this option. It is expected
* to be especially useful when one differentiates to reduce the
* index of a problem. This option must not be used when INFO(5) =
* 3 or 4.
*
* **** Do you want the code to solve the problem without
* invoking constraints?
* Yes - Set INFO(10)=0
* No - Set INFO(10)=k, where k is >0 if one wants
* to solve for k constraints using forward
* communication, and k < 0 to solve for -k
* constraints using reverse communication.
* Artificially require |k| .le. 1000.
*
* INFO(11) - SDASLX normally requires the initial T, Y, and YPRIME to
* be consistent. That is, you must have G(T,Y,YPRIME) = 0 at the
* initial time. If you do not know the initial derivative precisely,
* you can let SDASLX try to compute it.
*
* **** Are the initial T, Y, YPRIME consistent?
* Yes - Set INFO(11) = 0
* No - Set INFO(11) = 1, and set YPRIME to an initial
* approximation to YPRIME. (If you have no idea what
* YPRIME should be, set it to zero. Note that the
* initial Y should be such that there must exist a
* YPRIME so that G(T,Y,YPRIME) = 0.)
* No - Set INFO(11) = 2, and set YPRIME to an initial
* approximation to YPRIME. This option uses the
* Krogh/Hanson starting algorithm found in the report
* Solving Constrained Differential-Algebraic Systems
* Using Projections, 2008. This is the preferred
* option if the code solves for YPRIME.
*
* INFO(12) - SDASLX} normally allows up to 500 internal steps between
* each output points.
* Set INFO(12)=0 for the code to use up to 500 internal steps
* between output points.
* Set INFO(12)=k} and the code will use up to k internal steps
* between output points.
* INFO(13) - SDASLX} normally uses the smoothed step control algorithm
* developed by Gustaf Soderlind.
* Set INFO(13)=0 for the code to use the step control method
* of Soderlind.
* Set INFO(13)=1 and the code will use the original step
* control logic of DASSL.
*
* INFO(14:16) -- Not used currently by the code, but must be set to 0. */
/* RTOL,ATOL:INOUT These quantities represent relative and absolute
* error tolerances which you provide to indicate how accurately you
* wish the solution to be computed. You may choose them to be both
* scalars (INFO(2)=0) or else both vectors (INFO(2)=1). Caution:
* In Fortran 77, a scalar is not the same as an array of length 1.
* Some compilers may require to make the scalars into arrays of
* length 1. All components must be non-negative.
* The tolerances are used by the code in a local error test at each */
/* step which requires roughly that
* ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL
* for each vector component. (More specifically, a root-mean-square
* norm is used to measure the size of vectors, and the error test
* uses the magnitude of the solution at the beginning of the step.)
*
* The true (global) error is the difference between the true
* solution of the initial value problem and the computed
* approximation. Practically all present day codes, including this
* one, control the local error at each step and do not even attempt
* to control the global error directly. Usually, but not always,
* the true accuracy of the computed Y is comparable to the error
* tolerances. This code will usually, but not always, deliver a more
* accurate solution if you reduce the tolerances and integrate
* again. By comparing two such solutions you can get a fairly
* reliable idea of the true error in the solution at the bigger
* tolerances.
*
* Setting ATOL=0. results in a pure relative error test on that
* component. Setting RTOL=0. results in a pure absolute error test
* on that component. A mixed test with non-zero RTOL and ATOL
* corresponds roughly to a relative error test when the solution
* component is much bigger than ATOL and to an absolute error test
* when the solution component is smaller than the threshold ATOL.
*
* The code will not attempt to compute a solution at an accuracy
* unreasonable for the machine being used. It will advise you if
* you ask for too much accuracy and inform you as to the maximum
* accuracy it believes possible.
*
* IDID:OUT This scalar quantity is an indicator reporting what the
* code did. You must monitor this integer variable to decide what
* action to take next. Values returned are as follows.
* Task Completed or Ongoing
* 1 A step was successfully taken in the intermediate-output
* mode. The code has not yet reached TOUT.
* 2 The integration to TSTOP was successfully completed to
* T = TSTOP by stepping exactly to TSTOP.
* 3 The integration to TOUT was successfully completed T = TOUT by
* stepping past TOUT. Y(:)} is obtained by interpolation.
* YPRIME(:) is obtained by interpolation.
* 4 The integration has paused for reverse communication. Respond
* based on the values of IWORK(3).
* Task Interupted
* -1 IRES set to -2 by the user.
* -2 Accuracy requested exceeds machine precision. RTOL and ATOL
* have been increased.
* -3 There have been too many steps between output points.
* Quit or Restart Integration
* -4 No convergence due to IRES being set to -1.
* -5 A weight for computing error norms is .le. 0.
* -6 The error test has failed repeatedly.
* -7 Repeated failure of the corrector to converge.
* -8 The iteration matrix is singular.
* -9 Repeated corrector convergence failures, with singular
* matrices flagged.
* -10 Could not solve for the initial y'.
* Invalid input
* -11 An INFO entry has a value not allowed for that option.
* -12 The number of equations was set .le. 0.
* -13 The maximum order does not have a value from 1 to 5.
* -14 The size of RWORK is too small.
* -15 The size of IWORK is too small.
* -16 An entry of RTOL is < 0.
* -17 An entry of ATOL is < 0.
* -18 All entries of RTOL and ATOL are 0.
* -19 The value of TOUT is > TSTOP.
* -20 The maximum stepsize is set .le. 0.
* -21 The current TOUT is behind T.
* -22 The initial stepsize has been set to 0.
* -23 TOUT is too close to the starting T.
* -24 TSTOP is not consistent with the current T.
* -25 An illegal bandwidth.
* -26 The current T and TOUT are equal.
* -27 Constraints used with band matrices.
* -28 When solving constraints for a user defined matrix, IRES was
* not changed from 5.
* -29 Inconsistency in the constraints when starting.
* -30 During later integration steps, an inconsistency and rank
* deficiency were noted. Print amount constraint residual norm
* and linear system residual norm.
* -31 No appropriate action was taken after IDID was set < 0. */
/* RWORK:WORK A real work array of length LRW which provides the
* code with needed storage space. The data in RWORK is necessary
* for subsequent calls. You may find the following locations of
* interest.
* RWORK(3) The step size H to be attempted on the next step.
* RWORK(4) The current value of the independent variable, i.e.,
* the farthest point integration has reached. This will
* be different from T only when interpolation has been
* performed (IDID=3).
* RWORK(7) The stepsize used on the last successful step.
*
* LRW:IN The length of RWORK. Set it to the declared length of the
* RWORK array. You must have LRW .GE. 45+(MAXORD+4)*NEQ
* + NEQ**2 for the full (dense) JACOBIAN case
* (INFO(5) =0,1,2,7,8,11,12), or
* + (2*ML+MU+1)*NEQ for the banded user-defined JACOBIAN case */
/* (INFO(5)=3, 9, 13), or
* + (2*ML+MU+1)*NEQ+2*(NEQ/(ML+MU+1)+1) for the banded finite-
* difference-generated JACOBIAN (INFO(5)=4, 10, 14),
* + 0 For the user defined matrix case.
*
* IWORK:WORK An integer work array of length LIW which provides the
* code with needed storage space. As for RWORK, you may find the
* following locations of interest.
* IWORK(4) The order of the method to be tried on the next step.
* IWORK(5) The order of the method used on the last step.
* IWORK(8) The number of steps taken so far.
* IWORK(9) The number of calls to SDASF so far.
* IWORK(10) The number of evaluations of the matrix of partial
* derivatives needed so far.
* IWORK(11) The total number of error test failures so far.
* IWORK(12) Which contains the total number of convergence test
* failures so far. (Includes singular iteration matrix
* failures.)
*c
* LIW:IN The length of IWORK. Set it to the declared length of the
* IWORK array. You must have LIW .GE. 30+NEQ
*
* OPTIONALLY REPLACEABLE NORM ROUTINE:
*
* SDASLX uses a weighted norm SDASNM to measure the size of vectors
* such as the estimated error in each step. A FUNCTION subprogram
* REAL FUNCTION SDASNM(NEQ,V,WT,RWORK,IWORK)
* DIMENSION V(NEQ),WT(NEQ) */
/* is used to define this norm. Here, V is the vector whose norm is
* to be computed, and WT is a vector of weights. A SDASNM routine
* has been included with SDASLX which computes the weighted
* root-mean-square norm given by SDASNM = SQRT((1/NEQ) *
* SUM(V(I)/WT(I))**2) this norm is suitable for most problems. In
* some special cases, it may be more convenient and/or efficient to
* define your own norm by writing a function subprogram to be called
* instead of SDASNM. This should, however, be attempted only after
* careful consideration.
*
* ********************* Local Variables *******************************
*
* NOTE -- Parameter names are given followed by (ii name) where ii
* is the integer value of the parameter, and name is the array the
* parameter is used to reference.
*
* icnstr (10 info) Number of constraints.
* idb (6 info) Digits define rules for debugging. From right to left
* the digits define output for: entry to sdaslx, exit from sdaslx,
* after actions in sdasf, before actions id sdasf, actions inside
* sdastp, and what to do about printing iwork and rwork. The bigger
* the digit, the more the print; see above for details.
* ih0 (8 info) If set, the initial stepsize is stored in rwork(lh=3).
* imat (5 info) Defines how the Jacobian is computed and processed.
* If < 0, then computation of g and of all partials is done using
* reverse communication. Absolute value is used as follows:
* 1.Full matrix, derivatives computed using finite differences.
* 2 Full matrix, user computes derivatives.
* 3 Band matrix, ML=iwork(lml=1), MU=iwork(lmu=2) define the left
* and right bandwidths. Derivative approximated with differences.
* 4 As for 3, but with the derivatives computed with user supplied
* derivatives.
* 5 User defined matrix. User does all of the work!
* 6 As for 5, but the work is done using reverse communication.
* 7-10 As for 1 - 4, but user does linear algebra themselves inside
* sdasf.
* 11-14 As for 7-10, but reverse communication is used.
* imaxh (7 info) If Set, the maximum stepsize has been set by the
* user in rwork(lhmax=2)
* inityp (11 info) If set computations must be done to get
* consistent initial value for y'.
* iord (9 info) If set, gives the maximum integration order to use.
* iout (3 info) If set, gives output at the end of every step.
* istop (4 info) If set, then the user has set rwork(tstop=1) to a
* value that the integrator must not go past.
* itol (2 info) If set, then atol() and rtol() are arrays, else they
* are scalars.
* ixstep (12 info) If set, then iwork(mxstep=21) has been set to
* info(istop=4) giving the maximum number of steps between outputs.
* lalpha (11 rwork) Start of the alphas used in computing integ. coeff.
* lbeta (17 rwork) Start of the betas used in computing integ. coeff.
* lcj (1 rwork) The constant which multiplies the partials with
* respect to y' in the iteration matrix.
* lcjold (6 rwork) The previous value of rwork(lcj=5).
* lcnstr (10 iwork) The number of constraints.
* lctf (17 iwork) Counts the total number of various convergence
* failures.
* ldelt (4, iwwork) Always equal to ldelta below.
* ldelta (46 rwork) Place where residuals and corrections to y are
* stored.
* le (23 iwork) rwork(iwork(le=22) is where error estimates are stored.
* letf (16 iwork) Contains the number of error test failures so far.
* lgamma (23 rwork) Start of the gammas used in computing integ. coeff.
* lh (4 rwork) Stepsize to be used on the next step.
* lhmax (3 rwork) The largest stepsize permitted.
* lhmin (10 rwork) The minimum stepsize permitted.
* lhold (7 rwork) The previous stepsize.
* lipvt (31 iwork) Used for permutations when solving linear systems.
* lires (3 iwork) Value of IRES for reverse communication.
* ljcalc (19 iwork) This is set to -1 when a new Jacobian is needed,
* set to 0 when it is computed, and set to 1 when it has been used. */
/* lk (7 iwork) Contains the current integration order.
* lkold (8 iwork) Contains the previous integration order.
* lmat (9 iwork) Type of matrix, from info(imat).
* lml (1 iwork) Contains the left bandwith (if used).
* lmu (2 iwork) Contains the right bandwidth (if used).
* lmxord (6 iwork) Contains the maximum integration order allowed.
* lnjac (8 rwork) A large number here means we need a new Jacobian.
* Small number means that convergence is good.
* lnje (15 iwork) Counts the number of Jacobian evaluations so far.
* lnpd (18 iwork) This gives the space required to hold partial
* derivatives. Includes partials for constraints.
* lnre (14 iwork) Counts the number of evaluations of g.
* lns (11 iwork) Counts the number of steps without a stepsize change.
* lnst (13 iwork) Counts the number of steps taken so far.
* lnstl (12 iwork) Number of steps when last returned to user (not
* counting computations done in sdasf when not using reverse comm.)
* lphase (20 iwork) This is set to 0 initially allowing rapid order
* and stepsize increases. It is set to 1 if either the maximum
* order is reached, or the order is decreased or there has been some
* kind of failure.
* lphi (25 iwork) Location of difference tables in rwork.
* lpsi (29 rwork) Start of the psi's used in computing integ. coeff.
* lround (9 rwork) The machine epsilon.
* lsigma (35 rwork) Start of the sigmas used in error estimation.
* ltn (5 rwork) The point up to which we have integrated. Local
* variable tn in sdaslx contains this value.
* ltstop (2 rwork) Value that the integrator is not to go past.
* lwm (5 iwork) rwork(iwork(lwm=25)) is the start of matrix data.
* lwt (24 iwork) rwork(iwork(lwt=23)) is the start of weights used
* in computing error norms.
* mxstep (22 iwork) Contains the maximum number of steps to take
* between output points.
* nomat Logical variable set = .true. if matrix storage not used here.
* ntemp (26 iwork) Location relative to work(iwork(lwm=25)) of
* extra matrix data required to compute finite difference Jacobians.
* revloc (21 iwork) This is nonzero if reverse communication is in
* the process of being used. The low order 3 bits define the
* location to go to in the current routine, the next three for the
* routine that this one calls, and the next 3 bits for the routine
* called by that one, etc.
*
* ---------------------------------------------------------------------
*
****REFERENCES A DESCRIPTION OF DASSL: A DIFFERENTIAL/ALGEBRAIC
* SYSTEM SOLVER, L. R. PETZOLD, SAND82-8637,
* SANDIA NATIONAL LABORATORIES, SEPTEMBER 1982.
****ROUTINES CALLED R1MACH, SDAS1, SDASNM, SDASTP, SDASIN, sDASWT,
* XERMSG
****REVISION HISTORY (YYMMDD)
* 830315 DATE WRITTEN
* 880387 Code changes made. All common statements have been
* replaced by a DATA statement, which defines pointers into
* RWORK, and PARAMETER statements which define pointers
* into IWORK. As well the documentation has gone through
* grammatical changes.
* 881005 The prologue has been changed to mixed case.
* The subordinate routines had revision dates changed to
* this date, although the documentation for these routines
* is all upper case. No code changes.
* 890511 Code changes made. The DATA statement in the declaration
* section of SDASLX was replaced with a PARAMETER
* statement. Also the statement S = 100.E0 was removed
* from the top of the Newton iteration in SDASTP.
* The subordinate routines had revision dates changed to
* this date.
* 890517 The revision date syntax was replaced with the revision
* history syntax. Also the "DECK" comment was added to
* the top of all subroutines. These changes are consistent
* with new SLATEC guidelines.
* The subordinate routines had revision dates changed to
* this date. No code changes.
* 891013 Code changes made.
* Removed all occurrences of FLOAT or DBLE. All operations
* are now performed with "mixed-mode" arithmetic.
* Also, specific function names were replaced with generic
* function names to be consistent with new SLATEC guidelines.
* In particular:
* Replaced DSQRT with SQRT everywhere.
* Replaced DABS with ABS everywhere.
* Replaced DMIN1 with MIN everywhere.
* Replaced MIN0 with MIN everywhere.
* Replaced DMAX1 with MAX everywhere.
* Replaced MAX0 with MAX everywhere.
* Replaced DSIGN with SIGN everywhere.
* Also replaced REVISION DATE with REVISION HISTORY in all
* subordinate routines.
* 901004 Miscellaneous changes to prologue to complete conversion
* to SLATEC 4.0 format. No code changes. (F.N.Fritsch)
* 901009 Corrected GAMS classification code and converted subsidiary
* routines to 4.0 format. No code changes. (F.N.Fritsch)
* 901010 Converted XERRWV calls to XERMSG calls. (R.Clemens, AFWL)
* 901019 Code changes made.
* Merged SLATEC 4.0 changes with previous changes made
* by C. Ulrich. Below is a history of the changes made by
* C. Ulrich. (Changes in subsidiary routines are implied
* by this history)
* 891228 Bug was found and repaired inside the SDASLX
* and SDAS1 routines. SDAS1 was incorrectly
* returning the initial T with Y and YPRIME
* computed at T+H. The routine now returns T+H */
/* rather than the initial T.
* Cosmetic changes made to SDASTP.
* 900904 Three modifications were made to fix a bug (inside
* SDASLX) re interpolation for continuation calls and
* cases where TN is very close to TSTOP:
*
* 1) In testing for whether H is too large, just
* compare H to (TSTOP - TN), rather than
* (TSTOP - TN) * (1-4*UROUND), and set H to
* TSTOP - TN. This will force SDASTP to step
* exactly to TSTOP under certain situations
* (i.e. when H returned from SDASTP would otherwise
* take TN beyond TSTOP).
*
* 2) Inside the SDASTP loop, interpolate exactly to
* TSTOP if TN is very close to TSTOP (rather than
* interpolating to within roundoff of TSTOP).
*
* 3) Modified IDID description for IDID = 2 to say
* that the solution is returned by stepping exactly
* to TSTOP, rather than TOUT. (In some cases the
* solution is actually obtained by extrapolating
* over a distance near unit roundoff to TSTOP,
* but this small distance is deemed acceptable in
* these circumstances.)
* 901026 Added explicit declarations for all variables and minor
* cosmetic changes to prologue, removed unreferenced labels,
* and improved XERMSG calls. (FNF)
* 901030 Added ERROR MESSAGES section and reworked other sections to
* be of more uniform format. (FNF)
* 910624 Fixed minor bug related to HMAX (six lines after label
* 525). (LRP)
* 981110 Major revision on several codes in the package. Now
* use one external user routine. Also support reverse
* communication. Install JPL error processor in place of
* XERMSG. */
/****END PROLOGUE SDASLX
*
***End
*
* Declare arguments.
* */
/* Declare externals.
* */
/* Declare local variables.
* */
/* The following locations are reserved for future expansion:
* iwork(26:30), work(41:45),; and info(13:16)
*
* POINTERS INTO IWORK */
/* POINTERS INTO RWORK */
/* POINTERS INTO INFO */
/* ********* Error message text ***************
*[Last 2 letters of Param. name] [Text generating message.]
*AA SDASLX$B
*AB Evaluation terminated integration by returning IRES=-2.$N
* Current T = $F, Stepsize H = $F.$E
*AC Too much accuracy requested for machine precision.$N
* RTOL(:) and ATOL(:) replaced by appropriate values.$N
* Current T = $F. Last and next step BDF order =$I,$I.$E
*AD Did not reach TOUT in MXSTEPS steps.$N
* Current T = $F, TOUT = $1$F, MXSTEPS=$2$I.$E
*AE Corrector could not converge because evaluation returned $C
* IRES = -1.$N
* Current T = $F, Stepsize H = $F.$E
*AF After starting, an entry of WT(:) has become <= 0.0.$N
* Current T = $F.$E
*AG The error test failed repeatedly.$N
* Current T = $F, Stepsize H = $F. $N
* Last and next step BDF order =$I,$I.$E
*AH The corrector failed to converge repeatedly.$N
* Current T = $F, Stepsize H = $F. $N
* Last and next step BDF order =$I,$I.$E
*AI The iteration matrix is singular.$N
* Current T = $F, Stepsize H = $F.$E
*AJ Repeated corrector convergence failures with singular matrices $C
* flagged.$N
* Current T = $F, $C
* Last and next step BDF order =$I,$I.$E
*AK Could not initially solve G(t,y,y')=0 for y'. Consider using $C
* larger tolerances.$N
* Current T = $F, Stepsize H = $F. Last and next $C
* step BDF order =$I,$I.$E
*AL Setting INFO($4$I) = $I is not an allowed option.$E
*AM Number of equations, NEQ, is <= 0. NEQ = $12$I.$E
*AN MAXORD, maximum order of BDF used, must be >= 1 and <= 5. $C
* MAXORD = $3$I.$E
*AO Size of RWORK(:) is too small.$N
* Now have LRW = $6$I, but need LRW = $I.$E
*AP Size of IWORK(:) is too small.$N
* Now have LIW = $8$I, but need LIW = $I.$E
*AQ RTOL($4$I) is < 0.$E
*AR ATOL($4$I) is < 0.$E
*AS All entries of RTOL(:) and ATOL(:) = 0.0$E
*AT Your TOUT of $2$F is > than your TSTOP of $1$F whidh is $C
* not allowed$E
*AU Value of maximum stepsize HMAX = $3$F is <= 0.0.$E
*AV The current T = $F with H=$F is preceded by TOUT = $F.$E
*AW Value of initial stepsize H0 is = 0.$E
*AX With a starting T = $F, the TOUT = $1$F is too close.$E
*AY With the current T=$F and H=$F, the TSTOP=$2$F is inconsistent.$E
*AZ Lower bandwith is $10$I, upper bandwidth is $I and both should $C
* be .ge. 0, and < neq which is $I.$E
*BA The current T = $F and TOUT = $1$F are equal.$E
*BB Code does not handle constraints when using a band matrix.$E
*BC Flag not set for computing constraint correction as it must be $C
* with a user defined matrix and constraints.$E
*BD At the initial point constraints were not consistent.$E
*BE At T = $F, with H = $F, the constraints appear inconsistent. $C
* The norm of the residual with maximal perturbation of the $C
* solution is $F and we did not reduce the norm below $F.$E
*BF Last step terminated with IDID < 0. No appropriate action was $C
* taken. IDID = $I. Repeated occurrences of illegal input. $C
* Run terminated. Apparent infinite loop.$E
* $ D
*BG Adjustments to Y in order to get back on the constraints were:$N$B
* $
*BH Residuals in the constraints were:$N$B */
/*++ Save data by elements if ~.C. */
/* End of storing data by elements. */
/* End of data generated by pmess. */
/* Staging arrays for error printing of data: */
/* 2-4 5 6 7 8 9 10 11 12 13-16 */
/* The following parameters define the severity for stopping and
* printing of the associated error message . */
/****FIRST EXECUTABLE STATEMENT SDASLX */
if (Info[1] == 0)
{
/* Initialization
* ----------------------------------------------------------------------
* THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY.
* IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS.
* ---------------------------------------------------------------------- */
/* First time flag use in sdastp for step smoothing logic. */
Info[IFIRST] = 0;
/* FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO
* ARE valid. */
*idid = -11;
for (i = 2; i <= 16; i++)
{
if ((Info[i] > infobd[i-(2)]) || (Info[i] < 0))
{
if ((i != 5) || (Info[i] < -infobd[i-(2)]))
{
if (i == IORD)
*idid = -13;
goto L_500;
}
}
}
/* Reverse communication controller is cleared and set iwork(ldelt) */
Iwork[REVLOC] = 0;
Iwork[LDELT] = LDELTA;
Iwork[LKOLD] = 0;
Iwork[LK] = 1;
if (neq <= 0)
{
*idid = -12;
goto L_510;
}
/* CHECK LENGTH IWORK */
leniw = 30 + neq;
if (liw < leniw)
{
*idid = -15;
goto L_510;
}
/* Reserve space for constraints (if any) */
Iwork[LCNSTR] = Info[ICNSTR];
/* May be changed later if Jacobian is banded. */
ldd = neq + Iwork[LCNSTR];
/* May reset maximum number of steps between output points. */
Iwork[MXSTEP] = Info[IXSTEP];
if (Iwork[MXSTEP] <= 0)
Iwork[MXSTEP] = 510;
/* Set maximum order */
if (Info[IORD] == 0)
{
Iwork[LMXORD] = 5;
}
else
{
Iwork[LMXORD] = Info[IORD];
}
/* Compute lenpd,lenrw.check ml and mu. */
Iwork[LMAT] = Info[IMAT];
i = labs( Iwork[LMAT] );
nomat = (i == 5) || (i == 6);
lenrw = (Iwork[LMXORD] + 4 + Iwork[LCNSTR])*neq + 45;
if (i > 4)
{
lenpd = 0;
if (i <= 6)
goto L_70;
i -= 6;
if (i > 4)
i -= 4;
}
/* User is not taking care of linear algebra. */
if (i < 3)
{
/* The dense case */
lenpd = neq*(neq + Iwork[LCNSTR]);
lenrw += lenpd;
}
else
{
/* The banded case */
if (Iwork[LCNSTR] != 0)
{
*idid = -27;
goto L_510;
}
if ((((Iwork[LML] < 0) || (Iwork[LML] >= neq)) || (Iwork[LMU] <
0)) || (Iwork[LMU] >= neq))
{
*idid = -25;
goto L_510;
}
ldd = 2*Iwork[LML] + Iwork[LMU] + 1;
lenpd = ldd*neq;
if (i == 3)
{
/* Derivatives computed with finite differences. */
mband = Iwork[LML] + Iwork[LMU] + 1;
msave = (neq/mband) + 1;
lenrw += lenpd + 2*msave;
}
else
{
lenrw += lenpd;
}
}
/* CHECK LENGTH OF RWORK. */
L_70:
;
Iwork[LNPD] = lenpd;
if (lrw < lenrw)
{
*idid = -14;
goto L_510;
}
/* Check consistent logic for storage of derivatives
* and performing the linear algebra.
*
* CHECK TO SEE THAT TOUT IS DIFFERENT FROM T */
if (tout == *t)
{
*idid = -26;
goto L_510;
}
/* CHECK HMAX */
if (Info[IMAXH] != 0)
{
hmax = Rwork[LHMAX];
if (hmax <= 0.0e0)
{
*idid = -20;
goto L_510;
}
}
/* INITIALIZE COUNTERS */
Iwork[LNST] = 0;
Iwork[LNRE] = 0;
Iwork[LNJE] = 0;
Iwork[LNSTL] = 0;
*idid = 1;
Info[1] = 1;
/* Flag that we are not completely initialized. */
Rwork[LROUND] = 0.e0;
/* A call to the evaluator for setting up or initializing.
* If reverse communication is used this is either
* a benign call to the dummy routine or else to the user's code, */
ires = 0;
(*sdasf)( t, y, yprime, rwork, rwork, &ldd, &Rwork[LCJ], &ires,
rwork, iwork );
goto L_110;
}
/* ----------------------------------------------------------------------
* THIS BLOCK IS FOR CONTINUATION CALLS
* ONLY. HERE WE CHECK INFO(1), AND IF THE
* LAST STEP WAS INTERRUPTED WE CHECK WHETHER
* APPROPRIATE ACTION WAS TAKEN.
* ----------------------------------------------------------------------
* */
if (Info[1] != 1)
{
i = 1;
if (Info[1] != -1)
{
*idid = -11;
goto L_500;
}
/* IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED
* BY AN ERROR CONDITION FROM SDASTP, AND
* APPROPRIATE ACTION WAS NOT TAKEN. THIS
* IS A FATAL ERROR. */
Idat[1] = *idid;
*idid = -31;
goto L_520;
}
if (Iwork[REVLOC] > 0)
{
/* Using reverse communication */
locate = Iwork[REVLOC]%8;
Iwork[REVLOC] /= 8;
switch (locate)
{
case 1: goto L_180;
case 2: goto L_380;
}
}
Iwork[LNSTL] = Iwork[LNST];
/* ----------------------------------------------------------------------
* THIS BLOCK IS EXECUTED ON ALL CALLS.
* THE ERROR TOLERANCE PARAMETERS ARE
* CHECKED, AND THE WORK ARRAY POINTERS
* ARE SET.
* ----------------------------------------------------------------------
* */
L_110:
;
if (Info[IDB] != 0)
sdasdb( 0, neq, *t, y, yprime, info, rwork, iwork, *idid,
atol, rtol );
/* CHECK RTOL,ATOL */
nzflg = 0;
itemp = 1;
if (Info[ITOL] == 1)
itemp = neq;
for (i = 1; i <= itemp; i++)
{
rtoli = Rtol[i];
atoli = Atol[i];
if (rtoli > 0.0e0 || atoli > 0.0e0)
nzflg = 1;
if (rtoli < 0.0e0)
{
*idid = -16;
goto L_510;
}
if (atoli < 0.0e0)
{
*idid = -17;
goto L_510;
}
}
if (nzflg == 0)
{
*idid = -18;
goto L_510;
}
/* SET UP RWORK STORAGE.IWORK STORAGE IS FIXED
* IN DATA STATEMENT. */
Iwork[LE] = LDELTA + neq + Iwork[LCNSTR];
Iwork[LWT] = Iwork[LE] + neq;
Iwork[LPHI] = Iwork[LWT] + neq;
if (nomat)
{
Iwork[LWM] = lrw;
}
else
{
Iwork[LWM] = Iwork[LPHI] + (Iwork[LMXORD] + 1)*neq;
}
Iwork[NTEMP] = 1 + Iwork[LNPD];
if (Rwork[LROUND] != 0.e0)
goto L_220;
/* ----------------------------------------------------------------------
* THIS BLOCK IS EXECUTED ON THE INITIAL CALL
* ONLY. SET THE INITIAL STEP SIZE, AND
* THE ERROR WEIGHT VECTOR, AND PHI.
* COMPUTE INITIAL YPRIME, IF NECESSARY.
* ----------------------------------------------------------------------
* */
tn = *t;
*idid = 1;
/* SET ERROR WEIGHT VECTOR WT. Pass info() to routine for
* info(itol), info(ismoot) */
sdaswt( neq, info, rtol, atol, y, &Rwork[Iwork[LWT]], rwork, iwork );
for (i = 1; i <= neq; i++)
{
if (Rwork[Iwork[LWT] + i - 1] <= 0.0e0)
{
*idid = -5;
goto L_510;
}
}
/* COMPUTE UNIT ROUNDOFF AND HMIN */
Rwork[LROUND] = FLT_EPSILON;
Rwork[LHMIN] = 16.0e0*Rwork[LROUND]*fmaxf( fabsf( *t ), fabsf( tout ) );
/* CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH */
tdist = fabsf( tout - *t );
if (tdist < Rwork[LHMIN])
{
*idid = -23;
goto L_510;
}
/* CHECK HO, IF THIS WAS INPUT */
if (Info[IH0] != 0)
{
ho = Rwork[LH];
h = ho;
if ((tout - *t)*h < 0.0e0)
{
*idid = -21;
goto L_510;
}
if (h == 0.0e0)
{
*idid = -22;
goto L_510;
}
}
else
{
/* COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER
* SDASTP OR SDAS1, DEPENDING ON INFO(inityp) */
ho = 0.001e0*tdist;
ypnorm = sdasnm( neq, yprime, &Rwork[Iwork[LWT]], rwork, iwork );
if (ypnorm > 0.5e0/ho)
ho = 0.5e0/ypnorm;
ho = signf( ho, tout - *t );
}
/* ADJUST HO IF NECESSARY TO MEET HMAX BOUND */
if (Info[IMAXH] == 0)
goto L_160;
rh = fabsf( ho )/Rwork[LHMAX];
if (rh > 1.0e0)
ho /= rh;
/* COMPUTE TSTOP, IF APPLICABLE */
L_160:
if (Info[ISTOP] != 0)
{
tstop = Rwork[LTSTOP];
if ((tstop - *t)*ho < 0.0e0)
{
*idid = -24;
goto L_510;
}
if ((*t + ho - tstop)*ho > 0.0e0)
ho = tstop - *t;
if ((tstop - tout)*ho < 0.0e0)
{
*idid = -19;
goto L_510;
}
}
/* COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE */
Rwork[LTN] = tn;
/* Set initial error estimates to 0 (temporary??) */
for (i = Iwork[LE]; i <= (Iwork[LWT] - 1); i++)
{
Rwork[i] = 0.e0;
}
/* If info(inityp)==0, user has provided initial values of y
* and yprime. */
if (Info[INITYP] == 0)
goto L_190;
/* REVERSE ENTRY 1: */
L_180:
;
/* If info(inityp)==1, user has provided initial values of y
* and a defined value of yprime. Use Petzold's starting method.
* The preferred starting method uses info(inityp)==2. */
if (Info[INITYP] == 1)
{
sdas1( &Rwork[LTN], y, yprime, neq, &ldd, sdasf, info, &ho,
&Rwork[Iwork[LWT]], idid, &Rwork[Iwork[LPHI]], &Rwork[LDELTA],
&Rwork[Iwork[LE]], &Rwork[Iwork[LWM]], iwork, rwork );
}
else if (Info[INITYP] == 2)
{
/* Use algorithm of Krogh/Hanson based on Newton method for the
* starting values of YPRIME. */
printf("This should check that F(t,y,y') is small\n");
}
tn = Rwork[LTN];
/* See if reverse communication needed: */
if (Iwork[REVLOC] != 0)
{
Iwork[REVLOC] = 8*Iwork[REVLOC] + 1;
goto L_490;
}
if (*idid < 0)
goto L_510;
/* LOAD H WITH HO. STORE H IN RWORK(LH) */
L_190:
h = ho;
Rwork[LH] = h;
/* Flag that we are starting for sdastp. */
Iwork[REVLOC] = -1;
/* LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2) */
itemp = Iwork[LPHI] + neq;
for (i = 1; i <= neq; i++)
{
Rwork[Iwork[LPHI] + i - 1] = Y[i];
Rwork[itemp + i - 1] = h*Yprime[i];
}
goto L_320;
/* ------------------------------------------------------
* THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS
* PURPOSE IS TO CHECK STOP CONDITIONS BEFORE
* TAKING A STEP.
* ADJUST H IF NECESSARY TO MEET HMAX BOUND
* ------------------------------------------------------
* */
L_220:
;
done = FALSE;
tn = Rwork[LTN];
h = Rwork[LH];
if (Info[IMAXH] == 0)
goto L_230;
rh = fabsf( h )/Rwork[LHMAX];
if (rh > 1.0e0)
h /= rh;
L_230:
;
if (*t == tout)
{
*idid = -26;
goto L_510;
}
if ((*t - tout)*h > 0.0e0)
{
*idid = -21;
goto L_510;
}
if (Info[ISTOP] == 1)
goto L_260;
if (Info[IOUT] == 1)
goto L_240;
if ((tn - tout)*h < 0.0e0)
goto L_310;
sdasin( tn, tout, y, yprime, neq, Iwork[LKOLD], &Rwork[Iwork[LPHI]],
&Rwork[LPSI] );
*t = tout;
*idid = 3;
done = TRUE;
goto L_480;
L_240:
if ((tn - *t)*h <= 0.0e0)
goto L_310;
if ((tn - tout)*h > 0.0e0)
goto L_250;
sdasin( tn, tn, y, yprime, neq, Iwork[LKOLD], &Rwork[Iwork[LPHI]],
&Rwork[LPSI] );
*t = tn;
*idid = 1;
done = TRUE;
goto L_480;
L_250:
;
sdasin( tn, tout, y, yprime, neq, Iwork[LKOLD], &Rwork[Iwork[LPHI]],
&Rwork[LPSI] );
*t = tout;
*idid = 3;
done = TRUE;
goto L_480;
L_260:
if (Info[IOUT] == 1)
goto L_270;
tstop = Rwork[LTSTOP];
if ((tn - tstop)*h > 0.0e0)
{
*idid = -24;
goto L_510;
}
if ((tstop - tout)*h < 0.0e0)
{
*idid = -19;
goto L_510;
}
if ((tn - tout)*h < 0.0e0)
goto L_290;
sdasin( tn, tout, y, yprime, neq, Iwork[LKOLD], &Rwork[Iwork[LPHI]],
&Rwork[LPSI] );
*t = tout;
*idid = 3;
done = TRUE;
goto L_480;
L_270:
tstop = Rwork[LTSTOP];
if ((tn - tstop)*h > 0.0e0)
{
*idid = -24;
goto L_510;
}
if ((tstop - tout)*h < 0.0e0)
{
*idid = -19;
goto L_510;
}
if ((tn - *t)*h <= 0.0e0)
goto L_290;
if ((tn - tout)*h > 0.0e0)
goto L_280;
sdasin( tn, tn, y, yprime, neq, Iwork[LKOLD], &Rwork[Iwork[LPHI]],
&Rwork[LPSI] );
*t = tn;
*idid = 1;
done = TRUE;
goto L_310;
L_280:
;
sdasin( tn, tout, y, yprime, neq, Iwork[LKOLD], &Rwork[Iwork[LPHI]],
&Rwork[LPSI] );
*t = tout;
*idid = 3;
done = TRUE;
goto L_310;
L_290:
;
/* CHECK WHETHER WE ARE WITHIN ROUNDOFF OF TSTOP */
if (fabsf( tn - tstop ) > 100.0e0*Rwork[LROUND]*(fabsf( tn ) + fabsf( h )))
goto L_300;
sdasin( tn, tstop, y, yprime, neq, Iwork[LKOLD], &Rwork[Iwork[LPHI]],
&Rwork[LPSI] );
*idid = 2;
*t = tstop;
done = TRUE;
goto L_310;
L_300:
tnext = tn + h;
if ((tnext - tstop)*h <= 0.0e0)
goto L_310;
h = tstop - tn;
Rwork[LH] = h;
L_310:
if (done)
goto L_480;
/* ------------------------------------------------------
* THE NEXT BLOCK CONTAINS THE CALL TO THE
* ONE-STEP INTEGRATOR SDASTP.
* THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS.
* CHECK FOR TOO MANY STEPS.
* UPDATE WT.
* CHECK FOR TOO MUCH ACCURACY REQUESTED.
* COMPUTE MINIMUM STEPSIZE.
* ------------------------------------------------------
* */
L_320:
;
/* CHECK FOR TOO MANY STEPS */
if ((Iwork[LNST] - Iwork[LNSTL]) < Iwork[MXSTEP])
goto L_330;
*idid = -3;
goto L_390;
/* UPDATE WT. Pass info() to routine for
* info(itol), info(ismoot) */
L_330:
sdaswt( neq, info, rtol, atol, &Rwork[Iwork[LPHI]], &Rwork[Iwork[LWT]],
rwork, iwork );
for (i = 1; i <= neq; i++)
{
if (Rwork[i + Iwork[LWT] - 1] > 0.0e0)
goto L_340;
*idid = -5;
goto L_390;
L_340:
;
}
/* TEST FOR TOO MUCH ACCURACY REQUESTED. */
r = sdasnm( neq, &Rwork[Iwork[LPHI]], &Rwork[Iwork[LWT]], rwork,
iwork )*100.0e0*Rwork[LROUND];
if (r > 1.0e0)
{
/* MULTIPLY RTOL AND ATOL BY R AND RETURN */
if (Info[ITOL] == 0)
{
Rtol[1] *= r;
Atol[1] *= r;
}
else
{
for (i = 1; i <= neq; i++)
{
Rtol[i] *= r;
Atol[i] *= r;
}
}
*idid = -2;
goto L_510;
}
/* COMPUTE MINIMUM STEPSIZE */
Rwork[LHMIN] = 4.0e0*Rwork[LROUND]*fmaxf( fabsf( tn ), fabsf( tout ) );
/* TEST H VS. HMAX */
if (Info[IMAXH] != 0)
{
rh = fabsf( h )/Rwork[LHMAX];
if (rh > 1.0e0)
h /= rh;
}
/* REVERSE ENTRY 2: */
L_380:
;
sdastp( &Rwork[LTN], y, yprime, neq, &ldd, sdasf, info, &h, &Rwork[Iwork[LWT]],
idid, &Rwork[Iwork[LPHI]], &Rwork[LDELTA], &Rwork[Iwork[LE]],
&Rwork[Iwork[LWM]], iwork, rwork, &Rwork[LALPHA], &Rwork[LBETA],
&Rwork[LGAMMA], &Rwork[LPSI], &Rwork[LSIGMA], &Iwork[LK] );
tn = Rwork[LTN];
/* See if reverse communication needed: */
if (Iwork[REVLOC] != 0)
{
Iwork[REVLOC] = 8*Iwork[REVLOC] + 2;
goto L_490;
}
L_390:
if (*idid < 0)
goto L_510;
/* -------------------------------------------------------
* THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN
* FROM SDASTP (IDID=1). TEST FOR STOP CONDITIONS.
* -------------------------------------------------------
* */
if (Info[ISTOP] != 0)
goto L_420;
if (Info[IOUT] != 0)
goto L_410;
if ((tn - tout)*h < 0.0e0)
goto L_320;
L_400:
sdasin( tn, tout, y, yprime, neq, Iwork[LKOLD], &Rwork[Iwork[LPHI]],
&Rwork[LPSI] );
*idid = 3;
*t = tout;
goto L_480;
L_410:
if ((tn - tout)*h >= 0.0e0)
goto L_400;
*t = tn;
*idid = 1;
goto L_480;
L_420:
if (Info[IOUT] != 0)
goto L_450;
if ((tn - tout)*h >= 0.0e0)
goto L_400;
if (fabsf( tn - tstop ) <= 100.0e0*Rwork[LROUND]*(fabsf( tn ) +
fabsf( h )))
goto L_440;
tnext = tn + h;
if ((tnext - tstop)*h <= 0.0e0)
goto L_320;
h = tstop - tn;
goto L_320;
L_440:
sdasin( tn, tstop, y, yprime, neq, Iwork[LKOLD], &Rwork[Iwork[LPHI]],
&Rwork[LPSI] );
*idid = 2;
*t = tstop;
goto L_480;
L_450:
if ((tn - tout)*h >= 0.0e0)
goto L_400;
if (fabsf( tn - tstop ) <= 100.0e0*Rwork[LROUND]*(fabsf( tn ) +
fabsf( h )))
goto L_440;
*t = tn;
*idid = 1;
goto L_480;
/* -------------------------------------------------------
* ALL SUCCESSFUL RETURNS FROM SDASLX ARE MADE FROM
* THIS BLOCK.
* -------------------------------------------------------
* */
L_480:
;
Rwork[LH] = h;
if (Info[IDB] != 0)
sdasdb( 1, neq, *t, y, yprime, info, rwork, iwork, *idid,
atol, rtol );
return;
/* ----------------------------------------------------------------------
* THIS BLOCK HANDLES ALL REVERSE COMMUNICATION RETURNS
* ---------------------------------------------------------------------- */
L_490:
;
*t = Rwork[LTN];
*idid = 4;
return;
/* One central place to call SMESS to print error messages.: */
L_500:
;
Idat[5] = i;
Idat[6] = Info[i];
L_510:
;
Idat[1] = Iwork[LKOLD];
Idat[2] = Iwork[LK];
Idat[3] = Iwork[MXSTEP];
Idat[4] = Iwork[LMXORD];
Idat[7] = lrw;
Idat[8] = lenrw;
Idat[9] = liw;
Idat[10] = leniw;
Idat[11] = Iwork[1];
Idat[12] = Iwork[2];
Idat[13] = neq;
Idat[14] = *idid;
Edat[1] = tn;
Edat[2] = h;
if (*idid <= -29)
{
if (*idid == -29)
{
Mact[5] = MECONT;
}
else if (*idid == -30)
{
Edat[3] = Rwork[LDELTA];
Edat[4] = Rwork[LDELTA + 1];
}
}
else
{
Edat[3] = tout;
Edat[4] = hmax;
Edat[5] = tstop;
}
L_520:
;
Mact[2] = Maperr[-*idid];
Mact[3] = -*idid;
Mact[4] = Mloc[1 - *idid];
smess( mact, (char*)mtxtaa,232, idat, edat );
Info[1] = -1;
if (*idid == -29)
{
Mact[5] = MERET;
Mactv[5] = neq;
Mactv[6] = MECONT;
smess( mactv, (char*)mtxtab,67, idat, &Rwork[LDELTA] );
Mactv[5] = Info[10];
Mactv[6] = MERET;
smess( mactv, (char*)mtxtac,39, idat, &Rwork[LDELTA + neq] );
}
goto L_480;
/* ----------END OF SUBROUTINE SDASLX------------------------------------ */
} /* end of function */