/*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 #include "fcrt.h" #include "sdaslx.h" #include #include #include #include /* 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 */