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