COPASI API  4.16.103
Public Member Functions | Private Attributes | Static Private Attributes | List of all members
CLSODAR Class Reference

#include <CLSODAR.h>

Inheritance diagram for CLSODAR:
Inheritance graph
[legend]
Collaboration diagram for CLSODAR:
Collaboration graph
[legend]

Public Member Functions

 CLSODAR ()
 
C_INT operator() (evalF f, C_INT *neq, double *y, double *t, double *tout, C_INT *itol, double *rtol, double *atol, C_INT *itask, C_INT *istate, C_INT *iopt, double *rwork, C_INT *lrw, C_INT *iwork, C_INT *liw, evalJ jac, C_INT *jt, evalG g, C_INT *ng, C_INT *jroot)
 
 ~CLSODAR ()
 
- Public Member Functions inherited from CInternalSolver
void enablePrint (const bool &print=true)
 
void resetState ()
 
void saveState ()
 
void setOstream (std::ostream &os)
 
 ~CInternalSolver ()
 

Private Attributes

PJACmpPJAC
 
SLVSmpSLVS
 

Static Private Attributes

static const C_INT mord [2] = {12, 5}
 
static const C_INT mxhnl0 = 10
 
static const C_INT mxstp0 = 500
 

Additional Inherited Members

- Protected Member Functions inherited from CInternalSolver
 CInternalSolver ()
 
C_INT dintdy_ (double *t, const C_INT *k, double *yh, C_INT *nyh, double *dky, C_INT *iflag)
 
C_INT dprja_ (C_INT *neq, double *y, double *yh, C_INT *nyh, double *ewt, double *ftem, double *savf, double *wm, C_INT *iwm, evalF f, evalJ jac)
 
C_INT drchek_ (const C_INT *job, evalG g, C_INT *neq, double *y, double *yh, C_INT *nyh, double *g0, double *g1, double *gx, C_INT *jroot, C_INT *irt)
 
C_INT droots_ (C_INT *ng, double *hmin, C_INT *jflag, double *x0, double *x1, double *g0, double *g1, double *gx, double *x, C_INT *jroot)
 
C_INT dsolsy_ (double *wm, C_INT *iwm, double *x, double *tem)
 
C_INT dstoda_ (C_INT *neq, double *y, double *yh, C_INT *nyh, double *yh1, double *ewt, double *savf, double *acor, double *wm, C_INT *iwm, evalF f, evalJ jac, PJAC *pjac, SLVS *slvs)
 
- Protected Attributes inherited from CInternalSolver
dls001 mdls001_
 
dlsa01 mdlsa01_
 
dlsr01 mdlsr01_
 
State mState
 
Cxerrwd mxerrwd
 

Detailed Description

Definition at line 24 of file CLSODAR.h.

Constructor & Destructor Documentation

CLSODAR::CLSODAR ( )

Definition at line 100 of file CLSODAR.cpp.

References CInternalSolver::dprja_(), CInternalSolver::dsolsy_(), mpPJAC, and mpSLVS.

100  :
101  CInternalSolver(),
102  mpPJAC(NULL),
103  mpSLVS(NULL)
104 {
107 }
SLVS * mpSLVS
Definition: CLSODAR.h:53
C_INT dsolsy_(double *wm, C_INT *iwm, double *x, double *tem)
Definition: dsolsy.cpp:43
PJAC * mpPJAC
Definition: CLSODAR.h:52
C_INT dprja_(C_INT *neq, double *y, double *yh, C_INT *nyh, double *ewt, double *ftem, double *savf, double *wm, C_INT *iwm, evalF f, evalJ jac)
Definition: dprja.cpp:47
CLSODAR::~CLSODAR ( )

Definition at line 109 of file CLSODAR.cpp.

References mpPJAC, and mpSLVS.

110 {
111  if (mpPJAC != NULL)
112  {
113  delete mpPJAC; mpPJAC = NULL;
114  }
115 
116  if (mpSLVS != NULL)
117  {
118  delete mpSLVS; mpSLVS = NULL;
119  }
120 }
SLVS * mpSLVS
Definition: CLSODAR.h:53
PJAC * mpPJAC
Definition: CLSODAR.h:52

Member Function Documentation

C_INT CLSODAR::operator() ( evalF  f,
C_INT neq,
double *  y,
double *  t,
double *  tout,
C_INT itol,
double *  rtol,
double *  atol,
C_INT itask,
C_INT istate,
C_INT iopt,
double *  rwork,
C_INT lrw,
C_INT iwork,
C_INT liw,
evalJ  jac,
C_INT jt,
evalG  g,
C_INT ng,
C_INT jroot 
)

Definition at line 124 of file CLSODAR.cpp.

References c__0, c__1, c__10, c__101, c__102, c__103, c__104, c__105, c__106, c__107, c__11, c__12, c__13, c__14, c__15, c__16, c__17, c__18, c__19, c__2, c__20, c__201, c__202, c__203, c__204, c__205, c__206, c__207, c__21, c__22, c__23, c__24, c__25, c__26, c__27, c__28, c__29, c__3, c__30, c__303, c__31, c__32, c__4, c__40, c__5, c__50, c__6, c__60, c__7, c__8, c__9, c_b76, C_INT, d_sign(), dcopy_(), dewset_(), CInternalSolver::dintdy_(), dls001_1, dlsa01_1, dlsr01_1, dmnorm_(), CInternalSolver::drchek_(), CInternalSolver::dstoda_(), if(), max, min, mord, mpPJAC, mpSLVS, CInternalSolver::mxerrwd, mxhnl0, and mxstp0.

129 {
130  /* Initialized data */
131 
132  /* System generated locals */
133  C_INT i__1;
134  double d__1, d__2;
135 
136  /* Local variables */
137  C_INT i__;
138  double h0;
139  C_INT i1, i2;
140  double w0;
141  C_INT ml;
142  double rh;
143  C_INT mu;
144  double tp;
145  C_INT lf0;
146  double big;
147  C_INT kgo;
148  double ayi;
149  std::string msg;
150  double hmx;
151  C_INT irt;
152  double tol, sum;
153  C_INT len1, len2;
154  double hmax;
155  C_INT irfp;
156  bool ihit = false;
157  double ewti, size;
158  C_INT len1c, len1n, len1s, iflag;
159  double atoli;
160  C_INT leniw, lenwm = 0, lenyh, imxer;
161  double tcrit;
162  C_INT lenrw;
163  double rtoli, tdist, tolsf, tnext;
164  C_INT leniwc;
165  C_INT lenrwc;
166  C_INT lyhnew;
167 
168  /* ----------------------------------------------------------------------- */
169  /* This is the 12 November 2003 version of */
170  /* DLSODAR: Livermore Solver for Ordinary Differential Equations, with */
171  /* Automatic method switching for stiff and nonstiff problems, */
172  /* and with Root-finding. */
173 
174  /* This version is in double precision. */
175 
176  /* DLSODAR solves the initial value problem for stiff or nonstiff */
177  /* systems of first order ODEs, */
178  /* dy/dt = f(t,y) , or, in component form, */
179  /* dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ). */
180  /* At the same time, it locates the roots of any of a set of functions */
181  /* g(i) = g(i,t,y(1),...,y(NEQ)) (i = 1,...,ng). */
182 
183  /* This a variant version of the DLSODE package. It differs from it */
184  /* in two ways: */
185  /* (a) It switches automatically between stiff and nonstiff methods. */
186  /* This means that the user does not have to determine whether the */
187  /* problem is stiff or not, and the solver will automatically choose the */
188  /* appropriate method. It always starts with the nonstiff method. */
189  /* (b) It finds the root of at least one of a set of constraint */
190  /* functions g(i) of the independent and dependent variables. */
191  /* It finds only those roots for which some g(i), as a function */
192  /* of t, changes sign in the interval of integration. */
193  /* It then returns the solution at the root, if that occurs */
194  /* sooner than the specified stop condition, and otherwise returns */
195  /* the solution according the specified stop condition. */
196 
197  /* Authors: Alan C. Hindmarsh, */
198  /* Center for Applied Scientific Computing, L-561 */
199  /* Lawrence Livermore National Laboratory */
200  /* Livermore, CA 94551 */
201  /* and */
202  /* Linda R. Petzold */
203  /* Univ. of California at Santa Barbara */
204  /* Dept. of Computer Science */
205  /* Santa Barbara, CA 93106 */
206 
207  /* References: */
208  /* 1. Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE */
209  /* Solvers, in Scientific Computing, R. S. Stepleman et al. (Eds.), */
210  /* North-Holland, Amsterdam, 1983, pp. 55-64. */
211  /* 2. Linda R. Petzold, Automatic Selection of Methods for Solving */
212  /* Stiff and Nonstiff Systems of Ordinary Differential Equations, */
213  /* Siam J. Sci. Stat. Comput. 4 (1983), pp. 136-148. */
214  /* 3. Kathie L. Hiebert and Lawrence F. Shampine, Implicitly Defined */
215  /* Output Points for Solutions of ODEs, Sandia Report SAND80-0180, */
216  /* February 1980. */
217  /* ----------------------------------------------------------------------- */
218  /* Summary of Usage. */
219 
220  /* Communication between the user and the DLSODAR package, for normal */
221  /* situations, is summarized here. This summary describes only a subset */
222  /* of the full set of options available. See the full description for */
223  /* details, including alternative treatment of the Jacobian matrix, */
224  /* optional inputs and outputs, nonstandard options, and */
225  /* instructions for special situations. See also the example */
226  /* problem (with program and output) following this summary. */
227 
228  /* A. First provide a subroutine of the form: */
229  /* SUBROUTINE F (NEQ, T, Y, YDOT) */
230  /* DOUBLE PRECISION T, Y(*), YDOT(*) */
231  /* which supplies the vector function f by loading YDOT(i) with f(i). */
232 
233  /* B. Provide a subroutine of the form: */
234  /* SUBROUTINE G (NEQ, T, Y, NG, GOUT) */
235  /* DOUBLE PRECISION T, Y(*), GOUT(NG) */
236  /* which supplies the vector function g by loading GOUT(i) with */
237  /* g(i), the i-th constraint function whose root is sought. */
238 
239  /* C. Write a main program which calls Subroutine DLSODAR once for */
240  /* each point at which answers are desired. This should also provide */
241  /* for possible use of logical unit 6 for output of error messages by */
242  /* DLSODAR. On the first call to DLSODAR, supply arguments as follows: */
243  /* F = name of subroutine for right-hand side vector f. */
244  /* This name must be declared External in calling program. */
245  /* NEQ = number of first order ODEs. */
246  /* Y = array of initial values, of length NEQ. */
247  /* T = the initial value of the independent variable. */
248  /* TOUT = first point where output is desired (.ne. T). */
249  /* ITOL = 1 or 2 according as ATOL (below) is a scalar or array. */
250  /* RTOL = relative tolerance parameter (scalar). */
251  /* ATOL = absolute tolerance parameter (scalar or array). */
252  /* the estimated local error in y(i) will be controlled so as */
253  /* to be less than */
254  /* EWT(i) = RTOL*ABS(Y(i)) + ATOL if ITOL = 1, or */
255  /* EWT(i) = RTOL*ABS(Y(i)) + ATOL(i) if ITOL = 2. */
256  /* Thus the local error test passes if, in each component, */
257  /* either the absolute error is less than ATOL (or ATOL(i)), */
258  /* or the relative error is less than RTOL. */
259  /* Use RTOL = 0.0 for pure absolute error control, and */
260  /* use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative error */
261  /* control. Caution: actual (global) errors may exceed these */
262  /* local tolerances, so choose them conservatively. */
263  /* ITASK = 1 for normal computation of output values of y at t = TOUT. */
264  /* ISTATE = integer flag (input and output). Set ISTATE = 1. */
265  /* IOPT = 0 to indicate no optional inputs used. */
266  /* RWORK = real work array of length at least: */
267  /* 22 + NEQ * MAX(16, NEQ + 9) + 3*NG. */
268  /* See also Paragraph F below. */
269  /* LRW = declared length of RWORK (in user's dimension). */
270  /* IWORK = integer work array of length at least 20 + NEQ. */
271  /* LIW = declared length of IWORK (in user's dimension). */
272  /* JAC = name of subroutine for Jacobian matrix. */
273  /* Use a dummy name. See also Paragraph F below. */
274  /* JT = Jacobian type indicator. Set JT = 2. */
275  /* See also Paragraph F below. */
276  /* G = name of subroutine for constraint functions, whose */
277  /* roots are desired during the integration. */
278  /* This name must be declared External in calling program. */
279  /* NG = number of constraint functions g(i). If there are none, */
280  /* set NG = 0, and pass a dummy name for G. */
281  /* JROOT = integer array of length NG for output of root information. */
282  /* See next paragraph. */
283  /* Note that the main program must declare arrays Y, RWORK, IWORK, */
284  /* JROOT, and possibly ATOL. */
285 
286  /* D. The output from the first call (or any call) is: */
287  /* Y = array of computed values of y(t) vector. */
288  /* T = corresponding value of independent variable. This is */
289  /* TOUT if ISTATE = 2, or the root location if ISTATE = 3, */
290  /* or the farthest point reached if DLSODAR was unsuccessful. */
291  /* ISTATE = 2 or 3 if DLSODAR was successful, negative otherwise. */
292  /* 2 means no root was found, and TOUT was reached as desired. */
293  /* 3 means a root was found prior to reaching TOUT. */
294  /* -1 means excess work done on this call (perhaps wrong JT). */
295  /* -2 means excess accuracy requested (tolerances too small). */
296  /* -3 means illegal input detected (see printed message). */
297  /* -4 means repeated error test failures (check all inputs). */
298  /* -5 means repeated convergence failures (perhaps bad Jacobian */
299  /* supplied or wrong choice of JT or tolerances). */
300  /* -6 means error weight became zero during problem. (Solution */
301  /* component i vanished, and ATOL or ATOL(i) = 0.) */
302  /* -7 means work space insufficient to finish (see messages). */
303  /* -33 means One or more components of g have a root too near */
304  /* to the initial point. */
305  /* JROOT = array showing roots found if ISTATE = 3 on return. */
306  /* JROOT(i) = 1 if g(i) has a root at t, or 0 otherwise. */
307 
308  /* E. To continue the integration after a successful return, proceed */
309  /* as follows: */
310  /* (a) If ISTATE = 2 on return, reset TOUT and call DLSODAR again. */
311  /* (b) If ISTATE = 3 on return, reset ISTATE to 2, call DLSODAR again. */
312  /* In either case, no other parameters need be reset. */
313 
314  /* F. Note: If and when DLSODAR regards the problem as stiff, and */
315  /* switches methods accordingly, it must make use of the NEQ by NEQ */
316  /* Jacobian matrix, J = df/dy. For the sake of simplicity, the */
317  /* inputs to DLSODAR recommended in Paragraph C above cause DLSODAR to */
318  /* treat J as a full matrix, and to approximate it internally by */
319  /* difference quotients. Alternatively, J can be treated as a band */
320  /* matrix (with great potential reduction in the size of the RWORK */
321  /* array). Also, in either the full or banded case, the user can supply */
322  /* J in closed form, with a routine whose name is passed as the JAC */
323  /* argument. These alternatives are described in the paragraphs on */
324  /* RWORK, JAC, and JT in the full description of the call sequence below. */
325 
326  /* ----------------------------------------------------------------------- */
327  /* Example Problem. */
328 
329  /* The following is a simple example problem, with the coding */
330  /* needed for its solution by DLSODAR. The problem is from chemical */
331  /* kinetics, and consists of the following three rate equations: */
332  /* dy1/dt = -.04*y1 + 1.e4*y2*y3 */
333  /* dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2 */
334  /* dy3/dt = 3.e7*y2**2 */
335  /* on the interval from t = 0.0 to t = 4.e10, with initial conditions */
336  /* y1 = 1.0, y2 = y3 = 0. The problem is stiff. */
337  /* In addition, we want to find the values of t, y1, y2, and y3 at which */
338  /* (1) y1 reaches the value 1.e-4, and */
339  /* (2) y3 reaches the value 1.e-2. */
340 
341  /* The following coding solves this problem with DLSODAR, */
342  /* printing results at t = .4, 4., ..., 4.e10, and at the computed */
343  /* roots. It uses ITOL = 2 and ATOL much smaller for y2 than y1 or y3 */
344  /* because y2 has much smaller values. */
345  /* At the end of the run, statistical quantities of interest are */
346  /* printed (see optional outputs in the full description below). */
347 
348  /* EXTERNAL FEX, GEX */
349  /* DOUBLE PRECISION ATOL, RTOL, RWORK, T, TOUT, Y */
350  /* DIMENSION Y(3), ATOL(3), RWORK(76), IWORK(23), JROOT(2) */
351  /* NEQ = 3 */
352  /* Y(1) = 1. */
353  /* Y(2) = 0. */
354  /* Y(3) = 0. */
355  /* T = 0. */
356  /* TOUT = .4 */
357  /* ITOL = 2 */
358  /* RTOL = 1.D-4 */
359  /* ATOL(1) = 1.D-6 */
360  /* ATOL(2) = 1.D-10 */
361  /* ATOL(3) = 1.D-6 */
362  /* ITASK = 1 */
363  /* ISTATE = 1 */
364  /* IOPT = 0 */
365  /* LRW = 76 */
366  /* LIW = 23 */
367  /* JT = 2 */
368  /* NG = 2 */
369  /* DO 40 IOUT = 1,12 */
370  /* 10 CALL DLSODAR(FEX,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE, */
371  /* 1 IOPT,RWORK,LRW,IWORK,LIW,JDUM,JT,GEX,NG,JROOT) */
372  /* WRITE(6,20)T,Y(1),Y(2),Y(3) */
373  /* 20 FORMAT(' At t =',D12.4,' Y =',3D14.6) */
374  /* IF (ISTATE .LT. 0) GO TO 80 */
375  /* IF (ISTATE .EQ. 2) GO TO 40 */
376  /* WRITE(6,30)JROOT(1),JROOT(2) */
377  /* 30 FORMAT(5X,' The above line is a root, JROOT =',2I5) */
378  /* ISTATE = 2 */
379  /* GO TO 10 */
380  /* 40 TOUT = TOUT*10. */
381  /* WRITE(6,60)IWORK(11),IWORK(12),IWORK(13),IWORK(10), */
382  /* 1 IWORK(19),RWORK(15) */
383  /* 60 FORMAT(/' No. steps =',I4,' No. f-s =',I4,' No. J-s =',I4, */
384  /* 1 ' No. g-s =',I4/ */
385  /* 2 ' Method last used =',I2,' Last switch was at t =',D12.4) */
386  /* STOP */
387  /* 80 WRITE(6,90)ISTATE */
388  /* 90 FORMAT(///' Error halt.. ISTATE =',I3) */
389  /* STOP */
390  /* END */
391 
392  /* SUBROUTINE FEX (NEQ, T, Y, YDOT) */
393  /* DOUBLE PRECISION T, Y, YDOT */
394  /* DIMENSION Y(3), YDOT(3) */
395  /* YDOT(1) = -.04*Y(1) + 1.D4*Y(2)*Y(3) */
396  /* YDOT(3) = 3.D7*Y(2)*Y(2) */
397  /* YDOT(2) = -YDOT(1) - YDOT(3) */
398  /* RETURN */
399  /* END */
400 
401  /* SUBROUTINE GEX (NEQ, T, Y, NG, GOUT) */
402  /* DOUBLE PRECISION T, Y, GOUT */
403  /* DIMENSION Y(3), GOUT(2) */
404  /* GOUT(1) = Y(1) - 1.D-4 */
405  /* GOUT(2) = Y(3) - 1.D-2 */
406  /* RETURN */
407  /* END */
408 
409  /* The output of this program (on a CDC-7600 in single precision) */
410  /* is as follows: */
411 
412  /* At t = 2.6400e-01 y = 9.899653e-01 3.470563e-05 1.000000e-02 */
413  /* The above line is a root, JROOT = 0 1 */
414  /* At t = 4.0000e-01 Y = 9.851712e-01 3.386380e-05 1.479493e-02 */
415  /* At t = 4.0000e+00 Y = 9.055333e-01 2.240655e-05 9.444430e-02 */
416  /* At t = 4.0000e+01 Y = 7.158403e-01 9.186334e-06 2.841505e-01 */
417  /* At t = 4.0000e+02 Y = 4.505250e-01 3.222964e-06 5.494717e-01 */
418  /* At t = 4.0000e+03 Y = 1.831975e-01 8.941774e-07 8.168016e-01 */
419  /* At t = 4.0000e+04 Y = 3.898730e-02 1.621940e-07 9.610125e-01 */
420  /* At t = 4.0000e+05 Y = 4.936363e-03 1.984221e-08 9.950636e-01 */
421  /* At t = 4.0000e+06 Y = 5.161831e-04 2.065786e-09 9.994838e-01 */
422  /* At t = 2.0745e+07 Y = 1.000000e-04 4.000395e-10 9.999000e-01 */
423  /* The above line is a root, JROOT = 1 0 */
424  /* At t = 4.0000e+07 Y = 5.179817e-05 2.072032e-10 9.999482e-01 */
425  /* At t = 4.0000e+08 Y = 5.283401e-06 2.113371e-11 9.999947e-01 */
426  /* At t = 4.0000e+09 Y = 4.659031e-07 1.863613e-12 9.999995e-01 */
427  /* At t = 4.0000e+10 Y = 1.404280e-08 5.617126e-14 1.000000e+00 */
428 
429  /* No. steps = 361 No. f-s = 693 No. J-s = 64 No. g-s = 390 */
430  /* Method last used = 2 Last switch was at t = 6.0092e-03 */
431 
432  /* ----------------------------------------------------------------------- */
433  /* Full Description of User Interface to DLSODAR. */
434 
435  /* The user interface to DLSODAR consists of the following parts. */
436 
437  /* 1. The call sequence to Subroutine DLSODAR, which is a driver */
438  /* routine for the solver. This includes descriptions of both */
439  /* the call sequence arguments and of user-supplied routines. */
440  /* Following these descriptions is a description of */
441  /* optional inputs available through the call sequence, and then */
442  /* a description of optional outputs (in the work arrays). */
443 
444  /* 2. Descriptions of other routines in the DLSODAR package that may be */
445  /* (optionally) called by the user. These provide the ability to */
446  /* alter error message handling, save and restore the internal */
447  /* Common, and obtain specified derivatives of the solution y(t). */
448 
449  /* 3. Descriptions of Common blocks to be declared in overlay */
450  /* or similar environments, or to be saved when doing an interrupt */
451  /* of the problem and continued solution later. */
452 
453  /* 4. Description of a subroutine in the DLSODAR package, */
454  /* which the user may replace with his/her own version, if desired. */
455  /* this relates to the measurement of errors. */
456 
457  /* ----------------------------------------------------------------------- */
458  /* Part 1. Call Sequence. */
459 
460  /* The call sequence parameters used for input only are */
461  /* F, NEQ, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW, JAC, */
462  /* JT, G, and NG, */
463  /* that used only for output is JROOT, */
464  /* and those used for both input and output are */
465  /* Y, T, ISTATE. */
466  /* The work arrays RWORK and IWORK are also used for conditional and */
467  /* optional inputs and optional outputs. (The term output here refers */
468  /* to the return from Subroutine DLSODAR to the user's calling program.) */
469 
470  /* The legality of input parameters will be thoroughly checked on the */
471  /* initial call for the problem, but not checked thereafter unless a */
472  /* change in input parameters is flagged by ISTATE = 3 on input. */
473 
474  /* The descriptions of the call arguments are as follows. */
475 
476  /* F = the name of the user-supplied subroutine defining the */
477  /* ODE system. The system must be put in the first-order */
478  /* form dy/dt = f(t,y), where f is a vector-valued function */
479  /* of the scalar t and the vector y. Subroutine F is to */
480  /* compute the function f. It is to have the form */
481  /* SUBROUTINE F (NEQ, T, Y, YDOT) */
482  /* DOUBLE PRECISION T, Y(*), YDOT(*) */
483  /* where NEQ, T, and Y are input, and the array YDOT = f(t,y) */
484  /* is output. Y and YDOT are arrays of length NEQ. */
485  /* Subroutine F should not alter Y(1),...,Y(NEQ). */
486  /* F must be declared External in the calling program. */
487 
488  /* Subroutine F may access user-defined quantities in */
489  /* NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array */
490  /* (dimensioned in F) and/or Y has length exceeding NEQ(1). */
491  /* See the descriptions of NEQ and Y below. */
492 
493  /* If quantities computed in the F routine are needed */
494  /* externally to DLSODAR, an extra call to F should be made */
495  /* for this purpose, for consistent and accurate results. */
496  /* If only the derivative dy/dt is needed, use DINTDY instead. */
497 
498  /* NEQ = the size of the ODE system (number of first order */
499  /* ordinary differential equations). Used only for input. */
500  /* NEQ may be decreased, but not increased, during the problem. */
501  /* If NEQ is decreased (with ISTATE = 3 on input), the */
502  /* remaining components of Y should be left undisturbed, if */
503  /* these are to be accessed in F and/or JAC. */
504 
505  /* Normally, NEQ is a scalar, and it is generally referred to */
506  /* as a scalar in this user interface description. However, */
507  /* NEQ may be an array, with NEQ(1) set to the system size. */
508  /* (The DLSODAR package accesses only NEQ(1).) In either case, */
509  /* this parameter is passed as the NEQ argument in all calls */
510  /* to F, JAC, and G. Hence, if it is an array, locations */
511  /* NEQ(2),... may be used to store other integer data and pass */
512  /* it to F, JAC, and G. Each such subroutine must include */
513  /* NEQ in a Dimension statement in that case. */
514 
515  /* Y = a real array for the vector of dependent variables, of */
516  /* length NEQ or more. Used for both input and output on the */
517  /* first call (ISTATE = 1), and only for output on other calls. */
518  /* On the first call, Y must contain the vector of initial */
519  /* values. On output, Y contains the computed solution vector, */
520  /* evaluated at T. If desired, the Y array may be used */
521  /* for other purposes between calls to the solver. */
522 
523  /* This array is passed as the Y argument in all calls to F, */
524  /* JAC, and G. Hence its length may exceed NEQ, and locations */
525  /* Y(NEQ+1),... may be used to store other real data and */
526  /* pass it to F, JAC, and G. (The DLSODAR package accesses only */
527  /* Y(1),...,Y(NEQ).) */
528 
529  /* T = the independent variable. On input, T is used only on the */
530  /* first call, as the initial point of the integration. */
531  /* On output, after each call, T is the value at which a */
532  /* computed solution y is evaluated (usually the same as TOUT). */
533  /* If a root was found, T is the computed location of the */
534  /* root reached first, on output. */
535  /* On an error return, T is the farthest point reached. */
536 
537  /* TOUT = the next value of t at which a computed solution is desired. */
538  /* Used only for input. */
539 
540  /* When starting the problem (ISTATE = 1), TOUT may be equal */
541  /* to T for one call, then should .ne. T for the next call. */
542  /* For the initial T, an input value of TOUT .ne. T is used */
543  /* in order to determine the direction of the integration */
544  /* (i.e. the algebraic sign of the step sizes) and the rough */
545  /* scale of the problem. Integration in either direction */
546  /* (forward or backward in t) is permitted. */
547 
548  /* If ITASK = 2 or 5 (one-step modes), TOUT is ignored after */
549  /* the first call (i.e. the first call with TOUT .ne. T). */
550  /* Otherwise, TOUT is required on every call. */
551 
552  /* If ITASK = 1, 3, or 4, the values of TOUT need not be */
553  /* monotone, but a value of TOUT which backs up is limited */
554  /* to the current internal T interval, whose endpoints are */
555  /* TCUR - HU and TCUR (see optional outputs, below, for */
556  /* TCUR and HU). */
557 
558  /* ITOL = an indicator for the type of error control. See */
559  /* description below under ATOL. Used only for input. */
560 
561  /* RTOL = a relative error tolerance parameter, either a scalar or */
562  /* an array of length NEQ. See description below under ATOL. */
563  /* Input only. */
564 
565  /* ATOL = an absolute error tolerance parameter, either a scalar or */
566  /* an array of length NEQ. Input only. */
567 
568  /* The input parameters ITOL, RTOL, and ATOL determine */
569  /* the error control performed by the solver. The solver will */
570  /* control the vector E = (E(i)) of estimated local errors */
571  /* in y, according to an inequality of the form */
572  /* max-norm of (E(i)/EWT(i)) .le. 1, */
573  /* where EWT = (EWT(i)) is a vector of positive error weights. */
574  /* The values of RTOL and ATOL should all be non-negative. */
575  /* The following table gives the types (scalar/array) of */
576  /* RTOL and ATOL, and the corresponding form of EWT(i). */
577 
578  /* ITOL RTOL ATOL EWT(i) */
579  /* 1 scalar scalar RTOL*ABS(Y(i)) + ATOL */
580  /* 2 scalar array RTOL*ABS(Y(i)) + ATOL(i) */
581  /* 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL */
582  /* 4 array array RTOL(i)*ABS(Y(i)) + ATOL(i) */
583 
584  /* When either of these parameters is a scalar, it need not */
585  /* be dimensioned in the user's calling program. */
586 
587  /* If none of the above choices (with ITOL, RTOL, and ATOL */
588  /* fixed throughout the problem) is suitable, more general */
589  /* error controls can be obtained by substituting a */
590  /* user-supplied routine for the setting of EWT. */
591  /* See Part 4 below. */
592 
593  /* If global errors are to be estimated by making a repeated */
594  /* run on the same problem with smaller tolerances, then all */
595  /* components of RTOL and ATOL (i.e. of EWT) should be scaled */
596  /* down uniformly. */
597 
598  /* ITASK = an index specifying the task to be performed. */
599  /* input only. ITASK has the following values and meanings. */
600  /* 1 means normal computation of output values of y(t) at */
601  /* t = TOUT (by overshooting and interpolating). */
602  /* 2 means take one step only and return. */
603  /* 3 means stop at the first internal mesh point at or */
604  /* beyond t = TOUT and return. */
605  /* 4 means normal computation of output values of y(t) at */
606  /* t = TOUT but without overshooting t = TCRIT. */
607  /* TCRIT must be input as RWORK(1). TCRIT may be equal to */
608  /* or beyond TOUT, but not behind it in the direction of */
609  /* integration. This option is useful if the problem */
610  /* has a singularity at or beyond t = TCRIT. */
611  /* 5 means take one step, without passing TCRIT, and return. */
612  /* TCRIT must be input as RWORK(1). */
613 
614  /* Note: If ITASK = 4 or 5 and the solver reaches TCRIT */
615  /* (within roundoff), it will return T = TCRIT (exactly) to */
616  /* indicate this (unless ITASK = 4 and TOUT comes before TCRIT, */
617  /* in which case answers at t = TOUT are returned first). */
618 
619  /* ISTATE = an index used for input and output to specify the */
620  /* the state of the calculation. */
621 
622  /* On input, the values of ISTATE are as follows. */
623  /* 1 means this is the first call for the problem */
624  /* (initializations will be done). See note below. */
625  /* 2 means this is not the first call, and the calculation */
626  /* is to continue normally, with no change in any input */
627  /* parameters except possibly TOUT and ITASK. */
628  /* (If ITOL, RTOL, and/or ATOL are changed between calls */
629  /* with ISTATE = 2, the new values will be used but not */
630  /* tested for legality.) */
631  /* 3 means this is not the first call, and the */
632  /* calculation is to continue normally, but with */
633  /* a change in input parameters other than */
634  /* TOUT and ITASK. Changes are allowed in */
635  /* NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, JT, ML, MU, */
636  /* and any optional inputs except H0, MXORDN, and MXORDS. */
637  /* (See IWORK description for ML and MU.) */
638  /* In addition, immediately following a return with */
639  /* ISTATE = 3 (root found), NG and G may be changed. */
640  /* (But changing NG from 0 to .gt. 0 is not allowed.) */
641  /* Note: A preliminary call with TOUT = T is not counted */
642  /* as a first call here, as no initialization or checking of */
643  /* input is done. (Such a call is sometimes useful for the */
644  /* purpose of outputting the initial conditions.) */
645  /* Thus the first call for which TOUT .ne. T requires */
646  /* ISTATE = 1 on input. */
647 
648  /* On output, ISTATE has the following values and meanings. */
649  /* 1 means nothing was done; TOUT = t and ISTATE = 1 on input. */
650  /* 2 means the integration was performed successfully, and */
651  /* no roots were found. */
652  /* 3 means the integration was successful, and one or more */
653  /* roots were found before satisfying the stop condition */
654  /* specified by ITASK. See JROOT. */
655  /* -1 means an excessive amount of work (more than MXSTEP */
656  /* steps) was done on this call, before completing the */
657  /* requested task, but the integration was otherwise */
658  /* successful as far as T. (MXSTEP is an optional input */
659  /* and is normally 500.) To continue, the user may */
660  /* simply reset ISTATE to a value .gt. 1 and call again */
661  /* (the excess work step counter will be reset to 0). */
662  /* In addition, the user may increase MXSTEP to avoid */
663  /* this error return (see below on optional inputs). */
664  /* -2 means too much accuracy was requested for the precision */
665  /* of the machine being used. This was detected before */
666  /* completing the requested task, but the integration */
667  /* was successful as far as T. To continue, the tolerance */
668  /* parameters must be reset, and ISTATE must be set */
669  /* to 3. The optional output TOLSF may be used for this */
670  /* purpose. (Note: If this condition is detected before */
671  /* taking any steps, then an illegal input return */
672  /* (ISTATE = -3) occurs instead.) */
673  /* -3 means illegal input was detected, before taking any */
674  /* integration steps. See written message for details. */
675  /* Note: If the solver detects an infinite loop of calls */
676  /* to the solver with illegal input, it will cause */
677  /* the run to stop. */
678  /* -4 means there were repeated error test failures on */
679  /* one attempted step, before completing the requested */
680  /* task, but the integration was successful as far as T. */
681  /* The problem may have a singularity, or the input */
682  /* may be inappropriate. */
683  /* -5 means there were repeated convergence test failures on */
684  /* one attempted step, before completing the requested */
685  /* task, but the integration was successful as far as T. */
686  /* This may be caused by an inaccurate Jacobian matrix, */
687  /* if one is being used. */
688  /* -6 means EWT(i) became zero for some i during the */
689  /* integration. Pure relative error control (ATOL(i)=0.0) */
690  /* was requested on a variable which has now vanished. */
691  /* The integration was successful as far as T. */
692  /* -7 means the length of RWORK and/or IWORK was too small to */
693  /* proceed, but the integration was successful as far as T. */
694  /* This happens when DLSODAR chooses to switch methods */
695  /* but LRW and/or LIW is too small for the new method. */
696  /* -33 means One or more components of g have a root too near */
697  /* to the initial point. This happens when a root is constant */
698  /* and zero. */
699 
700  /* Note: Since the normal output value of ISTATE is 2, */
701  /* it does not need to be reset for normal continuation. */
702  /* Also, since a negative input value of ISTATE will be */
703  /* regarded as illegal, a negative output value requires the */
704  /* user to change it, and possibly other inputs, before */
705  /* calling the solver again. */
706 
707  /* IOPT = an integer flag to specify whether or not any optional */
708  /* inputs are being used on this call. Input only. */
709  /* The optional inputs are listed separately below. */
710  /* IOPT = 0 means no optional inputs are being used. */
711  /* Default values will be used in all cases. */
712  /* IOPT = 1 means one or more optional inputs are being used. */
713 
714  /* RWORK = a real array (double precision) for work space, and (in the */
715  /* first 20 words) for conditional and optional inputs and */
716  /* optional outputs. */
717  /* As DLSODAR switches automatically between stiff and nonstiff */
718  /* methods, the required length of RWORK can change during the */
719  /* problem. Thus the RWORK array passed to DLSODAR can either */
720  /* have a static (fixed) length large enough for both methods, */
721  /* or have a dynamic (changing) length altered by the calling */
722  /* program in response to output from DLSODAR. */
723 
724  /* --- Fixed Length Case --- */
725  /* If the RWORK length is to be fixed, it should be at least */
726  /* max (LRN, LRS), */
727  /* where LRN and LRS are the RWORK lengths required when the */
728  /* current method is nonstiff or stiff, respectively. */
729 
730  /* The separate RWORK length requirements LRN and LRS are */
731  /* as follows: */
732  /* If NEQ is constant and the maximum method orders have */
733  /* their default values, then */
734  /* LRN = 20 + 16*NEQ + 3*NG, */
735  /* LRS = 22 + 9*NEQ + NEQ**2 + 3*NG (JT = 1 or 2), */
736  /* LRS = 22 + 10*NEQ + (2*ML+MU)*NEQ + 3*NG (JT = 4 or 5). */
737  /* Under any other conditions, LRN and LRS are given by: */
738  /* LRN = 20 + NYH*(MXORDN+1) + 3*NEQ + 3*NG, */
739  /* LRS = 20 + NYH*(MXORDS+1) + 3*NEQ + LMAT + 3*NG, */
740  /* where */
741  /* NYH = the initial value of NEQ, */
742  /* MXORDN = 12, unless a smaller value is given as an */
743  /* optional input, */
744  /* MXORDS = 5, unless a smaller value is given as an */
745  /* optional input, */
746  /* LMAT = length of matrix work space: */
747  /* LMAT = NEQ**2 + 2 if JT = 1 or 2, */
748  /* LMAT = (2*ML + MU + 1)*NEQ + 2 if JT = 4 or 5. */
749 
750  /* --- Dynamic Length Case --- */
751  /* If the length of RWORK is to be dynamic, then it should */
752  /* be at least LRN or LRS, as defined above, depending on the */
753  /* current method. Initially, it must be at least LRN (since */
754  /* DLSODAR starts with the nonstiff method). On any return */
755  /* from DLSODAR, the optional output MCUR indicates the current */
756  /* method. If MCUR differs from the value it had on the */
757  /* previous return, or if there has only been one call to */
758  /* DLSODAR and MCUR is now 2, then DLSODAR has switched */
759  /* methods during the last call, and the length of RWORK */
760  /* should be reset (to LRN if MCUR = 1, or to LRS if */
761  /* MCUR = 2). (An increase in the RWORK length is required */
762  /* if DLSODAR returned ISTATE = -7, but not otherwise.) */
763  /* After resetting the length, call DLSODAR with ISTATE = 3 */
764  /* to signal that change. */
765 
766  /* LRW = the length of the array RWORK, as declared by the user. */
767  /* (This will be checked by the solver.) */
768 
769  /* IWORK = an integer array for work space. */
770  /* As DLSODAR switches automatically between stiff and nonstiff */
771  /* methods, the required length of IWORK can change during */
772  /* problem, between */
773  /* LIS = 20 + NEQ and LIN = 20, */
774  /* respectively. Thus the IWORK array passed to DLSODAR can */
775  /* either have a fixed length of at least 20 + NEQ, or have a */
776  /* dynamic length of at least LIN or LIS, depending on the */
777  /* current method. The comments on dynamic length under */
778  /* RWORK above apply here. Initially, this length need */
779  /* only be at least LIN = 20. */
780 
781  /* The first few words of IWORK are used for conditional and */
782  /* optional inputs and optional outputs. */
783 
784  /* The following 2 words in IWORK are conditional inputs: */
785  /* IWORK(1) = ML These are the lower and upper */
786  /* IWORK(2) = MU half-bandwidths, respectively, of the */
787  /* banded Jacobian, excluding the main diagonal. */
788  /* The band is defined by the matrix locations */
789  /* (i,j) with i-ML .le. j .le. i+MU. ML and MU */
790  /* must satisfy 0 .le. ML,MU .le. NEQ-1. */
791  /* These are required if JT is 4 or 5, and */
792  /* ignored otherwise. ML and MU may in fact be */
793  /* the band parameters for a matrix to which */
794  /* df/dy is only approximately equal. */
795 
796  /* LIW = the length of the array IWORK, as declared by the user. */
797  /* (This will be checked by the solver.) */
798 
799  /* Note: The base addresses of the work arrays must not be */
800  /* altered between calls to DLSODAR for the same problem. */
801  /* The contents of the work arrays must not be altered */
802  /* between calls, except possibly for the conditional and */
803  /* optional inputs, and except for the last 3*NEQ words of RWORK. */
804  /* The latter space is used for internal scratch space, and so is */
805  /* available for use by the user outside DLSODAR between calls, if */
806  /* desired (but not for use by F, JAC, or G). */
807 
808  /* JAC = the name of the user-supplied routine to compute the */
809  /* Jacobian matrix, df/dy, if JT = 1 or 4. The JAC routine */
810  /* is optional, but if the problem is expected to be stiff much */
811  /* of the time, you are encouraged to supply JAC, for the sake */
812  /* of efficiency. (Alternatively, set JT = 2 or 5 to have */
813  /* DLSODAR compute df/dy internally by difference quotients.) */
814  /* If and when DLSODAR uses df/dy, it treats this NEQ by NEQ */
815  /* matrix either as full (JT = 1 or 2), or as banded (JT = */
816  /* 4 or 5) with half-bandwidths ML and MU (discussed under */
817  /* IWORK above). In either case, if JT = 1 or 4, the JAC */
818  /* routine must compute df/dy as a function of the scalar t */
819  /* and the vector y. It is to have the form */
820  /* SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD) */
821  /* DOUBLE PRECISION T, Y(*), PD(NROWPD,*) */
822  /* where NEQ, T, Y, ML, MU, and NROWPD are input and the array */
823  /* PD is to be loaded with partial derivatives (elements of */
824  /* the Jacobian matrix) on output. PD must be given a first */
825  /* dimension of NROWPD. T and Y have the same meaning as in */
826  /* Subroutine F. */
827  /* In the full matrix case (JT = 1), ML and MU are */
828  /* ignored, and the Jacobian is to be loaded into PD in */
829  /* columnwise manner, with df(i)/dy(j) loaded into pd(i,j). */
830  /* In the band matrix case (JT = 4), the elements */
831  /* within the band are to be loaded into PD in columnwise */
832  /* manner, with diagonal lines of df/dy loaded into the rows */
833  /* of PD. Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j). */
834  /* ML and MU are the half-bandwidth parameters (see IWORK). */
835  /* The locations in PD in the two triangular areas which */
836  /* correspond to nonexistent matrix elements can be ignored */
837  /* or loaded arbitrarily, as they are overwritten by DLSODAR. */
838  /* JAC need not provide df/dy exactly. A crude */
839  /* approximation (possibly with a smaller bandwidth) will do. */
840  /* In either case, PD is preset to zero by the solver, */
841  /* so that only the nonzero elements need be loaded by JAC. */
842  /* Each call to JAC is preceded by a call to F with the same */
843  /* arguments NEQ, T, and Y. Thus to gain some efficiency, */
844  /* intermediate quantities shared by both calculations may be */
845  /* saved in a user Common block by F and not recomputed by JAC, */
846  /* if desired. Also, JAC may alter the Y array, if desired. */
847  /* JAC must be declared External in the calling program. */
848  /* Subroutine JAC may access user-defined quantities in */
849  /* NEQ(2),... and/or in Y(NEQ(1)+1),... if NEQ is an array */
850  /* (dimensioned in JAC) and/or Y has length exceeding NEQ(1). */
851  /* See the descriptions of NEQ and Y above. */
852 
853  /* JT = Jacobian type indicator. Used only for input. */
854  /* JT specifies how the Jacobian matrix df/dy will be */
855  /* treated, if and when DLSODAR requires this matrix. */
856  /* JT has the following values and meanings: */
857  /* 1 means a user-supplied full (NEQ by NEQ) Jacobian. */
858  /* 2 means an internally generated (difference quotient) full */
859  /* Jacobian (using NEQ extra calls to F per df/dy value). */
860  /* 4 means a user-supplied banded Jacobian. */
861  /* 5 means an internally generated banded Jacobian (using */
862  /* ML+MU+1 extra calls to F per df/dy evaluation). */
863  /* If JT = 1 or 4, the user must supply a Subroutine JAC */
864  /* (the name is arbitrary) as described above under JAC. */
865  /* If JT = 2 or 5, a dummy argument can be used. */
866 
867  /* G = the name of subroutine for constraint functions, whose */
868  /* roots are desired during the integration. It is to have */
869  /* the form */
870  /* SUBROUTINE G (NEQ, T, Y, NG, GOUT) */
871  /* DOUBLE PRECISION T, Y(*), GOUT(NG) */
872  /* where NEQ, T, Y, and NG are input, and the array GOUT */
873  /* is output. NEQ, T, and Y have the same meaning as in */
874  /* the F routine, and GOUT is an array of length NG. */
875  /* For i = 1,...,NG, this routine is to load into GOUT(i) */
876  /* the value at (T,Y) of the i-th constraint function g(i). */
877  /* DLSODAR will find roots of the g(i) of odd multiplicity */
878  /* (i.e. sign changes) as they occur during the integration. */
879  /* G must be declared External in the calling program. */
880 
881  /* Caution: Because of numerical errors in the functions */
882  /* g(i) due to roundoff and integration error, DLSODAR may */
883  /* return false roots, or return the same root at two or more */
884  /* nearly equal values of t. If such false roots are */
885  /* suspected, the user should consider smaller error tolerances */
886  /* and/or higher precision in the evaluation of the g(i). */
887 
888  /* If a root of some g(i) defines the end of the problem, */
889  /* the input to DLSODAR should nevertheless allow integration */
890  /* to a point slightly past that root, so that DLSODAR can */
891  /* locate the root by interpolation. */
892 
893  /* Subroutine G may access user-defined quantities in */
894  /* NEQ(2),... and Y(NEQ(1)+1),... if NEQ is an array */
895  /* (dimensioned in G) and/or Y has length exceeding NEQ(1). */
896  /* See the descriptions of NEQ and Y above. */
897 
898  /* NG = number of constraint functions g(i). If there are none, */
899  /* set NG = 0, and pass a dummy name for G. */
900 
901  /* JROOT = integer array of length NG. Used only for output. */
902  /* On a return with ISTATE = 3 (one or more roots found), */
903  /* JROOT(i) = 1 if g(i) has a root at T, or JROOT(i) = 0 if not. */
904  /* ----------------------------------------------------------------------- */
905  /* Optional Inputs. */
906 
907  /* The following is a list of the optional inputs provided for in the */
908  /* call sequence. (See also Part 2.) For each such input variable, */
909  /* this table lists its name as used in this documentation, its */
910  /* location in the call sequence, its meaning, and the default value. */
911  /* The use of any of these inputs requires IOPT = 1, and in that */
912  /* case all of these inputs are examined. A value of zero for any */
913  /* of these optional inputs will cause the default value to be used. */
914  /* Thus to use a subset of the optional inputs, simply preload */
915  /* locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, and */
916  /* then set those of interest to nonzero values. */
917 
918  /* Name Location Meaning and Default Value */
919 
920  /* H0 RWORK(5) the step size to be attempted on the first step. */
921  /* The default value is determined by the solver. */
922 
923  /* HMAX RWORK(6) the maximum absolute step size allowed. */
924  /* The default value is infinite. */
925 
926  /* HMIN RWORK(7) the minimum absolute step size allowed. */
927  /* The default value is 0. (This lower bound is not */
928  /* enforced on the final step before reaching TCRIT */
929  /* when ITASK = 4 or 5.) */
930 
931  /* IXPR IWORK(5) flag to generate extra printing at method switches. */
932  /* IXPR = 0 means no extra printing (the default). */
933  /* IXPR = 1 means print data on each switch. */
934  /* T, H, and NST will be printed on the same logical */
935  /* unit as used for error messages. */
936 
937  /* MXSTEP IWORK(6) maximum number of (internally defined) steps */
938  /* allowed during one call to the solver. */
939  /* The default value is 500. */
940 
941  /* MXHNIL IWORK(7) maximum number of messages printed (per problem) */
942  /* warning that T + H = T on a step (H = step size). */
943  /* This must be positive to result in a non-default */
944  /* value. The default value is 10. */
945 
946  /* MXORDN IWORK(8) the maximum order to be allowed for the nonstiff */
947  /* (Adams) method. The default value is 12. */
948  /* If MXORDN exceeds the default value, it will */
949  /* be reduced to the default value. */
950  /* MXORDN is held constant during the problem. */
951 
952  /* MXORDS IWORK(9) the maximum order to be allowed for the stiff */
953  /* (BDF) method. The default value is 5. */
954  /* If MXORDS exceeds the default value, it will */
955  /* be reduced to the default value. */
956  /* MXORDS is held constant during the problem. */
957  /* ----------------------------------------------------------------------- */
958  /* Optional Outputs. */
959 
960  /* As optional additional output from DLSODAR, the variables listed */
961  /* below are quantities related to the performance of DLSODAR */
962  /* which are available to the user. These are communicated by way of */
963  /* the work arrays, but also have internal mnemonic names as shown. */
964  /* Except where stated otherwise, all of these outputs are defined */
965  /* on any successful return from DLSODAR, and on any return with */
966  /* ISTATE = -1, -2, -4, -5, or -6. On an illegal input return */
967  /* (ISTATE = -3), they will be unchanged from their existing values */
968  /* (if any), except possibly for TOLSF, LENRW, and LENIW. */
969  /* On any error return, outputs relevant to the error will be defined, */
970  /* as noted below. */
971 
972  /* Name Location Meaning */
973 
974  /* HU RWORK(11) the step size in t last used (successfully). */
975 
976  /* HCUR RWORK(12) the step size to be attempted on the next step. */
977 
978  /* TCUR RWORK(13) the current value of the independent variable */
979  /* which the solver has actually reached, i.e. the */
980  /* current internal mesh point in t. On output, TCUR */
981  /* will always be at least as far as the argument */
982  /* T, but may be farther (if interpolation was done). */
983 
984  /* TOLSF RWORK(14) a tolerance scale factor, greater than 1.0, */
985  /* computed when a request for too much accuracy was */
986  /* detected (ISTATE = -3 if detected at the start of */
987  /* the problem, ISTATE = -2 otherwise). If ITOL is */
988  /* left unaltered but RTOL and ATOL are uniformly */
989  /* scaled up by a factor of TOLSF for the next call, */
990  /* then the solver is deemed likely to succeed. */
991  /* (The user may also ignore TOLSF and alter the */
992  /* tolerance parameters in any other way appropriate.) */
993 
994  /* TSW RWORK(15) the value of t at the time of the last method */
995  /* switch, if any. */
996 
997  /* NGE IWORK(10) the number of g evaluations for the problem so far. */
998 
999  /* NST IWORK(11) the number of steps taken for the problem so far. */
1000 
1001  /* NFE IWORK(12) the number of f evaluations for the problem so far. */
1002 
1003  /* NJE IWORK(13) the number of Jacobian evaluations (and of matrix */
1004  /* LU decompositions) for the problem so far. */
1005 
1006  /* NQU IWORK(14) the method order last used (successfully). */
1007 
1008  /* NQCUR IWORK(15) the order to be attempted on the next step. */
1009 
1010  /* IMXER IWORK(16) the index of the component of largest magnitude in */
1011  /* the weighted local error vector (E(i)/EWT(i)), */
1012  /* on an error return with ISTATE = -4 or -5. */
1013 
1014  /* LENRW IWORK(17) the length of RWORK actually required, assuming */
1015  /* that the length of RWORK is to be fixed for the */
1016  /* rest of the problem, and that switching may occur. */
1017  /* This is defined on normal returns and on an illegal */
1018  /* input return for insufficient storage. */
1019 
1020  /* LENIW IWORK(18) the length of IWORK actually required, assuming */
1021  /* that the length of IWORK is to be fixed for the */
1022  /* rest of the problem, and that switching may occur. */
1023  /* This is defined on normal returns and on an illegal */
1024  /* input return for insufficient storage. */
1025 
1026  /* MUSED IWORK(19) the method indicator for the last successful step: */
1027  /* 1 means Adams (nonstiff), 2 means BDF (stiff). */
1028 
1029  /* MCUR IWORK(20) the current method indicator: */
1030  /* 1 means Adams (nonstiff), 2 means BDF (stiff). */
1031  /* This is the method to be attempted */
1032  /* on the next step. Thus it differs from MUSED */
1033  /* only if a method switch has just been made. */
1034 
1035  /* The following two arrays are segments of the RWORK array which */
1036  /* may also be of interest to the user as optional outputs. */
1037  /* For each array, the table below gives its internal name, */
1038  /* its base address in RWORK, and its description. */
1039 
1040  /* Name Base Address Description */
1041 
1042  /* YH 21 + 3*NG the Nordsieck history array, of size NYH by */
1043  /* (NQCUR + 1), where NYH is the initial value */
1044  /* of NEQ. For j = 0,1,...,NQCUR, column j+1 */
1045  /* of YH contains HCUR**j/factorial(j) times */
1046  /* the j-th derivative of the interpolating */
1047  /* polynomial currently representing the solution, */
1048  /* evaluated at t = TCUR. */
1049 
1050  /* ACOR LACOR array of size NEQ used for the accumulated */
1051  /* (from Common corrections on each step, scaled on output */
1052  /* as noted) to represent the estimated local error in y */
1053  /* on the last step. This is the vector E in */
1054  /* the description of the error control. It is */
1055  /* defined only on a successful return from */
1056  /* DLSODAR. The base address LACOR is obtained by */
1057  /* including in the user's program the */
1058  /* following 2 lines: */
1059  /* COMMON /DLS001/ RLS(218), ILS(37) */
1060  /* LACOR = ILS(22) */
1061 
1062  /* ----------------------------------------------------------------------- */
1063  /* Part 2. Other Routines Callable. */
1064 
1065  /* The following are optional calls which the user may make to */
1066  /* gain additional capabilities in conjunction with DLSODAR. */
1067  /* (The routines XSETUN and XSETF are designed to conform to the */
1068  /* SLATEC error handling package.) */
1069 
1070  /* Form of Call Function */
1071  /* CALL XSETUN(LUN) Set the logical unit number, LUN, for */
1072  /* output of messages from DLSODAR, if */
1073  /* the default is not desired. */
1074  /* The default value of LUN is 6. */
1075 
1076  /* CALL XSETF(MFLAG) Set a flag to control the printing of */
1077  /* messages by DLSODAR. */
1078  /* MFLAG = 0 means do not print. (Danger: */
1079  /* This risks losing valuable information.) */
1080  /* MFLAG = 1 means print (the default). */
1081 
1082  /* Either of the above calls may be made at */
1083  /* any time and will take effect immediately. */
1084 
1085  /* CALL DSRCAR(RSAV,ISAV,JOB) saves and restores the contents of */
1086  /* the internal Common blocks used by */
1087  /* DLSODAR (see Part 3 below). */
1088  /* RSAV must be a real array of length 245 */
1089  /* or more, and ISAV must be an integer */
1090  /* array of length 55 or more. */
1091  /* JOB=1 means save Common into RSAV/ISAV. */
1092  /* JOB=2 means restore Common from RSAV/ISAV. */
1093  /* DSRCAR is useful if one is */
1094  /* interrupting a run and restarting */
1095  /* later, or alternating between two or */
1096  /* more problems solved with DLSODAR. */
1097 
1098  /* CALL DINTDY(,,,,,) Provide derivatives of y, of various */
1099  /* (see below) orders, at a specified point t, if */
1100  /* desired. It may be called only after */
1101  /* a successful return from DLSODAR. */
1102 
1103  /* The detailed instructions for using DINTDY are as follows. */
1104  /* The form of the call is: */
1105 
1106  /* LYH = 21 + 3*NG */
1107  /* CALL DINTDY (T, K, RWORK(LYH), NYH, DKY, IFLAG) */
1108 
1109  /* The input parameters are: */
1110 
1111  /* T = value of independent variable where answers are desired */
1112  /* (normally the same as the T last returned by DLSODAR). */
1113  /* For valid results, T must lie between TCUR - HU and TCUR. */
1114  /* (See optional outputs for TCUR and HU.) */
1115  /* K = integer order of the derivative desired. K must satisfy */
1116  /* 0 .le. K .le. NQCUR, where NQCUR is the current order */
1117  /* (see optional outputs). The capability corresponding */
1118  /* to K = 0, i.e. computing y(t), is already provided */
1119  /* by DLSODAR directly. Since NQCUR .ge. 1, the first */
1120  /* derivative dy/dt is always available with DINTDY. */
1121  /* LYH = 21 + 3*NG = base address in RWORK of the history array YH. */
1122  /* NYH = column length of YH, equal to the initial value of NEQ. */
1123 
1124  /* The output parameters are: */
1125 
1126  /* DKY = a real array of length NEQ containing the computed value */
1127  /* of the K-th derivative of y(t). */
1128  /* IFLAG = integer flag, returned as 0 if K and T were legal, */
1129  /* -1 if K was illegal, and -2 if T was illegal. */
1130  /* On an error return, a message is also written. */
1131  /* ----------------------------------------------------------------------- */
1132  /* Part 3. Common Blocks. */
1133 
1134  /* If DLSODAR is to be used in an overlay situation, the user */
1135  /* must declare, in the primary overlay, the variables in: */
1136  /* (1) the call sequence to DLSODAR, and */
1137  /* (2) the three internal Common blocks */
1138  /* /DLS001/ of length 255 (218 double precision words */
1139  /* followed by 37 integer words), */
1140  /* /DLSA01/ of length 31 (22 double precision words */
1141  /* followed by 9 integer words). */
1142  /* /DLSR01/ of length 7 (3 double precision words */
1143  /* followed by 4 integer words). */
1144 
1145  /* If DLSODAR is used on a system in which the contents of internal */
1146  /* Common blocks are not preserved between calls, the user should */
1147  /* declare the above Common blocks in the calling program to insure */
1148  /* that their contents are preserved. */
1149 
1150  /* If the solution of a given problem by DLSODAR is to be interrupted */
1151  /* and then later continued, such as when restarting an interrupted run */
1152  /* or alternating between two or more problems, the user should save, */
1153  /* following the return from the last DLSODAR call prior to the */
1154  /* interruption, the contents of the call sequence variables and the */
1155  /* internal Common blocks, and later restore these values before the */
1156  /* next DLSODAR call for that problem. To save and restore the Common */
1157  /* blocks, use Subroutine DSRCAR (see Part 2 above). */
1158 
1159  /* ----------------------------------------------------------------------- */
1160  /* Part 4. Optionally Replaceable Solver Routines. */
1161 
1162  /* Below is a description of a routine in the DLSODAR package which */
1163  /* relates to the measurement of errors, and can be */
1164  /* replaced by a user-supplied version, if desired. However, since such */
1165  /* a replacement may have a major impact on performance, it should be */
1166  /* done only when absolutely necessary, and only with great caution. */
1167  /* (Note: The means by which the package version of a routine is */
1168  /* superseded by the user's version may be system-dependent.) */
1169 
1170  /* (a) DEWSET. */
1171  /* The following subroutine is called just before each internal */
1172  /* integration step, and sets the array of error weights, EWT, as */
1173  /* described under ITOL/RTOL/ATOL above: */
1174  /* Subroutine DEWSET (NEQ, ITOL, RTOL, ATOL, YCUR, EWT) */
1175  /* where NEQ, ITOL, RTOL, and ATOL are as in the DLSODAR call sequence, */
1176  /* YCUR contains the current dependent variable vector, and */
1177  /* EWT is the array of weights set by DEWSET. */
1178 
1179  /* If the user supplies this subroutine, it must return in EWT(i) */
1180  /* (i = 1,...,NEQ) a positive quantity suitable for comparing errors */
1181  /* in y(i) to. The EWT array returned by DEWSET is passed to the */
1182  /* DMNORM routine, and also used by DLSODAR in the computation */
1183  /* of the optional output IMXER, and the increments for difference */
1184  /* quotient Jacobians. */
1185 
1186  /* In the user-supplied version of DEWSET, it may be desirable to use */
1187  /* the current values of derivatives of y. Derivatives up to order NQ */
1188  /* are available from the history array YH, described above under */
1189  /* optional outputs. In DEWSET, YH is identical to the YCUR array, */
1190  /* extended to NQ + 1 columns with a column length of NYH and scale */
1191  /* factors of H**j/factorial(j). On the first call for the problem, */
1192  /* given by NST = 0, NQ is 1 and H is temporarily set to 1.0. */
1193  /* NYH is the initial value of NEQ. The quantities NQ, H, and NST */
1194  /* can be obtained by including in DEWSET the statements: */
1195  /* DOUBLE PRECISION RLS */
1196  /* COMMON /DLS001/ RLS(218),ILS(37) */
1197  /* NQ = ILS(33) */
1198  /* NST = ILS(34) */
1199  /* H = RLS(212) */
1200  /* Thus, for example, the current value of dy/dt can be obtained as */
1201  /* YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is */
1202  /* unnecessary when NST = 0). */
1203  /* ----------------------------------------------------------------------- */
1204 
1205  /* ***REVISION HISTORY (YYYYMMDD) */
1206  /* 19811102 DATE WRITTEN */
1207  /* 19820126 Fixed bug in tests of work space lengths; */
1208  /* minor corrections in main prologue and comments. */
1209  /* 19820507 Fixed bug in RCHEK in setting HMING. */
1210  /* 19870330 Major update: corrected comments throughout; */
1211  /* removed TRET from Common; rewrote EWSET with 4 loops; */
1212  /* fixed t test in INTDY; added Cray directives in STODA; */
1213  /* in STODA, fixed DELP init. and logic around PJAC call; */
1214  /* combined routines to save/restore Common; */
1215  /* passed LEVEL = 0 in error message calls (except run abort). */
1216  /* 19970225 Fixed lines setting JSTART = -2 in Subroutine LSODAR. */
1217  /* 20010425 Major update: convert source lines to upper case; */
1218  /* added *DECK lines; changed from 1 to * in dummy dimensions; */
1219  /* changed names R1MACH/D1MACH to RUMACH/DUMACH; */
1220  /* renamed routines for uniqueness across single/double prec.; */
1221  /* converted intrinsic names to generic form; */
1222  /* removed ILLIN and NTREP (data loaded) from Common; */
1223  /* removed all 'own' variables from Common; */
1224  /* changed error messages to quoted strings; */
1225  /* replaced XERRWV/XERRWD with 1993 revised version; */
1226  /* converted prologues, comments, error messages to mixed case; */
1227  /* numerous corrections to prologues and internal comments. */
1228  /* 20010507 Converted single precision source to double precision. */
1229  /* 20010613 Revised excess accuracy test (to match rest of ODEPACK). */
1230  /* 20010808 Fixed bug in DPRJA (matrix in DBNORM call). */
1231  /* 20020502 Corrected declarations in descriptions of user routines. */
1232  /* 20031105 Restored 'own' variables to Common blocks, to enable */
1233  /* interrupt/restart feature. */
1234  /* 20031112 Added SAVE statements for data-loaded constants. */
1235 
1236  /* ----------------------------------------------------------------------- */
1237  /* Other routines in the DLSODAR package. */
1238 
1239  /* In addition to Subroutine DLSODAR, the DLSODAR package includes the */
1240  /* following subroutines and function routines: */
1241  /* DRCHEK does preliminary checking for roots, and serves as an */
1242  /* interface between Subroutine DLSODAR and Subroutine DROOTS. */
1243  /* DROOTS finds the leftmost root of a set of functions. */
1244  /* DINTDY computes an interpolated value of the y vector at t = TOUT. */
1245  /* DSTODA is the core integrator, which does one step of the */
1246  /* integration and the associated error control. */
1247  /* DCFODE sets all method coefficients and test constants. */
1248  /* DPRJA computes and preprocesses the Jacobian matrix J = df/dy */
1249  /* and the Newton iteration matrix P = I - h*l0*J. */
1250  /* DSOLSY manages solution of linear system in chord iteration. */
1251  /* DEWSET sets the error weight vector EWT before each step. */
1252  /* DMNORM computes the weighted max-norm of a vector. */
1253  /* DFNORM computes the norm of a full matrix consistent with the */
1254  /* weighted max-norm on vectors. */
1255  /* DBNORM computes the norm of a band matrix consistent with the */
1256  /* weighted max-norm on vectors. */
1257  /* DSRCAR is a user-callable routine to save and restore */
1258  /* the contents of the internal Common blocks. */
1259  /* DGEFA and DGESL are routines from LINPACK for solving full */
1260  /* systems of linear algebraic equations. */
1261  /* DGBFA and DGBSL are routines from LINPACK for solving banded */
1262  /* linear systems. */
1263  /* DCOPY is one of the basic linear algebra modules (BLAS). */
1264  /* DUMACH computes the unit roundoff in a machine-independent manner. */
1265  /* XERRWD, XSETUN, XSETF, IXSAV, and IUMACH handle the printing of all */
1266  /* error messages and warnings. XERRWD is machine-dependent. */
1267  /* Note: DMNORM, DFNORM, DBNORM, DUMACH, IXSAV, and IUMACH are */
1268  /* function routines. All the others are subroutines. */
1269 
1270  /* ----------------------------------------------------------------------- */
1271  /* ----------------------------------------------------------------------- */
1272  /* The following three internal Common blocks contain */
1273  /* (a) variables which are local to any subroutine but whose values must */
1274  /* be preserved between calls to the routine ("own" variables), and */
1275  /* (b) variables which are communicated between subroutines. */
1276  /* The block DLS001 is declared in subroutines DLSODAR, DINTDY, DSTODA, */
1277  /* DPRJA, and DSOLSY. */
1278  /* The block DLSA01 is declared in subroutines DLSODAR, DSTODA, DPRJA. */
1279  /* The block DLSR01 is declared in subroutines DLSODAR, DRCHEK, DROOTS. */
1280  /* Groups of variables are replaced by dummy arrays in the Common */
1281  /* declarations in routines where those variables are not used. */
1282  /* ----------------------------------------------------------------------- */
1283 
1284  /* Parameter adjustments */
1285  --neq;
1286  --y;
1287  --rtol;
1288  --atol;
1289  --rwork;
1290  --iwork;
1291  --jroot;
1292 
1293  /* Function Body */
1294 
1295  /* ----------------------------------------------------------------------- */
1296 
1297  /* Block A. */
1298 
1299  /* This code block is executed on every call. */
1300 
1301  /* It tests ISTATE and ITASK for legality and branches appropriately. */
1302 
1303  /* If ISTATE .gt. 1 but the flag INIT shows that initialization has */
1304 
1305  /* not yet been done, an error return occurs. */
1306 
1307  /* If ISTATE = 1 and TOUT = T, return immediately. */
1308 
1309  /* ----------------------------------------------------------------------- */
1310  if (*istate < 1 || *istate > 3)
1311  {
1312  goto L601;
1313  }
1314 
1315  if (*itask < 1 || *itask > 5)
1316  {
1317  goto L602;
1318  }
1319 
1320  dlsr01_1.itaskc = *itask;
1321 
1322  if (*istate == 1)
1323  {
1324  goto L10;
1325  }
1326 
1327  if (dls001_1.init == 0)
1328  {
1329  goto L603;
1330  }
1331 
1332  if (*istate == 2)
1333  {
1334  goto L200;
1335  }
1336 
1337  goto L20;
1338 L10:
1339  dls001_1.init = 0;
1340 
1341  if (*tout == *t)
1342  {
1343  return 0;
1344  }
1345 
1346  /* ----------------------------------------------------------------------- */
1347  /* Block B. */
1348  /* The next code block is executed for the initial call (ISTATE = 1), */
1349  /* or for a continuation call with parameter changes (ISTATE = 3). */
1350  /* It contains checking of all inputs and various initializations. */
1351 
1352  /* First check legality of the non-optional inputs NEQ, ITOL, IOPT, */
1353  /* JT, ML, MU, and NG. */
1354  /* ----------------------------------------------------------------------- */
1355 L20:
1356 
1357  if (neq[1] <= 0)
1358  {
1359  goto L604;
1360  }
1361 
1362  if (*istate == 1)
1363  {
1364  goto L25;
1365  }
1366 
1367  if (neq[1] > dls001_1.n)
1368  {
1369  goto L605;
1370  }
1371 
1372 L25:
1373  dls001_1.n = neq[1];
1374 
1375  if (*itol < 1 || *itol > 4)
1376  {
1377  goto L606;
1378  }
1379 
1380  if (*iopt < 0 || *iopt > 1)
1381  {
1382  goto L607;
1383  }
1384 
1385  if (*jt == 3 || *jt < 1 || *jt > 5)
1386  {
1387  goto L608;
1388  }
1389 
1390  dlsa01_1.jtyp = *jt;
1391 
1392  if (*jt <= 2)
1393  {
1394  goto L30;
1395  }
1396 
1397  ml = iwork[1];
1398  mu = iwork[2];
1399 
1400  if (ml < 0 || ml >= dls001_1.n)
1401  {
1402  goto L609;
1403  }
1404 
1405  if (mu < 0 || mu >= dls001_1.n)
1406  {
1407  goto L610;
1408  }
1409 
1410 L30:
1411 
1412  if (*ng < 0)
1413  {
1414  goto L630;
1415  }
1416 
1417  if (*istate == 1)
1418  {
1419  goto L35;
1420  }
1421 
1422  if (dlsr01_1.irfnd == 0 && *ng != dlsr01_1.ngc)
1423  {
1424  goto L631;
1425  }
1426 
1427 L35:
1428  dlsr01_1.ngc = *ng;
1429 
1430  /* Next process and check the optional inputs. -------------------------- */
1431  if (*iopt == 1)
1432  {
1433  goto L40;
1434  }
1435 
1436  dlsa01_1.ixpr = 0;
1437  dls001_1.mxstep = mxstp0;
1438  dls001_1.mxhnil = mxhnl0;
1439  dls001_1.hmxi = 0.;
1440  dls001_1.hmin = 0.;
1441 
1442  if (*istate != 1)
1443  {
1444  goto L60;
1445  }
1446 
1447  h0 = 0.;
1448  dlsa01_1.mxordn = mord[0];
1449  dlsa01_1.mxords = mord[1];
1450  goto L60;
1451 L40:
1452  dlsa01_1.ixpr = iwork[5];
1453 
1454  if (dlsa01_1.ixpr < 0 || dlsa01_1.ixpr > 1)
1455  {
1456  goto L611;
1457  }
1458 
1459  dls001_1.mxstep = iwork[6];
1460 
1461  if (dls001_1.mxstep < 0)
1462  {
1463  goto L612;
1464  }
1465 
1466  if (dls001_1.mxstep == 0)
1467  {
1468  dls001_1.mxstep = mxstp0;
1469  }
1470 
1471  dls001_1.mxhnil = iwork[7];
1472 
1473  if (dls001_1.mxhnil < 0)
1474  {
1475  goto L613;
1476  }
1477 
1478  if (dls001_1.mxhnil == 0)
1479  {
1480  dls001_1.mxhnil = mxhnl0;
1481  }
1482 
1483  if (*istate != 1)
1484  {
1485  goto L50;
1486  }
1487 
1488  h0 = rwork[5];
1489  dlsa01_1.mxordn = iwork[8];
1490 
1491  if (dlsa01_1.mxordn < 0)
1492  {
1493  goto L628;
1494  }
1495 
1496  if (dlsa01_1.mxordn == 0)
1497  {
1498  dlsa01_1.mxordn = 100;
1499  }
1500 
1501  dlsa01_1.mxordn = std::min(dlsa01_1.mxordn, mord[0]);
1502  dlsa01_1.mxords = iwork[9];
1503 
1504  if (dlsa01_1.mxords < 0)
1505  {
1506  goto L629;
1507  }
1508 
1509  if (dlsa01_1.mxords == 0)
1510  {
1511  dlsa01_1.mxords = 100;
1512  }
1513 
1514  dlsa01_1.mxords = std::min(dlsa01_1.mxords, mord[1]);
1515 
1516  if ((*tout - *t) * h0 < 0.)
1517  {
1518  goto L614;
1519  }
1520 
1521 L50:
1522  hmax = rwork[6];
1523 
1524  if (hmax < 0.)
1525  {
1526  goto L615;
1527  }
1528 
1529  dls001_1.hmxi = 0.;
1530 
1531  if (hmax > 0.)
1532  {
1533  dls001_1.hmxi = 1. / hmax;
1534  }
1535 
1536  dls001_1.hmin = rwork[7];
1537 
1538  if (dls001_1.hmin < 0.)
1539  {
1540  goto L616;
1541  }
1542 
1543  /* ----------------------------------------------------------------------- */
1544  /* Set work array pointers and check lengths LRW and LIW. */
1545  /* If ISTATE = 1, METH is initialized to 1 here to facilitate the */
1546  /* checking of work space lengths. */
1547  /* Pointers to segments of RWORK and IWORK are named by prefixing L to */
1548  /* the name of the segment. E.g., the segment YH starts at RWORK(LYH). */
1549  /* Segments of RWORK (in order) are denoted G0, G1, GX, YH, WM, */
1550  /* EWT, SAVF, ACOR. */
1551  /* If the lengths provided are insufficient for the current method, */
1552  /* an error return occurs. This is treated as illegal input on the */
1553  /* first call, but as a problem interruption with ISTATE = -7 on a */
1554  /* continuation call. If the lengths are sufficient for the current */
1555  /* method but not for both methods, a warning message is sent. */
1556  /* ----------------------------------------------------------------------- */
1557 L60:
1558 
1559  if (*istate == 1)
1560  {
1561  dls001_1.meth = 1;
1562  }
1563 
1564  if (*istate == 1)
1565  {
1566  dls001_1.nyh = dls001_1.n;
1567  }
1568 
1569  dlsr01_1.lg0 = 21;
1570  dlsr01_1.lg1 = dlsr01_1.lg0 + *ng;
1571  dlsr01_1.lgx = dlsr01_1.lg1 + *ng;
1572  lyhnew = dlsr01_1.lgx + *ng;
1573 
1574  if (*istate == 1)
1575  {
1576  dls001_1.lyh = lyhnew;
1577  }
1578 
1579  if (lyhnew == dls001_1.lyh)
1580  {
1581  goto L62;
1582  }
1583 
1584  /* If ISTATE = 3 and NG was changed, shift YH to its new location. ------ */
1585  lenyh = dls001_1.l * dls001_1.nyh;
1586 
1587  if (*lrw < lyhnew - 1 + lenyh)
1588  {
1589  goto L62;
1590  }
1591 
1592  i1 = 1;
1593 
1594  if (lyhnew > dls001_1.lyh)
1595  {
1596  i1 = -1;
1597  }
1598 
1599  dcopy_(&lenyh, &rwork[dls001_1.lyh], &i1, &rwork[lyhnew], &i1);
1600  dls001_1.lyh = lyhnew;
1601 L62:
1602  len1n = lyhnew - 1 + (dlsa01_1.mxordn + 1) * dls001_1.nyh;
1603  len1s = lyhnew - 1 + (dlsa01_1.mxords + 1) * dls001_1.nyh;
1604  dls001_1.lwm = len1s + 1;
1605 
1606  if (*jt <= 2)
1607  {
1608  lenwm = dls001_1.n * dls001_1.n + 2;
1609  }
1610 
1611  if (*jt >= 4)
1612  {
1613  lenwm = ((ml << 1) + mu + 1) * dls001_1.n + 2;
1614  }
1615 
1616  len1s += lenwm;
1617  len1c = len1n;
1618 
1619  if (dls001_1.meth == 2)
1620  {
1621  len1c = len1s;
1622  }
1623 
1624  len1 = std::max(len1n, len1s);
1625  len2 = dls001_1.n * 3;
1626  lenrw = len1 + len2;
1627  lenrwc = len1c + len2;
1628  iwork[17] = lenrw;
1629  dls001_1.liwm = 1;
1630  leniw = dls001_1.n + 20;
1631  leniwc = 20;
1632 
1633  if (dls001_1.meth == 2)
1634  {
1635  leniwc = leniw;
1636  }
1637 
1638  iwork[18] = leniw;
1639 
1640  if (*istate == 1 && *lrw < lenrwc)
1641  {
1642  goto L617;
1643  }
1644 
1645  if (*istate == 1 && *liw < leniwc)
1646  {
1647  goto L618;
1648  }
1649 
1650  if (*istate == 3 && *lrw < lenrwc)
1651  {
1652  goto L550;
1653  }
1654 
1655  if (*istate == 3 && *liw < leniwc)
1656  {
1657  goto L555;
1658  }
1659 
1660  dls001_1.lewt = len1 + 1;
1661  dlsa01_1.insufr = 0;
1662 
1663  if (*lrw >= lenrw)
1664  {
1665  goto L65;
1666  }
1667 
1668  dlsa01_1.insufr = 2;
1669  dls001_1.lewt = len1c + 1;
1670  msg = "DLSODAR- Warning.. RWORK length is sufficient for now, but ";
1671  mxerrwd(msg, &c__60, &c__103, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
1672  c_b76, (C_INT)60);
1673  msg = " may not be later. Integration will proceed anyway. ";
1674  mxerrwd(msg, &c__60, &c__103, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
1675  c_b76, (C_INT)60);
1676  msg = " Length needed is LENRW = I1, while LRW = I2.";
1677  mxerrwd(msg, &c__50, &c__103, &c__0, &c__2, &lenrw, lrw, &c__0, &c_b76, &
1678  c_b76, (C_INT)60);
1679 L65:
1680  dls001_1.lsavf = dls001_1.lewt + dls001_1.n;
1681  dls001_1.lacor = dls001_1.lsavf + dls001_1.n;
1682  dlsa01_1.insufi = 0;
1683 
1684  if (*liw >= leniw)
1685  {
1686  goto L70;
1687  }
1688 
1689  dlsa01_1.insufi = 2;
1690  msg = "DLSODAR- Warning.. IWORK length is sufficient for now, but ";
1691  mxerrwd(msg, &c__60, &c__104, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
1692  c_b76, (C_INT)60);
1693  msg = " may not be later. Integration will proceed anyway. ";
1694  mxerrwd(msg, &c__60, &c__104, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
1695  c_b76, (C_INT)60);
1696  msg = " Length needed is LENIW = I1, while LIW = I2.";
1697  mxerrwd(msg, &c__50, &c__104, &c__0, &c__2, &leniw, liw, &c__0, &c_b76, &
1698  c_b76, (C_INT)60);
1699 L70:
1700  /* Check RTOL and ATOL for legality. ------------------------------------ */
1701  rtoli = rtol[1];
1702  atoli = atol[1];
1703  i__1 = dls001_1.n;
1704 
1705  for (i__ = 1; i__ <= i__1; ++i__)
1706  {
1707  if (*itol >= 3)
1708  {
1709  rtoli = rtol[i__];
1710  }
1711 
1712  if (*itol == 2 || *itol == 4)
1713  {
1714  atoli = atol[i__];
1715  }
1716 
1717  if (rtoli < 0.)
1718  {
1719  goto L619;
1720  }
1721 
1722  if (atoli < 0.)
1723  {
1724  goto L620;
1725  }
1726 
1727  /* L75: */
1728  }
1729 
1730  if (*istate == 1)
1731  {
1732  goto L100;
1733  }
1734 
1735  /* if ISTATE = 3, set flag to signal parameter changes to DSTODA. ------- */
1736  dls001_1.jstart = -1;
1737 
1738  if (dls001_1.n == dls001_1.nyh)
1739  {
1740  goto L200;
1741  }
1742 
1743  /* NEQ was reduced. zero part of yh to avoid undefined references. ----- */
1744  i1 = dls001_1.lyh + dls001_1.l * dls001_1.nyh;
1745  i2 = dls001_1.lyh + (dls001_1.maxord + 1) * dls001_1.nyh - 1;
1746 
1747  if (i1 > i2)
1748  {
1749  goto L200;
1750  }
1751 
1752  i__1 = i2;
1753 
1754  for (i__ = i1; i__ <= i__1; ++i__)
1755  {
1756  /* L95: */
1757  rwork[i__] = 0.;
1758  }
1759 
1760  goto L200;
1761  /* ----------------------------------------------------------------------- */
1762  /* Block C. */
1763  /* The next block is for the initial call only (ISTATE = 1). */
1764  /* It contains all remaining initializations, the initial call to F, */
1765  /* and the calculation of the initial step size. */
1766  /* The error weights in EWT are inverted after being loaded. */
1767  /* ----------------------------------------------------------------------- */
1768 L100:
1769  dls001_1.uround = std::numeric_limits< C_FLOAT64 >::epsilon();
1770  dls001_1.tn = *t;
1771  dlsa01_1.tsw = *t;
1772  dls001_1.maxord = dlsa01_1.mxordn;
1773 
1774  if (*itask != 4 && *itask != 5)
1775  {
1776  goto L110;
1777  }
1778 
1779  tcrit = rwork[1];
1780 
1781  if ((tcrit - *tout) * (*tout - *t) < 0.)
1782  {
1783  goto L625;
1784  }
1785 
1786  if (h0 != 0. && (*t + h0 - tcrit) * h0 > 0.)
1787  {
1788  h0 = tcrit - *t;
1789  }
1790 
1791 L110:
1792  dls001_1.jstart = 0;
1793  dls001_1.nhnil = 0;
1794  dls001_1.nst = 0;
1795  dls001_1.nje = 0;
1796  dls001_1.nslast = 0;
1797  dls001_1.hu = 0.;
1798  dls001_1.nqu = 0;
1799  dlsa01_1.mused = 0;
1800  dls001_1.miter = 0;
1801  dls001_1.ccmax = .3;
1802  dls001_1.maxcor = 3;
1803  dls001_1.msbp = 20;
1804  dls001_1.mxncf = 10;
1805  /* Initial call to F. (LF0 points to YH(*,2).) ------------------------- */
1806  lf0 = dls001_1.lyh + dls001_1.nyh;
1807  (*f)(&neq[1], t, &y[1], &rwork[lf0]);
1808  dls001_1.nfe = 1;
1809  /* Load the initial value vector in YH. --------------------------------- */
1810  i__1 = dls001_1.n;
1811 
1812  for (i__ = 1; i__ <= i__1; ++i__)
1813  {
1814  /* L115: */
1815  rwork[i__ + dls001_1.lyh - 1] = y[i__];
1816  }
1817 
1818  /* Load and invert the EWT array. (H is temporarily set to 1.0.) ------- */
1819  dls001_1.nq = 1;
1820  dls001_1.h__ = 1.;
1821  dewset_(&dls001_1.n, itol, &rtol[1], &atol[1], &rwork[dls001_1.lyh], &
1822  rwork[dls001_1.lewt]);
1823  i__1 = dls001_1.n;
1824 
1825  for (i__ = 1; i__ <= i__1; ++i__)
1826  {
1827  if (rwork[i__ + dls001_1.lewt - 1] <= 0.)
1828  {
1829  goto L621;
1830  }
1831 
1832  /* L120: */
1833  rwork[i__ + dls001_1.lewt - 1] = 1. / rwork[i__ + dls001_1.lewt - 1];
1834  }
1835 
1836  /* ----------------------------------------------------------------------- */
1837  /* The coding below computes the step size, H0, to be attempted on the */
1838  /* first step, unless the user has supplied a value for this. */
1839  /* First check that TOUT - T differs significantly from zero. */
1840  /* A scalar tolerance quantity TOL is computed, as MAX(RTOL(i)) */
1841  /* if this is positive, or MAX(ATOL(i)/ABS(Y(i))) otherwise, adjusted */
1842  /* so as to be between 100*UROUND and 1.0E-3. */
1843  /* Then the computed value H0 is given by: */
1844 
1845  /* H0**(-2) = 1./(TOL * w0**2) + TOL * (norm(F))**2 */
1846 
1847  /* where w0 = MAX (ABS(T), ABS(TOUT)), */
1848 
1849  /* F = the initial value of the vector f(t,y), and */
1850 
1851  /* norm() = the weighted vector norm used throughout, given by */
1852 
1853  /* the DMNORM function routine, and weighted by the */
1854 
1855  /* tolerances initially loaded into the EWT array. */
1856 
1857  /* The sign of H0 is inferred from the initial values of TOUT and T. */
1858 
1859  /* ABS(H0) is made .le. ABS(TOUT-T) in any case. */
1860 
1861  /* ----------------------------------------------------------------------- */
1862  if (h0 != 0.)
1863  {
1864  goto L180;
1865  }
1866 
1867  tdist = (d__1 = *tout - *t, fabs(d__1));
1868  /* Computing MAX */
1869  d__1 = fabs(*t), d__2 = fabs(*tout);
1870  w0 = std::max(d__1, d__2);
1871 
1872  if (tdist < dls001_1.uround * 2. * w0)
1873  {
1874  goto L622;
1875  }
1876 
1877  tol = rtol[1];
1878 
1879  if (*itol <= 2)
1880  {
1881  goto L140;
1882  }
1883 
1884  i__1 = dls001_1.n;
1885 
1886  for (i__ = 1; i__ <= i__1; ++i__)
1887  {
1888  /* L130: */
1889  /* Computing MAX */
1890  d__1 = tol, d__2 = rtol[i__];
1891  tol = std::max(d__1, d__2);
1892  }
1893 
1894 L140:
1895 
1896  if (tol > 0.)
1897  {
1898  goto L160;
1899  }
1900 
1901  atoli = atol[1];
1902  i__1 = dls001_1.n;
1903 
1904  for (i__ = 1; i__ <= i__1; ++i__)
1905  {
1906  if (*itol == 2 || *itol == 4)
1907  {
1908  atoli = atol[i__];
1909  }
1910 
1911  ayi = (d__1 = y[i__], fabs(d__1));
1912 
1913  if (ayi != 0.)
1914  {
1915  /* Computing MAX */
1916  d__1 = tol, d__2 = atoli / ayi;
1917  tol = std::max(d__1, d__2);
1918  }
1919 
1920  /* L150: */
1921  }
1922 
1923 L160:
1924  /* Computing MAX */
1925  d__1 = tol, d__2 = dls001_1.uround * 100.;
1926  tol = std::max(d__1, d__2);
1927  tol = std::min(tol, .001);
1928  sum = dmnorm_(&dls001_1.n, &rwork[lf0], &rwork[dls001_1.lewt]);
1929  /* Computing 2nd power */
1930  d__1 = sum;
1931  sum = 1. / (tol * w0 * w0) + tol * (d__1 * d__1);
1932  h0 = 1. / sqrt(sum);
1933  h0 = std::min(h0, tdist);
1934  d__1 = *tout - *t;
1935  h0 = d_sign(h0, d__1);
1936  /* Adjust H0 if necessary to meet HMAX bound. --------------------------- */
1937 L180:
1938  rh = fabs(h0) * dls001_1.hmxi;
1939 
1940  if (rh > 1.)
1941  {
1942  h0 /= rh;
1943  }
1944 
1945  /* Load H with H0 and scale YH(*,2) by H0. ------------------------------ */
1946  dls001_1.h__ = h0;
1947  i__1 = dls001_1.n;
1948 
1949  for (i__ = 1; i__ <= i__1; ++i__)
1950  {
1951  /* L190: */
1952  rwork[i__ + lf0 - 1] = h0 * rwork[i__ + lf0 - 1];
1953  }
1954 
1955  /* Check for a zero of g at T. ------------------------------------------ */
1956  dlsr01_1.irfnd = 0;
1957  dlsr01_1.toutc = *tout;
1958 
1959  if (dlsr01_1.ngc == 0)
1960  {
1961  goto L270;
1962  }
1963 
1964  drchek_(&c__1, (evalG)g, &neq[1], &y[1], &rwork[dls001_1.lyh], &
1965  dls001_1.nyh, &rwork[dlsr01_1.lg0], &rwork[dlsr01_1.lg1], &rwork[
1966  dlsr01_1.lgx], &jroot[1], &irt);
1967 
1968  if (irt == 0)
1969  {
1970  goto L270;
1971  }
1972 
1973  goto L632;
1974  /* ----------------------------------------------------------------------- */
1975  /* Block D. */
1976  /* The next code block is for continuation calls only (ISTATE = 2 or 3) */
1977  /* and is to check stop conditions before taking a step. */
1978  /* First, DRCHEK is called to check for a root within the last step */
1979  /* taken, other than the last root found there, if any. */
1980  /* If ITASK = 2 or 5, and y(TN) has not yet been returned to the user */
1981  /* because of an intervening root, return through Block G. */
1982  /* ----------------------------------------------------------------------- */
1983 L200:
1984  dls001_1.nslast = dls001_1.nst;
1985 
1986  irfp = dlsr01_1.irfnd;
1987 
1988  if (dlsr01_1.ngc == 0)
1989  {
1990  goto L205;
1991  }
1992 
1993  if (*itask == 1 || *itask == 4)
1994  {
1995  dlsr01_1.toutc = *tout;
1996  }
1997 
1998  drchek_(&c__2, (evalG)g, &neq[1], &y[1], &rwork[dls001_1.lyh], &
1999  dls001_1.nyh, &rwork[dlsr01_1.lg0], &rwork[dlsr01_1.lg1], &rwork[
2000  dlsr01_1.lgx], &jroot[1], &irt);
2001 
2002  if (irt != 1)
2003  {
2004  goto L205;
2005  }
2006 
2007  dlsr01_1.irfnd = 1;
2008  *istate = 3;
2009  *t = dlsr01_1.t0;
2010  goto L425;
2011 L205:
2012  dlsr01_1.irfnd = 0;
2013 
2014  if (irfp == 1 && dlsr01_1.tlast != dls001_1.tn && *itask == 2)
2015  {
2016  goto L400;
2017  }
2018 
2019  switch (*itask)
2020  {
2021  case 1: goto L210;
2022  case 2: goto L250;
2023  case 3: goto L220;
2024  case 4: goto L230;
2025  case 5: goto L240;
2026  }
2027 
2028 L210:
2029 
2030  if ((dls001_1.tn - *tout) * dls001_1.h__ < 0.)
2031  {
2032  goto L250;
2033  }
2034 
2035  dintdy_(tout, &c__0, &rwork[dls001_1.lyh], &dls001_1.nyh, &y[1], &iflag);
2036 
2037  if (iflag != 0)
2038  {
2039  goto L627;
2040  }
2041 
2042  *t = *tout;
2043  goto L420;
2044 L220:
2045  tp = dls001_1.tn - dls001_1.hu * (dls001_1.uround * 100. + 1.);
2046 
2047  if ((tp - *tout) * dls001_1.h__ > 0.)
2048  {
2049  goto L623;
2050  }
2051 
2052  if ((dls001_1.tn - *tout) * dls001_1.h__ < 0.)
2053  {
2054  goto L250;
2055  }
2056 
2057  *t = dls001_1.tn;
2058  goto L400;
2059 L230:
2060  tcrit = rwork[1];
2061 
2062  if ((dls001_1.tn - tcrit) * dls001_1.h__ > 0.)
2063  {
2064  goto L624;
2065  }
2066 
2067  if ((tcrit - *tout) * dls001_1.h__ < 0.)
2068  {
2069  goto L625;
2070  }
2071 
2072  if ((dls001_1.tn - *tout) * dls001_1.h__ < 0.)
2073  {
2074  goto L245;
2075  }
2076 
2077  dintdy_(tout, &c__0, &rwork[dls001_1.lyh], &dls001_1.nyh, &y[1], &iflag);
2078 
2079  if (iflag != 0)
2080  {
2081  goto L627;
2082  }
2083 
2084  *t = *tout;
2085  goto L420;
2086 L240:
2087  tcrit = rwork[1];
2088 
2089  if ((dls001_1.tn - tcrit) * dls001_1.h__ > 0.)
2090  {
2091  goto L624;
2092  }
2093 
2094 L245:
2095  hmx = fabs(dls001_1.tn) + fabs(dls001_1.h__);
2096  ihit = (d__1 = dls001_1.tn - tcrit, fabs(d__1)) <= dls001_1.uround * 100. *
2097  hmx;
2098 
2099  if (ihit)
2100  {
2101  *t = tcrit;
2102  }
2103 
2104  if (irfp == 1 && dlsr01_1.tlast != dls001_1.tn && *itask == 5)
2105  {
2106  goto L400;
2107  }
2108 
2109  if (ihit)
2110  {
2111  goto L400;
2112  }
2113 
2114  tnext = dls001_1.tn + dls001_1.h__ * (dls001_1.uround * 4. + 1.);
2115 
2116  if ((tnext - tcrit) * dls001_1.h__ <= 0.)
2117  {
2118  goto L250;
2119  }
2120 
2121  dls001_1.h__ = (tcrit - dls001_1.tn) * (1. - dls001_1.uround * 4.);
2122 
2123  if (*istate == 2 && dls001_1.jstart >= 0)
2124  {
2125  dls001_1.jstart = -2;
2126  }
2127 
2128  /* ----------------------------------------------------------------------- */
2129  /* Block E. */
2130  /* The next block is normally executed for all calls and contains */
2131  /* the call to the one-step core integrator DSTODA. */
2132 
2133  /* This is a looping point for the integration steps. */
2134 
2135  /* First check for too many steps being taken, update EWT (if not at */
2136  /* start of problem), check for too much accuracy being requested, and */
2137  /* check for H below the roundoff level in T. */
2138  /* ----------------------------------------------------------------------- */
2139 L250:
2140 
2141  if (dls001_1.meth == dlsa01_1.mused)
2142  {
2143  goto L255;
2144  }
2145 
2146  if (dlsa01_1.insufr == 1)
2147  {
2148  goto L550;
2149  }
2150 
2151  if (dlsa01_1.insufi == 1)
2152  {
2153  goto L555;
2154  }
2155 
2156 L255:
2157 
2158  if (dls001_1.nst - dls001_1.nslast >= dls001_1.mxstep)
2159  {
2160  goto L500;
2161  }
2162 
2163  dewset_(&dls001_1.n, itol, &rtol[1], &atol[1], &rwork[dls001_1.lyh], &
2164  rwork[dls001_1.lewt]);
2165  i__1 = dls001_1.n;
2166 
2167  for (i__ = 1; i__ <= i__1; ++i__)
2168  {
2169  if (rwork[i__ + dls001_1.lewt - 1] <= 0.)
2170  {
2171  goto L510;
2172  }
2173 
2174  /* L260: */
2175  rwork[i__ + dls001_1.lewt - 1] = 1. / rwork[i__ + dls001_1.lewt - 1];
2176  }
2177 
2178 L270:
2179  tolsf = dls001_1.uround * dmnorm_(&dls001_1.n, &rwork[dls001_1.lyh], &
2180  rwork[dls001_1.lewt]);
2181 
2182  if (tolsf <= 1.)
2183  {
2184  goto L280;
2185  }
2186 
2187  tolsf *= 2.;
2188 
2189  if (dls001_1.nst == 0)
2190  {
2191  goto L626;
2192  }
2193 
2194  goto L520;
2195 L280:
2196 
2197  if (dls001_1.tn + dls001_1.h__ != dls001_1.tn)
2198  {
2199  goto L290;
2200  }
2201 
2202  ++dls001_1.nhnil;
2203 
2204  if (dls001_1.nhnil > dls001_1.mxhnil)
2205  {
2206  goto L290;
2207  }
2208 
2209  msg = "DLSODAR- Warning..Internal T(=R1) and H(=R2) are ";
2210  mxerrwd(msg, &c__50, &c__101, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2211  c_b76, (C_INT)60);
2212  msg = " such that in the machine, T + H = T on the next step ";
2213  mxerrwd(msg, &c__60, &c__101, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2214  c_b76, (C_INT)60);
2215  msg = " (H = step size). Solver will continue anyway.";
2216  mxerrwd(msg, &c__50, &c__101, &c__0, &c__0, &c__0, &c__0, &c__2, &
2217  dls001_1.tn, &dls001_1.h__, (C_INT)60);
2218 
2219  if (dls001_1.nhnil < dls001_1.mxhnil)
2220  {
2221  goto L290;
2222  }
2223 
2224  msg = "DLSODAR- Above warning has been issued I1 times. ";
2225  mxerrwd(msg, &c__50, &c__102, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2226  c_b76, (C_INT)60);
2227  msg = " It will not be issued again for this problem.";
2228  mxerrwd(msg, &c__50, &c__102, &c__0, &c__1, &dls001_1.mxhnil, &c__0, &
2229  c__0, &c_b76, &c_b76, (C_INT)60);
2230 L290:
2231  /* ----------------------------------------------------------------------- */
2232  /* CALL DSTODA(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,DPRJA,DSOLSY) */
2233  /* ----------------------------------------------------------------------- */
2234  dstoda_(&neq[1], &y[1], &rwork[dls001_1.lyh], &dls001_1.nyh, &rwork[
2235  dls001_1.lyh], &rwork[dls001_1.lewt], &rwork[dls001_1.lsavf], &
2236  rwork[dls001_1.lacor], &rwork[dls001_1.lwm], &iwork[dls001_1.liwm]
2237  , f, jac, mpPJAC, mpSLVS);
2238  kgo = 1 - dls001_1.kflag;
2239 
2240  switch (kgo)
2241  {
2242  case 1: goto L300;
2243  case 2: goto L530;
2244  case 3: goto L540;
2245  }
2246 
2247  /* ----------------------------------------------------------------------- */
2248  /* Block F. */
2249  /* The following block handles the case of a successful return from the */
2250  /* core integrator (KFLAG = 0). */
2251  /* If a method switch was just made, record TSW, reset MAXORD, */
2252  /* set JSTART to -1 to signal DSTODA to complete the switch, */
2253  /* and do extra printing of data if IXPR = 1. */
2254  /* Then call DRCHEK to check for a root within the last step. */
2255  /* Then, if no root was found, check for stop conditions. */
2256  /* ----------------------------------------------------------------------- */
2257 L300:
2258  dls001_1.init = 1;
2259 
2260  if (dls001_1.meth == dlsa01_1.mused)
2261  {
2262  goto L310;
2263  }
2264 
2265  dlsa01_1.tsw = dls001_1.tn;
2266  dls001_1.maxord = dlsa01_1.mxordn;
2267 
2268  if (dls001_1.meth == 2)
2269  {
2270  dls001_1.maxord = dlsa01_1.mxords;
2271  }
2272 
2273  if (dls001_1.meth == 2)
2274  {
2275  rwork[dls001_1.lwm] = sqrt(dls001_1.uround);
2276  }
2277 
2278  dlsa01_1.insufr = std::min(dlsa01_1.insufr, (C_INT)1);
2279  dlsa01_1.insufi = std::min(dlsa01_1.insufi, (C_INT)1);
2280  dls001_1.jstart = -1;
2281 
2282  if (dlsa01_1.ixpr == 0)
2283  {
2284  goto L310;
2285  }
2286 
2287  if (dls001_1.meth == 2)
2288  {
2289  msg = "DLSODAR- A switch to the BDF (stiff) method has occurred";
2290  mxerrwd(msg, &c__60, &c__105, &c__0, &c__0, &c__0, &c__0, &c__0, &
2291  c_b76, &c_b76, (C_INT)60);
2292  }
2293 
2294  if (dls001_1.meth == 1)
2295  {
2296  msg = "DLSODAR- A switch to the Adams (nonstiff) method occurred";
2297  mxerrwd(msg, &c__60, &c__106, &c__0, &c__0, &c__0, &c__0, &c__0, &
2298  c_b76, &c_b76, (C_INT)60);
2299  }
2300 
2301  msg = " at T = R1, tentative step size H = R2, step NST = I1";
2302  mxerrwd(msg, &c__60, &c__107, &c__0, &c__1, &dls001_1.nst, &c__0, &c__2, &
2303  dls001_1.tn, &dls001_1.h__, (C_INT)60);
2304 L310:
2305 
2306  if (dlsr01_1.ngc == 0)
2307  {
2308  goto L315;
2309  }
2310 
2311  drchek_(&c__3, (evalG)g, &neq[1], &y[1], &rwork[dls001_1.lyh], &
2312  dls001_1.nyh, &rwork[dlsr01_1.lg0], &rwork[dlsr01_1.lg1], &rwork[
2313  dlsr01_1.lgx], &jroot[1], &irt);
2314 
2315  if (irt != 1)
2316  {
2317  goto L315;
2318  }
2319 
2320  dlsr01_1.irfnd = 1;
2321  *istate = 3;
2322  *t = dlsr01_1.t0;
2323  goto L425;
2324 L315:
2325 
2326  switch (*itask)
2327  {
2328  case 1: goto L320;
2329  case 2: goto L400;
2330  case 3: goto L330;
2331  case 4: goto L340;
2332  case 5: goto L350;
2333  }
2334 
2335  /* ITASK = 1. If TOUT has been reached, interpolate. ------------------- */
2336 L320:
2337 
2338  if ((dls001_1.tn - *tout) * dls001_1.h__ < 0.)
2339  {
2340  goto L250;
2341  }
2342 
2343  dintdy_(tout, &c__0, &rwork[dls001_1.lyh], &dls001_1.nyh, &y[1], &iflag);
2344  *t = *tout;
2345  goto L420;
2346  /* ITASK = 3. Jump to exit if TOUT was reached. ------------------------ */
2347 L330:
2348 
2349  if ((dls001_1.tn - *tout) * dls001_1.h__ >= 0.)
2350  {
2351  goto L400;
2352  }
2353 
2354  goto L250;
2355  /* ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary. */
2356 L340:
2357 
2358  if ((dls001_1.tn - *tout) * dls001_1.h__ < 0.)
2359  {
2360  goto L345;
2361  }
2362 
2363  dintdy_(tout, &c__0, &rwork[dls001_1.lyh], &dls001_1.nyh, &y[1], &iflag);
2364  *t = *tout;
2365  goto L420;
2366 L345:
2367  hmx = fabs(dls001_1.tn) + fabs(dls001_1.h__);
2368  ihit = (d__1 = dls001_1.tn - tcrit, fabs(d__1)) <= dls001_1.uround * 100. *
2369  hmx;
2370 
2371  if (ihit)
2372  {
2373  goto L400;
2374  }
2375 
2376  tnext = dls001_1.tn + dls001_1.h__ * (dls001_1.uround * 4. + 1.);
2377 
2378  if ((tnext - tcrit) * dls001_1.h__ <= 0.)
2379  {
2380  goto L250;
2381  }
2382 
2383  dls001_1.h__ = (tcrit - dls001_1.tn) * (1. - dls001_1.uround * 4.);
2384 
2385  if (dls001_1.jstart >= 0)
2386  {
2387  dls001_1.jstart = -2;
2388  }
2389 
2390  goto L250;
2391  /* ITASK = 5. See if TCRIT was reached and jump to exit. --------------- */
2392 L350:
2393  hmx = fabs(dls001_1.tn) + fabs(dls001_1.h__);
2394  ihit = (d__1 = dls001_1.tn - tcrit, fabs(d__1)) <= dls001_1.uround * 100. *
2395  hmx;
2396  /* ----------------------------------------------------------------------- */
2397  /* Block G. */
2398  /* The following block handles all successful returns from DLSODAR. */
2399  /* If ITASK .ne. 1, Y is loaded from YH and T is set accordingly. */
2400  /* ISTATE is set to 2, and the optional outputs are loaded into the */
2401  /* work arrays before returning. */
2402  /* ----------------------------------------------------------------------- */
2403 L400:
2404  i__1 = dls001_1.n;
2405 
2406  for (i__ = 1; i__ <= i__1; ++i__)
2407  {
2408  /* L410: */
2409  y[i__] = rwork[i__ + dls001_1.lyh - 1];
2410  }
2411 
2412  *t = dls001_1.tn;
2413 
2414  if (*itask != 4 && *itask != 5)
2415  {
2416  goto L420;
2417  }
2418 
2419  if (ihit)
2420  {
2421  *t = tcrit;
2422  }
2423 
2424 L420:
2425  *istate = 2;
2426 L425:
2427  rwork[11] = dls001_1.hu;
2428  rwork[12] = dls001_1.h__;
2429  rwork[13] = dls001_1.tn;
2430  rwork[15] = dlsa01_1.tsw;
2431  iwork[11] = dls001_1.nst;
2432  iwork[12] = dls001_1.nfe;
2433  iwork[13] = dls001_1.nje;
2434  iwork[14] = dls001_1.nqu;
2435  iwork[15] = dls001_1.nq;
2436  iwork[19] = dlsa01_1.mused;
2437  iwork[20] = dls001_1.meth;
2438  iwork[10] = dlsr01_1.nge;
2439  dlsr01_1.tlast = *t;
2440  return 0;
2441  /* ----------------------------------------------------------------------- */
2442  /* Block H. */
2443  /* The following block handles all unsuccessful returns other than */
2444  /* those for illegal input. First the error message routine is called. */
2445  /* If there was an error test or convergence test failure, IMXER is set. */
2446  /* Then Y is loaded from YH and T is set to TN. */
2447  /* The optional outputs are loaded into the work arrays before returning. */
2448  /* ----------------------------------------------------------------------- */
2449  /* The maximum number of steps was taken before reaching TOUT. ---------- */
2450 L500:
2451  msg = "DLSODAR- At current T (=R1), MXSTEP (=I1) steps ";
2452  mxerrwd(msg, &c__50, &c__201, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2453  c_b76, (C_INT)60);
2454  msg = " taken on this call before reaching TOUT ";
2455  mxerrwd(msg, &c__50, &c__201, &c__0, &c__1, &dls001_1.mxstep, &c__0, &
2456  c__1, &dls001_1.tn, &c_b76, (C_INT)60);
2457  *istate = -1;
2458  goto L580;
2459  /* EWT(i) .le. 0.0 for some i (not at start of problem). ---------------- */
2460 L510:
2461  ewti = rwork[dls001_1.lewt + i__ - 1];
2462  msg = "DLSODAR- At T(=R1), EWT(I1) has become R2 .le. 0.";
2463  mxerrwd(msg, &c__50, &c__202, &c__0, &c__1, &i__, &c__0, &c__2, &
2464  dls001_1.tn, &ewti, (C_INT)60);
2465  *istate = -6;
2466  goto L580;
2467  /* Too much accuracy requested for machine precision. ------------------- */
2468 L520:
2469  msg = "DLSODAR- At T (=R1), too much accuracy requested ";
2470  mxerrwd(msg, &c__50, &c__203, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2471  c_b76, (C_INT)60);
2472  msg = " for precision of machine.. See TOLSF (=R2) ";
2473  mxerrwd(msg, &c__50, &c__203, &c__0, &c__0, &c__0, &c__0, &c__2, &
2474  dls001_1.tn, &tolsf, (C_INT)60);
2475  rwork[14] = tolsf;
2476  *istate = -2;
2477  goto L580;
2478  /* KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. ----- */
2479 L530:
2480  msg = "DLSODAR- At T(=R1), step size H(=R2), the error ";
2481  mxerrwd(msg, &c__50, &c__204, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2482  c_b76, (C_INT)60);
2483  msg = " test failed repeatedly or with ABS(H) = HMIN";
2484  mxerrwd(msg, &c__50, &c__204, &c__0, &c__0, &c__0, &c__0, &c__2, &
2485  dls001_1.tn, &dls001_1.h__, (C_INT)60);
2486  *istate = -4;
2487  goto L560;
2488  /* KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ---- */
2489 L540:
2490  msg = "DLSODAR- At T (=R1) and step size H (=R2), the ";
2491  mxerrwd(msg, &c__50, &c__205, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2492  c_b76, (C_INT)60);
2493  msg = " corrector convergence failed repeatedly ";
2494  mxerrwd(msg, &c__50, &c__205, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2495  c_b76, (C_INT)60);
2496  msg = " or with ABS(H) = HMIN ";
2497  mxerrwd(msg, &c__30, &c__205, &c__0, &c__0, &c__0, &c__0, &c__2, &
2498  dls001_1.tn, &dls001_1.h__, (C_INT)60);
2499  *istate = -5;
2500  goto L560;
2501  /* RWORK length too small to proceed. ----------------------------------- */
2502 L550:
2503  msg = "DLSODAR- At current T(=R1), RWORK length too small";
2504  mxerrwd(msg, &c__50, &c__206, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2505  c_b76, (C_INT)60);
2506  msg = " to proceed. The integration was otherwise successful.";
2507  mxerrwd(msg, &c__60, &c__206, &c__0, &c__0, &c__0, &c__0, &c__1, &
2508  dls001_1.tn, &c_b76, (C_INT)60);
2509  *istate = -7;
2510  goto L580;
2511  /* IWORK length too small to proceed. ----------------------------------- */
2512 L555:
2513  msg = "DLSODAR- At current T(=R1), IWORK length too small";
2514  mxerrwd(msg, &c__50, &c__207, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2515  c_b76, (C_INT)60);
2516  msg = " to proceed. The integration was otherwise successful.";
2517  mxerrwd(msg, &c__60, &c__207, &c__0, &c__0, &c__0, &c__0, &c__1, &
2518  dls001_1.tn, &c_b76, (C_INT)60);
2519  *istate = -7;
2520  goto L580;
2521  /* Compute IMXER if relevant. ------------------------------------------- */
2522 L560:
2523  big = 0.;
2524  imxer = 1;
2525  i__1 = dls001_1.n;
2526 
2527  for (i__ = 1; i__ <= i__1; ++i__)
2528  {
2529  size = (d__1 = rwork[i__ + dls001_1.lacor - 1] * rwork[i__ +
2530  dls001_1.lewt - 1], fabs(d__1));
2531 
2532  if (big >= size)
2533  {
2534  goto L570;
2535  }
2536 
2537  big = size;
2538  imxer = i__;
2539 L570:
2540  ;
2541  }
2542 
2543  iwork[16] = imxer;
2544  /* Set Y vector, T, and optional outputs. ------------------------------- */
2545 L580:
2546  i__1 = dls001_1.n;
2547 
2548  for (i__ = 1; i__ <= i__1; ++i__)
2549  {
2550  /* L590: */
2551  y[i__] = rwork[i__ + dls001_1.lyh - 1];
2552  }
2553 
2554  *t = dls001_1.tn;
2555  rwork[11] = dls001_1.hu;
2556  rwork[12] = dls001_1.h__;
2557  rwork[13] = dls001_1.tn;
2558  rwork[15] = dlsa01_1.tsw;
2559  iwork[11] = dls001_1.nst;
2560  iwork[12] = dls001_1.nfe;
2561  iwork[13] = dls001_1.nje;
2562  iwork[14] = dls001_1.nqu;
2563  iwork[15] = dls001_1.nq;
2564  iwork[19] = dlsa01_1.mused;
2565  iwork[20] = dls001_1.meth;
2566  iwork[10] = dlsr01_1.nge;
2567  dlsr01_1.tlast = *t;
2568  return 0;
2569  /* ----------------------------------------------------------------------- */
2570  /* Block I. */
2571  /* The following block handles all error returns due to illegal input */
2572  /* (ISTATE = -3), as detected before calling the core integrator. */
2573  /* First the error message routine is called. If the illegal input */
2574  /* is a negative ISTATE, the run is aborted (apparent infinite loop). */
2575  /* ----------------------------------------------------------------------- */
2576 L601:
2577  msg = "DLSODAR- ISTATE(=I1) illegal.";
2578  mxerrwd(msg, &c__30, &c__1, &c__0, &c__1, istate, &c__0, &c__0, &c_b76, &
2579  c_b76, (C_INT)60);
2580 
2581  if (*istate < 0)
2582  {
2583  goto L800;
2584  }
2585 
2586  goto L700;
2587 L602:
2588  msg = "DLSODAR- ITASK (=I1) illegal.";
2589  mxerrwd(msg, &c__30, &c__2, &c__0, &c__1, itask, &c__0, &c__0, &c_b76, &
2590  c_b76, (C_INT)60);
2591  goto L700;
2592 L603:
2593  msg = "DLSODAR- ISTATE.gt.1 but DLSODAR not initialized.";
2594  mxerrwd(msg, &c__50, &c__3, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2595  c_b76, (C_INT)60);
2596  goto L700;
2597 L604:
2598  msg = "DLSODAR- NEQ (=I1) .lt. 1 ";
2599  mxerrwd(msg, &c__30, &c__4, &c__0, &c__1, &neq[1], &c__0, &c__0, &c_b76, &
2600  c_b76, (C_INT)60);
2601  goto L700;
2602 L605:
2603  msg = "DLSODAR- ISTATE = 3 and NEQ increased (I1 to I2).";
2604  mxerrwd(msg, &c__50, &c__5, &c__0, &c__2, &dls001_1.n, &neq[1], &c__0, &
2605  c_b76, &c_b76, (C_INT)60);
2606  goto L700;
2607 L606:
2608  msg = "DLSODAR- ITOL (=I1) illegal. ";
2609  mxerrwd(msg, &c__30, &c__6, &c__0, &c__1, itol, &c__0, &c__0, &c_b76, &
2610  c_b76, (C_INT)60);
2611  goto L700;
2612 L607:
2613  msg = "DLSODAR- IOPT (=I1) illegal. ";
2614  mxerrwd(msg, &c__30, &c__7, &c__0, &c__1, iopt, &c__0, &c__0, &c_b76, &
2615  c_b76, (C_INT)60);
2616  goto L700;
2617 L608:
2618  msg = "DLSODAR- JT (=I1) illegal. ";
2619  mxerrwd(msg, &c__30, &c__8, &c__0, &c__1, jt, &c__0, &c__0, &c_b76, &
2620  c_b76, (C_INT)60);
2621  goto L700;
2622 L609:
2623  msg = "DLSODAR- ML (=I1) illegal: .lt.0 or .ge.NEQ (=I2)";
2624  mxerrwd(msg, &c__50, &c__9, &c__0, &c__2, &ml, &neq[1], &c__0, &c_b76, &
2625  c_b76, (C_INT)60);
2626  goto L700;
2627 L610:
2628  msg = "DLSODAR- MU (=I1) illegal: .lt.0 or .ge.NEQ (=I2)";
2629  mxerrwd(msg, &c__50, &c__10, &c__0, &c__2, &mu, &neq[1], &c__0, &c_b76, &
2630  c_b76, (C_INT)60);
2631  goto L700;
2632 L611:
2633  msg = "DLSODAR- IXPR (=I1) illegal. ";
2634  mxerrwd(msg, &c__30, &c__11, &c__0, &c__1, &dlsa01_1.ixpr, &c__0, &c__0, &
2635  c_b76, &c_b76, (C_INT)60);
2636  goto L700;
2637 L612:
2638  msg = "DLSODAR- MXSTEP (=I1) .lt. 0 ";
2639  mxerrwd(msg, &c__30, &c__12, &c__0, &c__1, &dls001_1.mxstep, &c__0, &c__0,
2640  &c_b76, &c_b76, (C_INT)60);
2641  goto L700;
2642 L613:
2643  msg = "DLSODAR- MXHNIL (=I1) .lt. 0 ";
2644  mxerrwd(msg, &c__30, &c__13, &c__0, &c__1, &dls001_1.mxhnil, &c__0, &c__0,
2645  &c_b76, &c_b76, (C_INT)60);
2646  goto L700;
2647 L614:
2648  msg = "DLSODAR- TOUT (=R1) behind T (=R2) ";
2649  mxerrwd(msg, &c__40, &c__14, &c__0, &c__0, &c__0, &c__0, &c__2, tout, t, (
2650  C_INT)60);
2651  msg = " Integration direction is given by H0 (=R1) ";
2652  mxerrwd(msg, &c__50, &c__14, &c__0, &c__0, &c__0, &c__0, &c__1, &h0, &
2653  c_b76, (C_INT)60);
2654  goto L700;
2655 L615:
2656  msg = "DLSODAR- HMAX (=R1) .lt. 0.0 ";
2657  mxerrwd(msg, &c__30, &c__15, &c__0, &c__0, &c__0, &c__0, &c__1, &hmax, &
2658  c_b76, (C_INT)60);
2659  goto L700;
2660 L616:
2661  msg = "DLSODAR- HMIN (=R1) .lt. 0.0 ";
2662  mxerrwd(msg, &c__30, &c__16, &c__0, &c__0, &c__0, &c__0, &c__1, &
2663  dls001_1.hmin, &c_b76, (C_INT)60);
2664  goto L700;
2665 L617:
2666  msg = "DLSODAR- RWORK length needed, LENRW(=I1), exceeds LRW(=I2) ";
2667  mxerrwd(msg, &c__60, &c__17, &c__0, &c__2, &lenrw, lrw, &c__0, &c_b76, &
2668  c_b76, (C_INT)60);
2669  goto L700;
2670 L618:
2671  msg = "DLSODAR- IWORK length needed, LENIW(=I1), exceeds LIW(=I2) ";
2672  mxerrwd(msg, &c__60, &c__18, &c__0, &c__2, &leniw, liw, &c__0, &c_b76, &
2673  c_b76, (C_INT)60);
2674  goto L700;
2675 L619:
2676  msg = "DLSODAR- RTOL(I1) is R1 .lt. 0.0 ";
2677  mxerrwd(msg, &c__40, &c__19, &c__0, &c__1, &i__, &c__0, &c__1, &rtoli, &
2678  c_b76, (C_INT)60);
2679  goto L700;
2680 L620:
2681  msg = "DLSODAR- ATOL(I1) is R1 .lt. 0.0 ";
2682  mxerrwd(msg, &c__40, &c__20, &c__0, &c__1, &i__, &c__0, &c__1, &atoli, &
2683  c_b76, (C_INT)60);
2684  goto L700;
2685 L621:
2686  ewti = rwork[dls001_1.lewt + i__ - 1];
2687  msg = "DLSODAR- EWT(I1) is R1 .le. 0.0 ";
2688  mxerrwd(msg, &c__40, &c__21, &c__0, &c__1, &i__, &c__0, &c__1, &ewti, &
2689  c_b76, (C_INT)60);
2690  goto L700;
2691 L622:
2692  msg = "DLSODAR- TOUT(=R1) too close to T(=R2) to start integration.";
2693  mxerrwd(msg, &c__60, &c__22, &c__0, &c__0, &c__0, &c__0, &c__2, tout, t, (
2694  C_INT)60);
2695  goto L700;
2696 L623:
2697  msg = "DLSODAR- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ";
2698  mxerrwd(msg, &c__60, &c__23, &c__0, &c__1, itask, &c__0, &c__2, tout, &tp,
2699  (C_INT)60);
2700  goto L700;
2701 L624:
2702  msg = "DLSODAR- ITASK = 4 or 5 and TCRIT (=R1) behind TCUR (=R2) ";
2703  mxerrwd(msg, &c__60, &c__24, &c__0, &c__0, &c__0, &c__0, &c__2, &tcrit, &
2704  dls001_1.tn, (C_INT)60);
2705  goto L700;
2706 L625:
2707  msg = "DLSODAR- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ";
2708  mxerrwd(msg, &c__60, &c__25, &c__0, &c__0, &c__0, &c__0, &c__2, &tcrit,
2709  tout, (C_INT)60);
2710  goto L700;
2711 L626:
2712  msg = "DLSODAR- At start of problem, too much accuracy ";
2713  mxerrwd(msg, &c__50, &c__26, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2714  c_b76, (C_INT)60);
2715  msg = " requested for precision of machine.. See TOLSF (=R1) ";
2716  mxerrwd(msg, &c__60, &c__26, &c__0, &c__0, &c__0, &c__0, &c__1, &tolsf, &
2717  c_b76, (C_INT)60);
2718  rwork[14] = tolsf;
2719  goto L700;
2720 L627:
2721  msg = "DLSODAR- Trouble in DINTDY. ITASK = I1, TOUT = R1";
2722  mxerrwd(msg, &c__50, &c__27, &c__0, &c__1, itask, &c__0, &c__1, tout, &
2723  c_b76, (C_INT)60);
2724  goto L700;
2725 L628:
2726  msg = "DLSODAR- MXORDN (=I1) .lt. 0 ";
2727  mxerrwd(msg, &c__30, &c__28, &c__0, &c__1, &dlsa01_1.mxordn, &c__0, &c__0,
2728  &c_b76, &c_b76, (C_INT)60);
2729  goto L700;
2730 L629:
2731  msg = "DLSODAR- MXORDS (=I1) .lt. 0 ";
2732  mxerrwd(msg, &c__30, &c__29, &c__0, &c__1, &dlsa01_1.mxords, &c__0, &c__0,
2733  &c_b76, &c_b76, (C_INT)60);
2734  goto L700;
2735 L630:
2736  msg = "DLSODAR- NG (=I1) .lt. 0 ";
2737  mxerrwd(msg, &c__30, &c__30, &c__0, &c__1, ng, &c__0, &c__0, &c_b76, &
2738  c_b76, (C_INT)60);
2739  goto L700;
2740 L631:
2741  msg = "DLSODAR- NG changed (from I1 to I2) illegally, ";
2742  mxerrwd(msg, &c__50, &c__31, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2743  c_b76, (C_INT)60);
2744  msg = " i.e. not immediately after a root was found.";
2745  mxerrwd(msg, &c__50, &c__31, &c__0, &c__2, &dlsr01_1.ngc, ng, &c__0, &
2746  c_b76, &c_b76, (C_INT)60);
2747  goto L700;
2748 L632:
2749  msg = "DLSODAR- One or more components of g has a root ";
2750  mxerrwd(msg, &c__50, &c__32, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2751  c_b76, (C_INT)60);
2752  msg = " too near to the initial point. ";
2753  mxerrwd(msg, &c__40, &c__32, &c__0, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2754  c_b76, (C_INT)60);
2755  *istate = -33;
2756  return 0;
2757 
2758 L700:
2759  *istate = -3;
2760  return 0;
2761 
2762 L800:
2763  msg = "DLSODAR- Run aborted.. apparent infinite loop. ";
2764  mxerrwd(msg, &c__50, &c__303, &c__2, &c__0, &c__0, &c__0, &c__0, &c_b76, &
2765  c_b76, (C_INT)60);
2766  return 0;
2767  /* ----------------------- End of Subroutine DLSODAR --------------------- */
2768 } /* dlsodar_ */
C_INT dstoda_(C_INT *neq, double *y, double *yh, C_INT *nyh, double *yh1, double *ewt, double *savf, double *acor, double *wm, C_INT *iwm, evalF f, evalJ jac, PJAC *pjac, SLVS *slvs)
Definition: dstoda.cpp:45
static const C_INT c__24
Definition: CLSODAR.cpp:68
#define C_INT
Definition: copasi.h:115
static const C_INT c__15
Definition: CLSODAR.cpp:59
static const C_INT c__50
Definition: CLSODAR.cpp:78
void(* evalG)(const C_INT *, const double *, const double *, const C_INT *, double *)
Definition: common.h:32
static const C_INT mxstp0
Definition: CLSODAR.h:55
static const C_INT c__40
Definition: CLSODAR.cpp:77
static const C_INT c__16
Definition: CLSODAR.cpp:60
static const C_INT c__105
Definition: CLSODAR.cpp:84
int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
SLVS * mpSLVS
Definition: CLSODAR.h:53
static const C_INT c__0
Definition: CLSODAR.cpp:44
static const C_INT c__202
Definition: CLSODAR.cpp:88
#define dls001_1
Definition: CLSODAR.cpp:38
static const C_INT c__27
Definition: CLSODAR.cpp:71
static const C_INT c__103
Definition: CLSODAR.cpp:82
static const C_INT c__303
Definition: CLSODAR.cpp:94
C_INT dintdy_(double *t, const C_INT *k, double *yh, C_INT *nyh, double *dky, C_INT *iflag)
Definition: dintdy.cpp:53
static const C_INT c__30
Definition: CLSODAR.cpp:74
static const C_INT c__205
Definition: CLSODAR.cpp:91
static const C_INT c__17
Definition: CLSODAR.cpp:61
static const C_INT c__106
Definition: CLSODAR.cpp:85
static const C_INT c__23
Definition: CLSODAR.cpp:67
static const C_INT c__8
Definition: CLSODAR.cpp:52
C_INT dewset_(C_INT *n, C_INT *itol, double *rtol, double *atol, double *ycur, double *ewt)
Definition: dewset.cpp:30
double d_sign(const double &a, const double &b)
Definition: CLSODA.cpp:2471
static const C_INT c__12
Definition: CLSODAR.cpp:56
static const C_INT c__3
Definition: CLSODAR.cpp:47
#define dlsa01_1
Definition: CLSODAR.cpp:39
static const C_INT c__204
Definition: CLSODAR.cpp:90
static const C_INT c__4
Definition: CLSODAR.cpp:48
static const C_INT c__101
Definition: CLSODAR.cpp:80
static const C_INT c__25
Definition: CLSODAR.cpp:69
static const C_INT c__1
Definition: CLSODAR.cpp:45
static const C_INT c__6
Definition: CLSODAR.cpp:50
static const C_INT c__9
Definition: CLSODAR.cpp:53
static const C_INT c__102
Definition: CLSODAR.cpp:81
static const C_INT c__7
Definition: CLSODAR.cpp:51
#define dlsr01_1
Definition: CLSODAR.cpp:40
static const C_INT c__11
Definition: CLSODAR.cpp:55
PJAC * mpPJAC
Definition: CLSODAR.h:52
static const C_INT c__207
Definition: CLSODAR.cpp:93
static const C_INT c__206
Definition: CLSODAR.cpp:92
static const C_INT c__32
Definition: CLSODAR.cpp:76
static const C_INT c__28
Definition: CLSODAR.cpp:72
static const C_INT c__60
Definition: CLSODAR.cpp:79
static const C_INT c__18
Definition: CLSODAR.cpp:62
static const C_INT c__22
Definition: CLSODAR.cpp:66
static const C_INT c__107
Definition: CLSODAR.cpp:86
static const C_INT c__29
Definition: CLSODAR.cpp:73
static const C_INT c__19
Definition: CLSODAR.cpp:63
static const C_INT c__20
Definition: CLSODAR.cpp:64
static const double c_b76
Definition: CLSODAR.cpp:42
static const C_INT mxhnl0
Definition: CLSODAR.h:56
double dmnorm_(C_INT *n, double *v, double *w)
Definition: dmnorm.cpp:35
static const C_INT c__201
Definition: CLSODAR.cpp:87
static const C_INT c__203
Definition: CLSODAR.cpp:89
static const C_INT c__5
Definition: CLSODAR.cpp:49
static const C_INT c__31
Definition: CLSODAR.cpp:75
static const C_INT mord[2]
Definition: CLSODAR.h:57
static const C_INT c__2
Definition: CLSODAR.cpp:46
if(!yymsg) yymsg
static const C_INT c__13
Definition: CLSODAR.cpp:57
static const C_INT c__26
Definition: CLSODAR.cpp:70
C_INT drchek_(const C_INT *job, evalG g, C_INT *neq, double *y, double *yh, C_INT *nyh, double *g0, double *g1, double *gx, C_INT *jroot, C_INT *irt)
Definition: drcheck.cpp:46
static const C_INT c__10
Definition: CLSODAR.cpp:54
static const C_INT c__21
Definition: CLSODAR.cpp:65
#define min(a, b)
Definition: f2c.h:175
static const C_INT c__104
Definition: CLSODAR.cpp:83
static const C_INT c__14
Definition: CLSODAR.cpp:58
#define max(a, b)
Definition: f2c.h:176

Member Data Documentation

const C_INT CLSODAR::mord = {12, 5}
staticprivate

Definition at line 57 of file CLSODAR.h.

Referenced by operator()().

PJAC* CLSODAR::mpPJAC
private

Definition at line 52 of file CLSODAR.h.

Referenced by CLSODAR(), operator()(), and ~CLSODAR().

SLVS* CLSODAR::mpSLVS
private

Definition at line 53 of file CLSODAR.h.

Referenced by CLSODAR(), operator()(), and ~CLSODAR().

const C_INT CLSODAR::mxhnl0 = 10
staticprivate

Definition at line 56 of file CLSODAR.h.

Referenced by operator()().

const C_INT CLSODAR::mxstp0 = 500
staticprivate

Definition at line 55 of file CLSODAR.h.

Referenced by operator()().


The documentation for this class was generated from the following files: