Porting continues.
[mipav.git] / mipav / src / gov / nih / mipav / model / algorithms / CVODES.java
1 package gov.nih.mipav.model.algorithms;\r
2 \r
3 import java.io.RandomAccessFile;\r
4 \r
5 import gov.nih.mipav.view.MipavUtil;\r
6 import gov.nih.mipav.view.Preferences;\r
7 \r
8 public abstract class CVODES {\r
9         // This is a port of code from CVODES, a stiff and non-stiff ODE solver, from C to Java\r
10         /*Copyright (c) 2002-2016, Lawrence Livermore National Security. \r
11         Produced at the Lawrence Livermore National Laboratory.\r
12         Written by A.C. Hindmarsh, D.R. Reynolds, R. Serban, C.S. Woodward,\r
13         S.D. Cohen, A.G. Taylor, S. Peles, L.E. Banks, and D. Shumaker.\r
14         LLNL-CODE-667205    (ARKODE)\r
15         UCRL-CODE-155951    (CVODE)\r
16         UCRL-CODE-155950    (CVODES)\r
17         UCRL-CODE-155952    (IDA)\r
18         UCRL-CODE-237203    (IDAS)\r
19         LLNL-CODE-665877    (KINSOL)\r
20         All rights reserved. \r
21 \r
22         This file is part of SUNDIALS.  For details, \r
23         see http://computation.llnl.gov/projects/sundials\r
24 \r
25         Redistribution and use in source and binary forms, with or without\r
26         modification, are permitted provided that the following conditions\r
27         are met:\r
28 \r
29         1. Redistributions of source code must retain the above copyright\r
30         notice, this list of conditions and the disclaimer below.\r
31 \r
32         2. Redistributions in binary form must reproduce the above copyright\r
33         notice, this list of conditions and the disclaimer (as noted below)\r
34         in the documentation and/or other materials provided with the\r
35         distribution.\r
36 \r
37         3. Neither the name of the LLNS/LLNL nor the names of its contributors\r
38         may be used to endorse or promote products derived from this software\r
39         without specific prior written permission.\r
40 \r
41         THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \r
42         "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT \r
43         LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS \r
44         FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL \r
45         LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF \r
46         ENERGY OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, \r
47         SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED \r
48         TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, \r
49         DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY \r
50         THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT \r
51         (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE \r
52         OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\r
53 \r
54         Additional BSD Notice\r
55         ---------------------\r
56         1. This notice is required to be provided under our contract with\r
57         the U.S. Department of Energy (DOE). This work was produced at\r
58         Lawrence Livermore National Laboratory under Contract \r
59         No. DE-AC52-07NA27344 with the DOE.\r
60 \r
61         2. Neither the United States Government nor Lawrence Livermore \r
62         National Security, LLC nor any of their employees, makes any warranty, \r
63         express or implied, or assumes any liability or responsibility for the\r
64         accuracy, completeness, or usefulness of any */\r
65         \r
66         /* Basic CVODES constants */\r
67 \r
68         final int ADAMS_Q_MAX = 12;      /* max value of q for lmm == ADAMS    */\r
69         final int BDF_Q_MAX = 5;      /* max value of q for lmm == BDF      */\r
70         final int  Q_MAX = ADAMS_Q_MAX;  /* max value of q for either lmm      */\r
71         final int L_MAX  = (Q_MAX+1);    /* max value of L for either lmm      */\r
72         final int NUM_TESTS = 5;      /* number of error test quantities    */\r
73         \r
74     double DBL_EPSILON;\r
75     double UNIT_ROUNDOFF;\r
76 \r
77 \r
78 \r
79          // lmm:   The user of the CVODES package specifies whether to use\r
80          // the CV_ADAMS or CV_BDF (backward differentiation formula)\r
81          // linear multistep method. The BDF method is recommended\r
82          // for stiff problems, and the CV_ADAMS method is recommended\r
83          // for nonstiff problems.\r
84         final int CV_ADAMS = 1;\r
85         final int CV_BDF = 2;\r
86         \r
87         //iter:  At each internal time step, a nonlinear equation must\r
88         //        be solved. The user can specify either CV_FUNCTIONAL\r
89         //        iteration, which does not require linear algebra, or a\r
90         //        CV_NEWTON iteration, which requires the solution of linear\r
91         //        systems. In the CV_NEWTON case, the user also specifies a\r
92         //        CVODE linear solver. CV_NEWTON is recommended in case of\r
93         //        stiff problems.\r
94         final int CV_FUNCTIONAL = 1;\r
95         final int CV_NEWTON = 2;\r
96         \r
97         //itask: The itask input parameter to CVode indicates the job\r
98         //        of the solver for the next user step. The CV_NORMAL\r
99         //        itask is to have the solver take internal steps until\r
100         //        it has reached or just passed the user specified tout\r
101         //        parameter. The solver then interpolates in order to\r
102         //        return an approximate value of y(tout). The CV_ONE_STEP\r
103         //        option tells the solver to just take one internal step\r
104         //        and return the solution at the point reached by that step.\r
105         final int CV_NORMAL = 1;\r
106         final int CV_ONE_STEP = 2;\r
107         \r
108         //ism:   This parameter specifies the sensitivity corrector type\r
109         //        to be used. In the CV_SIMULTANEOUS case, the nonlinear\r
110         //        systems for states and all sensitivities are solved\r
111         //        simultaneously. In the CV_STAGGERED case, the nonlinear\r
112         //        system for states is solved first and then, the\r
113         //        nonlinear systems for all sensitivities are solved\r
114         //        at the same time. Finally, in the CV_STAGGERED1 approach\r
115         //        all nonlinear systems are solved in a sequence.\r
116         final int CV_SIMULTANEOUS = 1;\r
117         final int CV_STAGGERED = 2;\r
118         final int CV_STAGGERED1 = 3;\r
119 \r
120 \r
121         \r
122         /*\r
123          * Control constants for type of sensitivity RHS\r
124          * ---------------------------------------------\r
125          */\r
126 \r
127         final int CV_ONESENS = 1;\r
128         final int CV_ALLSENS = 2;\r
129 \r
130         \r
131         /*\r
132          * Control constants for tolerances\r
133          * --------------------------------\r
134          */\r
135 \r
136         final int CV_NN = 0;\r
137         final int CV_SS = 1;\r
138         final int CV_SV = 2;\r
139         final int CV_WF = 3;\r
140         final int CV_EE = 4;\r
141         \r
142         /* DQtype */\r
143         final int CV_CENTERED = 1;\r
144         final int CV_FORWARD = 2;\r
145         \r
146         /* \r
147          * ----------------------------------------\r
148          * CVODES return flags\r
149          * ----------------------------------------\r
150          */\r
151 \r
152         final int CV_SUCCESS = 0;\r
153         final int CV_TSTOP_RETURN = 1;\r
154         final int CV_ROOT_RETURN = 2;\r
155 \r
156         final int CV_WARNING = 99;\r
157 \r
158         final int CV_TOO_MUCH_WORK = -1;\r
159         final int CV_TOO_MUCH_ACC = -2;\r
160         final int CV_ERR_FAILURE = -3;\r
161         final int CV_CONV_FAILURE = -4;\r
162 \r
163         final int CV_LINIT_FAIL = -5;\r
164         final int CV_LSETUP_FAIL = -6;\r
165         final int CV_LSOLVE_FAIL = -7;\r
166         final int CV_RHSFUNC_FAIL = -8;\r
167         final int CV_FIRST_RHSFUNC_ERR = -9;\r
168         final int CV_REPTD_RHSFUNC_ERR = -10;\r
169         final int CV_UNREC_RHSFUNC_ERR = -11;\r
170         final int  CV_RTFUNC_FAIL = -12;\r
171 \r
172         final int CV_MEM_FAIL = -20;\r
173         final int CV_MEM_NULL = -21;\r
174         final int CV_ILL_INPUT = -22;\r
175         final int CV_NO_MALLOC = -23;\r
176         final int CV_BAD_K = -24;\r
177         final int CV_BAD_T = -25;\r
178         final int CV_BAD_DKY = -26;\r
179         final int CV_TOO_CLOSE = -27;\r
180 \r
181         final int CV_NO_QUAD = -30;\r
182         final int CV_QRHSFUNC_FAIL = -31;\r
183         final int CV_FIRST_QRHSFUNC_ERR = -32;\r
184         final int CV_REPTD_QRHSFUNC_ERR = -33;\r
185         final int CV_UNREC_QRHSFUNC_ERR = -34;\r
186 \r
187         final int CV_NO_SENS = -40;\r
188         final int CV_SRHSFUNC_FAIL = -41;\r
189         final int CV_FIRST_SRHSFUNC_ERR = -42;\r
190         final int CV_REPTD_SRHSFUNC_ERR = -43;\r
191         final int CV_UNREC_SRHSFUNC_ERR = -44;\r
192 \r
193         final int CV_BAD_IS = -45;\r
194 \r
195         final int CV_NO_QUADSENS = -50;\r
196         final int CV_QSRHSFUNC_FAIL = -51;\r
197         final int CV_FIRST_QSRHSFUNC_ERR = -52;\r
198         final int CV_REPTD_QSRHSFUNC_ERR = -53;\r
199         final int CV_UNREC_QSRHSFUNC_ERR = -54;\r
200 \r
201         final double HMIN_DEFAULT = 0.0;    /* hmin default value     */\r
202         final double HMAX_INV_DEFAULT = 0.0;    /* hmax_inv default value */\r
203         final int MXHNIL_DEFAULT =  10;             /* mxhnil default value   */\r
204         final long MXSTEP_DEFAULT = 500L;            /* mxstep default value   */\r
205         final int NLS_MAXCOR = 3;\r
206         final int MXNCF = 10;\r
207         final int MXNEF = 7;\r
208         final double CORTES = 0.1;\r
209         final double ZERO = 0.0;\r
210         final double TINY = 1.0E-10;\r
211         final double PT1 = 0.1;\r
212         final double POINT2 = 0.2;\r
213         final double FOURTH = 0.25;\r
214         final double HALF = 0.5;\r
215         final double H_BIAS = HALF;\r
216         final double ONE = 1.0;\r
217         final double TWO = 2.0;\r
218         final double THREE = 3.0;\r
219         final double FOUR = 4.0;\r
220         final double FIVE = 5.0;\r
221         final double TWELVE = 12.0;\r
222         final double HUNDRED = 100.0;\r
223         final double ETAMXF = 0.2;\r
224         final double ETAMX1 = 10000.0;\r
225         final double ETAMX2 = 10.0;\r
226         final double ETAMX3 = 10.0;\r
227         final double HUB_FACTOR = 0.1;\r
228         final double HLB_FACTOR = 100.0;\r
229         final double FUZZ_FACTOR = 100.0;\r
230         final int MAX_ITERS = 4;\r
231         // CRDOWN constant used in the estimation of the convergence rate (crate)\r
232         //        of the iterates for the nonlinear equation\r
233         final double CRDOWN = 0.3;\r
234         // DGMAX  iter == CV_NEWTON, |gamma/gammap-1| > DGMAX => call lsetup\r
235     final double DGMAX = 0.3;\r
236     // RDIV declare divergence if ratio del/delp > RDIV\r
237     final double RDIV = TWO;\r
238     // MSBP  max no. of steps between lsetup calls\r
239     final int MSBP = 20;\r
240     // ONEPSM (1+epsilon) used in testing if the step size is below its bound\r
241     final double ONEPSM = 1.000001;\r
242     // THRESH      if eta < THRESH reject a change in step size or order\r
243     final double THRESH = 1.5;\r
244     // SMALL_NST   nst > SMALL_NST => use ETAMX3 \r
245     final int SMALL_NST = 10;\r
246     final double ETAMIN = 0.1;\r
247     // MXNEF1      max no. of error test failures before forcing a reduction of order\r
248     // SMALL_NEF   if an error failure occurs and SMALL_NEF <= nef <= MXNEF1, then\r
249     //             reset eta =  SUNMIN(eta, ETAMXF)\r
250     final double MXNEF1 = 3;\r
251     final double SMALL_NEF = 2;\r
252 \r
253     final double ETACF = 0.25;\r
254     final double ADDON  = 0.000001;\r
255     final double BIAS1 = 6.0;\r
256     final double BIAS2 = 6.0;\r
257     final double BIAS3 = 10.0;\r
258     // LONG_WAIT   number of steps to wait before considering an order change when\r
259     //             q==1 and MXNEF1 error test failures have occurred\r
260     final int LONG_WAIT = 10;\r
261 \r
262 \r
263         final int RTFOUND = +1;\r
264         final int CLOSERT = +3;\r
265 \r
266         /*\r
267          * Control constants for sensitivity DQ\r
268          * ------------------------------------\r
269          */\r
270 \r
271         final int CENTERED1 = +1;\r
272         final int CENTERED2 = +2;\r
273         final int FORWARD1 = +3;\r
274         final int FORWARD2 = +4;\r
275 \r
276         \r
277         final String MSG_TIME = "t = %lg";\r
278         final String MSG_TIME_H = "t = %lg and h = %lg";\r
279         final String MSG_TIME_INT = "t = %lg is not between tcur - hu = %lg and tcur = %lg.";\r
280         final String MSG_TIME_TOUT = "tout = %lg";\r
281         final String MSG_TIME_TSTOP = "tstop = %lg";\r
282 \r
283         \r
284         /* Initialization and I/O error messages */\r
285 \r
286         final String MSGCV_NO_MEM = "cvode_mem = NULL illegal.";\r
287         final String MSGCV_CVMEM_FAIL = "Allocation of cvode_mem failed.";\r
288         final String MSGCV_MEM_FAIL = "A memory request failed.";\r
289         final String MSGCV_BAD_LMM = "Illegal value for lmm. The legal values are CV_ADAMS and CV_BDF.";\r
290         final String MSGCV_BAD_ITER = "Illegal value for iter. The legal values are CV_FUNCTIONAL and CV_NEWTON.";\r
291         final String MSGCV_NO_MALLOC = "Attempt to call before CVodeInit.";\r
292         final String MSGCV_NEG_MAXORD = "maxord <= 0 illegal.";\r
293         final String MSGCV_BAD_MAXORD = "Illegal attempt to increase maximum method order.";\r
294         final String MSGCV_SET_SLDET = "Attempt to use stability limit detection with the CV_ADAMS method illegal.";\r
295         final String MSGCV_NEG_HMIN = "hmin < 0 illegal.";\r
296         final String MSGCV_NEG_HMAX = "hmax < 0 illegal.";\r
297         final String MSGCV_BAD_HMIN_HMAX = "Inconsistent step size limits: hmin > hmax.";\r
298         final String MSGCV_BAD_RELTOL = "reltol < 0 illegal.";\r
299         final String MSGCV_BAD_ABSTOL = "abstol has negative component(s) (illegal).";\r
300         final String MSGCV_NULL_ABSTOL = "abstol = NULL illegal.";\r
301         final String MSGCV_NULL_Y0 = "y0 = NULL illegal.";\r
302         final String MSGCV_NULL_F = "f = NULL illegal.";\r
303         final String MSGCV_NULL_G = "g = NULL illegal.";\r
304         final String MSGCV_BAD_NVECTOR = "A required vector operation is not implemented.";\r
305         final String MSGCV_BAD_K = "Illegal value for k.";\r
306         final String MSGCV_NULL_DKY = "dky = NULL illegal.";\r
307         final String MSGCV_BAD_T = "Illegal value for t.";\r
308         final String MSGCV_NO_ROOT = "Rootfinding was not initialized.";\r
309 \r
310         final String MSGCV_NO_QUAD  = "Quadrature integration not activated.";\r
311         final String MSGCV_BAD_ITOLQ = "Illegal value for itolQ. The legal values are CV_SS and CV_SV.";\r
312         final String MSGCV_NULL_ABSTOLQ = "abstolQ = NULL illegal.";\r
313         final String MSGCV_BAD_RELTOLQ = "reltolQ < 0 illegal.";\r
314         final String MSGCV_BAD_ABSTOLQ = "abstolQ has negative component(s) (illegal).";  \r
315 \r
316         final String MSGCV_SENSINIT_2 = "Sensitivity analysis already initialized.";\r
317         final String MSGCV_NO_SENSI  = "Forward sensitivity analysis not activated.";\r
318         final String MSGCV_BAD_ITOLS = "Illegal value for itolS. The legal values are CV_SS, CV_SV, and CV_EE.";\r
319         final String MSGCV_NULL_ABSTOLS = "abstolS = NULL illegal.";\r
320         final String MSGCV_BAD_RELTOLS = "reltolS < 0 illegal.";\r
321         final String MSGCV_BAD_ABSTOLS = "abstolS has negative component(s) (illegal).";  \r
322         final String MSGCV_BAD_PBAR = "pbar has zero component(s) (illegal).";\r
323         final String MSGCV_BAD_PLIST = "plist has negative component(s) (illegal).";\r
324         final String MSGCV_BAD_NS = "NS <= 0 illegal.";\r
325         final String MSGCV_NULL_YS0 = "yS0 = NULL illegal.";\r
326         final String MSGCV_BAD_ISM = "Illegal value for ism. Legal values are: CV_SIMULTANEOUS, CV_STAGGERED and CV_STAGGERED1.";\r
327         final String MSGCV_BAD_IFS = "Illegal value for ifS. Legal values are: CV_ALLSENS and CV_ONESENS.";\r
328         final String MSGCV_BAD_ISM_IFS = "Illegal ism = CV_STAGGERED1 for CVodeSensInit.";\r
329         final String MSGCV_BAD_IS = "Illegal value for is.";\r
330         final String MSGCV_NULL_DKYA = "dkyA = NULL illegal.";\r
331         final String MSGCV_BAD_DQTYPE = "Illegal value for DQtype. Legal values are: CV_CENTERED and CV_FORWARD.";\r
332         final String MSGCV_BAD_DQRHO = "DQrhomax < 0 illegal.";\r
333 \r
334         final String MSGCV_BAD_ITOLQS = "Illegal value for itolQS. The legal values are CV_SS, CV_SV, and CV_EE.";\r
335         final String MSGCV_NULL_ABSTOLQS = "abstolQS = NULL illegal.";\r
336         final String MSGCV_BAD_RELTOLQS = "reltolQS < 0 illegal.";\r
337         final String MSGCV_BAD_ABSTOLQS = "abstolQS has negative component(s) (illegal).";  \r
338         final String MSGCV_NO_QUADSENSI = "Forward sensitivity analysis for quadrature variables not activated.";\r
339         final String MSGCV_NULL_YQS0 = "yQS0 = NULL illegal.";\r
340 \r
341         /* CVode Error Messages */\r
342 \r
343         final String MSGCV_NO_TOL = "No integration tolerances have been specified.";\r
344         final String MSGCV_LSOLVE_NULL = "The linear solver's solve routine is NULL.";\r
345         final String MSGCV_YOUT_NULL = "yout = NULL illegal.";\r
346         final String MSGCV_TRET_NULL = "tret = NULL illegal.";\r
347         final String MSGCV_BAD_EWT = "Initial ewt has component(s) equal to zero (illegal).";\r
348         final String MSGCV_EWT_NOW_BAD = "At " + MSG_TIME + ", a component of ewt has become <= 0.";\r
349         final String MSGCV_BAD_ITASK = "Illegal value for itask.";\r
350         final String MSGCV_BAD_H0 = "h0 and tout - t0 inconsistent.";\r
351         final String MSGCV_BAD_TOUT = "Trouble interpolating at " + MSG_TIME_TOUT + ". tout too far back in direction of integration";\r
352         final String MSGCV_EWT_FAIL = "The user-provide EwtSet function failed.";\r
353         final String MSGCV_EWT_NOW_FAIL = "At " + MSG_TIME + ", the user-provide EwtSet function failed.";\r
354         final String MSGCV_LINIT_FAIL = "The linear solver's init routine failed.";\r
355         final String MSGCV_HNIL_DONE = "The above warning has been issued mxhnil times and will not be issued again for this problem.";\r
356         final String MSGCV_TOO_CLOSE = "tout too close to t0 to start integration.";\r
357         final String MSGCV_MAX_STEPS = "At " + MSG_TIME + ", mxstep steps taken before reaching tout.";\r
358         final String MSGCV_TOO_MUCH_ACC = "At " + MSG_TIME + ", too much accuracy requested.";\r
359         final String MSGCV_HNIL = "Internal " + MSG_TIME_H + " are such that t + h = t on the next step. The solver will continue anyway.";\r
360         final String MSGCV_ERR_FAILS = "At " + MSG_TIME_H + ", the error test failed repeatedly or with |h| = hmin.";\r
361         final String MSGCV_CONV_FAILS = "At " + MSG_TIME_H + ", the corrector convergence test failed repeatedly or with |h| = hmin.";\r
362         final String MSGCV_SETUP_FAILED = "At " + MSG_TIME + ", the setup routine failed in an unrecoverable manner.";\r
363         final String MSGCV_SOLVE_FAILED = "At " + MSG_TIME + ", the solve routine failed in an unrecoverable manner.";\r
364         final String MSGCV_RHSFUNC_FAILED = "At " + MSG_TIME + ", the right-hand side routine failed in an unrecoverable manner.";\r
365         final String MSGCV_RHSFUNC_UNREC = "At " + MSG_TIME + ", the right-hand side failed in a recoverable manner, but no recovery is possible.";\r
366         final String MSGCV_RHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable right-hand side function errors.";\r
367         final String MSGCV_RHSFUNC_FIRST = "The right-hand side routine failed at the first call.";\r
368         final String MSGCV_RTFUNC_FAILED = "At " + MSG_TIME + ", the rootfinding routine failed in an unrecoverable manner.";\r
369         final String MSGCV_CLOSE_ROOTS = "Root found at and very near " + MSG_TIME + ".";\r
370         final String MSGCV_BAD_TSTOP = "The value " + MSG_TIME_TSTOP + " is behind current " + MSG_TIME + " in the direction of integration.";\r
371         final String MSGCV_INACTIVE_ROOTS = "At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.";\r
372 \r
373         final String MSGCV_NO_TOLQ = "No integration tolerances for quadrature variables have been specified.";\r
374         final String MSGCV_BAD_EWTQ = "Initial ewtQ has component(s) equal to zero (illegal).";\r
375         final String MSGCV_EWTQ_NOW_BAD = "At " + MSG_TIME + ", a component of ewtQ has become <= 0.";\r
376         final String MSGCV_QRHSFUNC_FAILED = "At " + MSG_TIME + ", the quadrature right-hand side routine failed in an unrecoverable manner.";\r
377         final String MSGCV_QRHSFUNC_UNREC = "At " + MSG_TIME + ", the quadrature right-hand side failed in a recoverable manner, but no recovery is possible.";\r
378         final String MSGCV_QRHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable quadrature right-hand side function errors.";\r
379         final String MSGCV_QRHSFUNC_FIRST = "The quadrature right-hand side routine failed at the first call.";\r
380 \r
381         final String MSGCV_NO_TOLS = "No integration tolerances for sensitivity variables have been specified.";\r
382         final String MSGCV_NULL_P = "p = NULL when using internal DQ for sensitivity RHS illegal.";\r
383         final String MSGCV_BAD_EWTS = "Initial ewtS has component(s) equal to zero (illegal).";\r
384         final String MSGCV_EWTS_NOW_BAD = "At " + MSG_TIME + ", a component of ewtS has become <= 0.";\r
385         final String MSGCV_SRHSFUNC_FAILED = "At " + MSG_TIME + ", the sensitivity right-hand side routine failed in an unrecoverable manner.";\r
386         final String MSGCV_SRHSFUNC_UNREC = "At " + MSG_TIME + ", the sensitivity right-hand side failed in a recoverable manner, but no recovery is possible.";\r
387         final String MSGCV_SRHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable sensitivity right-hand side function errors.";\r
388         final String MSGCV_SRHSFUNC_FIRST = "The sensitivity right-hand side routine failed at the first call.";\r
389 \r
390         final String MSGCV_NULL_FQ = "CVODES is expected to use DQ to evaluate the RHS of quad. sensi., but quadratures were not initialized.";\r
391         final String MSGCV_NO_TOLQS = "No integration tolerances for quadrature sensitivity variables have been specified.";\r
392         final String MSGCV_BAD_EWTQS = "Initial ewtQS has component(s) equal to zero (illegal).";\r
393         final String MSGCV_EWTQS_NOW_BAD = "At " + MSG_TIME + ", a component of ewtQS has become <= 0.";\r
394         final String MSGCV_QSRHSFUNC_FAILED = "At " + MSG_TIME + ", the quadrature sensitivity right-hand side routine failed in an unrecoverable manner.";\r
395         final String MSGCV_QSRHSFUNC_UNREC = "At " + MSG_TIME + ", the quadrature sensitivity right-hand side failed in a recoverable manner, but no recovery is possible.";\r
396         final String MSGCV_QSRHSFUNC_REPTD = "At " + MSG_TIME + " repeated recoverable quadrature sensitivity right-hand side function errors.";\r
397         final String MSGCV_QSRHSFUNC_FIRST = "The quadrature sensitivity right-hand side routine failed at the first call.";\r
398 \r
399         /* \r
400          * =================================================================\r
401          *   C V O D E A    E R R O R    M E S S A G E S\r
402          * =================================================================\r
403          */\r
404 \r
405         final String MSGCV_NO_ADJ = "Illegal attempt to call before calling CVodeAdjMalloc.";\r
406         final String MSGCV_BAD_STEPS = "Steps nonpositive illegal.";\r
407         final String MSGCV_BAD_INTERP = "Illegal value for interp.";\r
408         final String MSGCV_BAD_WHICH  = "Illegal value for which.";\r
409         final String MSGCV_NO_BCK = "No backward problems have been defined yet.";\r
410         final String MSGCV_NO_FWD = "Illegal attempt to call before calling CVodeF.";\r
411         final String MSGCV_BAD_TB0 = "The initial time tB0 for problem %d is outside the interval over which the forward problem was solved.";\r
412         final String MSGCV_BAD_SENSI = "At least one backward problem requires sensitivities, but they were not stored for interpolation.";\r
413         final String MSGCV_BAD_ITASKB = "Illegal value for itaskB. Legal values are CV_NORMAL and CV_ONE_STEP.";\r
414         final String MSGCV_BAD_TBOUT  = "The final time tBout is outside the interval over which the forward problem was solved.";\r
415         final String MSGCV_BACK_ERROR  = "Error occured while integrating backward problem # %d"; \r
416         final String MSGCV_BAD_TINTERP = "Bad t = %g for interpolation.";\r
417         final String MSGCV_WRONG_INTERP = "This function cannot be called for the specified interp type.";\r
418         \r
419         final int CVDLS_SUCCESS = 0;\r
420         final int CVDLS_MEM_NULL = -1;\r
421         final int CVDLS_LMEM_NULL = -2;\r
422         final int CVDLS_ILL_INPUT = -3;\r
423         final int  CVDLS_MEM_FAIL = -4;\r
424 \r
425         /* Additional last_flag values */\r
426 \r
427         final int CVDLS_JACFUNC_UNRECVR = -5;\r
428         final int CVDLS_JACFUNC_RECVR = -6;\r
429         final int CVDLS_SUNMAT_FAIL = -7;\r
430 \r
431         /* Return values for the adjoint module */\r
432 \r
433         final int CVDLS_NO_ADJ = -101;\r
434         final int CVDLS_LMEMB_NULL = -102;\r
435 \r
436         final String MSGD_CVMEM_NULL = "Integrator memory is NULL.";\r
437         final String MSGD_BAD_NVECTOR = "A required vector operation is not implemented.";\r
438         final String MSGD_BAD_SIZES = "Illegal bandwidth parameter(s). Must have 0 <=  ml, mu <= N-1.";\r
439         final String MSGD_MEM_FAIL = "A memory request failed.";\r
440         final String MSGD_LMEM_NULL ="Linear solver memory is NULL.";\r
441         final String MSGD_JACFUNC_FAILED = "The Jacobian routine failed in an unrecoverable manner.";\r
442         final String MSGD_MATCOPY_FAILED = "The SUNMatCopy routine failed in an unrecoverable manner.";\r
443         final String MSGD_MATZERO_FAILED = "The SUNMatZero routine failed in an unrecoverable manner.";\r
444         final String MSGD_MATSCALEADDI_FAILED = "The SUNMatScaleAddI routine failed in an unrecoverable manner.";\r
445 \r
446         final String MSGD_NO_ADJ = "Illegal attempt to call before calling CVodeAdjMalloc.";\r
447         final String MSGD_BAD_WHICH = "Illegal value for which.";\r
448         final String MSGD_LMEMB_NULL = "Linear solver memory is NULL for the backward integration.";\r
449         final String MSGD_BAD_TINTERP = "Bad t for interpolation.";\r
450         \r
451         final int DO_ERROR_TEST = +2;\r
452         final int PREDICT_AGAIN = +3;\r
453 \r
454         final int CONV_FAIL = +4; \r
455         final int TRY_AGAIN = +5;\r
456 \r
457         final int FIRST_CALL = +6;\r
458         final int PREV_CONV_FAIL = +7;\r
459         final int PREV_ERR_FAIL = +8;\r
460 \r
461         final int RHSFUNC_RECVR = +9;\r
462 \r
463         final int QRHSFUNC_RECVR = +11;\r
464         final int SRHSFUNC_RECVR = +12;\r
465         final int QSRHSFUNC_RECVR = +13;\r
466 \r
467         /*\r
468          * -----------------------------------------------------------------\r
469          * Communication between CVODE and a CVODE Linear Solver\r
470          * -----------------------------------------------------------------\r
471          * convfail (input to cv_lsetup)\r
472          *\r
473          * CV_NO_FAILURES : Either this is the first cv_setup call for this\r
474          *                  step, or the local error test failed on the\r
475          *                  previous attempt at this step (but the Newton\r
476          *                  iteration converged).\r
477          *\r
478          * CV_FAIL_BAD_J  : This value is passed to cv_lsetup if\r
479          *\r
480          *                  (a) The previous Newton corrector iteration\r
481          *                      did not converge and the linear solver's\r
482          *                      setup routine indicated that its Jacobian-\r
483          *                      related data is not current\r
484          *                                   or\r
485          *                  (b) During the previous Newton corrector\r
486          *                      iteration, the linear solver's solve routine\r
487          *                      failed in a recoverable manner and the\r
488          *                      linear solver's setup routine indicated that\r
489          *                      its Jacobian-related data is not current.\r
490          *\r
491          * CV_FAIL_OTHER  : During the current internal step try, the\r
492          *                  previous Newton iteration failed to converge\r
493          *                  even though the linear solver was using current\r
494          *                  Jacobian-related data.\r
495          * -----------------------------------------------------------------\r
496          */\r
497         /* Constants for convfail (input to cv_lsetup) */\r
498         final int CV_NO_FAILURES = 0;\r
499     final int CV_FAIL_BAD_J = 1;\r
500         final int CV_FAIL_OTHER = 2;\r
501         \r
502         /*-----------------------------------------------------------------\r
503           CVDLS solver constants\r
504           -----------------------------------------------------------------\r
505           CVD_MSBJ   maximum number of steps between Jacobian evaluations\r
506           CVD_DGMAX  maximum change in gamma between Jacobian evaluations\r
507           -----------------------------------------------------------------*/\r
508 \r
509         final int CVD_MSBJ = 50;\r
510         final double CVD_DGMAX = 0.2;\r
511 \r
512 \r
513         /*-----------------------------------------------------------------\r
514          * IV. SUNLinearSolver error codes\r
515          * ---------------------------------------------------------------\r
516          */\r
517 \r
518         final int SUNLS_SUCCESS = 0;  /* successful/converged          */\r
519 \r
520         final int SUNLS_MEM_NULL = -1;  /* mem argument is NULL          */\r
521         final int SUNLS_ILL_INPUT = -2;  /* illegal function input        */\r
522         final int SUNLS_MEM_FAIL = -3;  /* failed memory access          */\r
523         final int SUNLS_ATIMES_FAIL_UNREC = -4;  /* atimes unrecoverable failure  */\r
524         final int SUNLS_PSET_FAIL_UNREC = -5;  /* pset unrecoverable failure    */\r
525         final int SUNLS_PSOLVE_FAIL_UNREC = -6;  /* psolve unrecoverable failure  */\r
526         final int SUNLS_PACKAGE_FAIL_UNREC = -7;  /* external package unrec. fail  */\r
527         final int SUNLS_GS_FAIL = -8;  /* Gram-Schmidt failure          */        \r
528         final int SUNLS_QRSOL_FAIL = -9;  /* QRsol found singular R        */\r
529 \r
530         final int SUNLS_RES_REDUCED = 1;  /* nonconv. solve, resid reduced */\r
531         final int SUNLS_CONV_FAIL = 2;  /* nonconvergent solve           */\r
532         final int SUNLS_ATIMES_FAIL_REC = 3;  /* atimes failed recoverably     */\r
533         final int SUNLS_PSET_FAIL_REC = 4;  /* pset failed recoverably       */\r
534         final int SUNLS_PSOLVE_FAIL_REC = 5;  /* psolve failed recoverably     */\r
535         final int SUNLS_PACKAGE_FAIL_REC = 6;  /* external package recov. fail  */\r
536         final int SUNLS_QRFACT_FAIL = 7;  /* QRfact found singular matrix  */\r
537         final int SUNLS_LUFACT_FAIL = 8;  /* LUfact found singular matrix  */\r
538 \r
539         \r
540         public enum SUNLinearSolver_Type{\r
541                   SUNLINEARSOLVER_DIRECT,\r
542                   SUNLINEARSOLVER_ITERATIVE,\r
543                   SUNLINEARSOLVER_CUSTOM\r
544                 };\r
545 \r
546     final int cvDlsDQJac = -1;\r
547     \r
548     // For cv_mem.cv_efun_select\r
549     final int cvEwtSet_select = 1;\r
550     // For cv_mem.cv_linit_select\r
551     final int cvDlsInitialize_select = 1;\r
552     // For cv_mem.cv_lsolve_select\r
553     final int cvDlsSolve_select = 1;\r
554     // For cv_mem.cv_lfree_select\r
555     final int cvDlsFree_select = 1;\r
556     // For cv_mem.cv_setup_select\r
557     final int cvDlsSetup_select = 1;\r
558 \r
559     final int cvsRoberts_dns = 1;\r
560     int problem = cvsRoberts_dns;\r
561         \r
562         \r
563         public CVODES() {\r
564                 // eps returns the distance from 1.0 to the next larger double-precision\r
565                 // number, that is, eps = 2^-52.\r
566                 // epsilon = D1MACH(4)\r
567                 // Machine epsilon is the smallest positive epsilon such that\r
568                 // (1.0 + epsilon) != 1.0.\r
569                 // epsilon = 2**(1 - doubleDigits) = 2**(1 - 53) = 2**(-52)\r
570                 // epsilon = 2.2204460e-16\r
571                 // epsilon is called the largest relative spacing\r
572                 DBL_EPSILON = 1.0;\r
573                 double neweps = 1.0;\r
574 \r
575                 while (true) {\r
576 \r
577                         if (1.0 == (1.0 + neweps)) {\r
578                                 break;\r
579                         } else {\r
580                                 DBL_EPSILON = neweps;\r
581                                 neweps = neweps / 2.0;\r
582                         }\r
583                 } // while(true)\r
584                 \r
585             UNIT_ROUNDOFF = DBL_EPSILON;\r
586             if (problem == cvsRoberts_dns) {\r
587                 runcvsRoberts_dns();\r
588             }\r
589             return;\r
590         }\r
591         \r
592         private void runcvsRoberts_dns() {\r
593                 /* -----------------------------------------------------------------\r
594                  * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and\r
595                  *                Radu Serban @ LLNL\r
596                  * -----------------------------------------------------------------\r
597                  * Example problem:\r
598                  * \r
599                  * The following is a simple example problem, with the coding\r
600                  * needed for its solution by CVODE. The problem is from\r
601                  * chemical kinetics, and consists of the following three rate\r
602                  * equations:         \r
603                  *    dy1/dt = -.04*y1 + 1.e4*y2*y3\r
604                  *    dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2\r
605                  *    dy3/dt = 3.e7*(y2)^2\r
606                  * on the interval from t = 0.0 to t = 4.e10, with initial\r
607                  * conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.\r
608                  * While integrating the system, we also use the rootfinding\r
609                  * feature to find the points at which y1 = 1e-4 or at which\r
610                  * y3 = 0.01. This program solves the problem with the BDF method,\r
611                  * Newton iteration with the SUNDENSE dense linear solver, and a\r
612                  * user-supplied Jacobian routine.\r
613                  * It uses a scalar relative tolerance and a vector absolute\r
614                  * tolerance. Output is printed in decades from t = .4 to t = 4.e10.\r
615                  * Run statistics (optional outputs) are printed at the end.\r
616                  * -----------------------------------------------------------------*/\r
617                 \r
618                 /** Problem Constants */\r
619                 final int NEQ = 3; // Number of equations\r
620                 final double Y1 = 1.0; // Initial y components\r
621                 final double Y2 = 0.0;\r
622                 final double Y3 = 0.0;\r
623                 final double RTOL = 1.0E-4; // scalar relative tolerance\r
624                 final double ATOL1 = 1.0E-8; // vector absolute tolerance components\r
625                 final double ATOL2 = 1.0E-14;\r
626                 final double ATOL3 = 1.0E-6;\r
627                 final double T0 = 0.0; // initial time\r
628                 final double T1 = 0.4; // first output time\r
629                 final double TMULT = 10.0; // output time factor\r
630                 final int NOUT = 12; // number of output times\r
631                 int f = cvsRoberts_dns;\r
632                 int g = cvsRoberts_dns;\r
633                 int Jac = cvsRoberts_dns;\r
634                 \r
635                 // Set initial conditions\r
636                 NVector y = new NVector();\r
637                 NVector abstol = new NVector();\r
638                 double yr[] = new double[]{Y1,Y2,Y3};\r
639                 N_VNew_Serial(y, NEQ);\r
640                 y.setData(yr);\r
641                 double reltol = RTOL; // Set the scalar relative tolerance\r
642                 // Set the vector absolute tolerance\r
643                 double abstolr[] = new double[]{ATOL1,ATOL2,ATOL3};\r
644                 N_VNew_Serial(abstol,NEQ);\r
645                 abstol.setData(abstolr);\r
646                 CVodeMemRec cvode_mem;\r
647                 int flag;\r
648                 int flagr;\r
649                 double A[][];\r
650                 SUNLinearSolver LS;\r
651                 double tout;\r
652                 int iout;\r
653                 double t[] = new double[1];\r
654                 int rootsfound[] = new int[2];\r
655                 \r
656                 // Call CVodeCreate to create the solver memory and specify the\r
657                 // Backward Differentiation Formula and the use of a Newton\r
658                 // iteration\r
659                 cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);\r
660                 if (cvode_mem == null) {\r
661                     return;     \r
662                 }\r
663                 \r
664                 // Call CVodeInit to initialize the integrator memory and specify the\r
665                 // user's right hand side function in y' = f(t,y), the initial time T0, and\r
666             // the initial dependent variable vector y\r
667                 flag = CVodeInit(cvode_mem, f, T0, y);\r
668                 if (flag != CV_SUCCESS) {\r
669                         return;\r
670                 }\r
671                 \r
672                 // Call CVodeSVtolerances to specify the scalar relative tolerance\r
673                 // and vector absolute tolerances\r
674                 flag = CVodeSVtolerances(cvode_mem, reltol, abstol);\r
675                 if (flag != CV_SUCCESS) {\r
676                         return;\r
677                 }\r
678                 \r
679                 // Call CVRootInit to specify the root function g with 2 components\r
680                 flag = CVodeRootInit(cvode_mem, 2, g);\r
681                 if (flag != CV_SUCCESS) {\r
682                         return;\r
683                 }\r
684                 \r
685                 // Create dense SUNMATRIX for use in linear solver\r
686                 // indexed by columns stacked on top of each other\r
687                 try {\r
688                     A = new double[NEQ][NEQ];\r
689                 }\r
690                 catch (OutOfMemoryError e) {\r
691                     MipavUtil.displayError("Out of memory error trying to create A");\r
692                     return;\r
693                 }\r
694                 \r
695                 // Create dense SUNLinearSolver object for use by CVode\r
696                 LS = SUNDenseLinearSolver(y, A);\r
697                 if (LS == null) {\r
698                         return;\r
699                 }\r
700                 \r
701                 // Call CVDlsSetLinearSolver to attach the matrix and linear solver to CVode\r
702                 flag = CVDlsSetLinearSolver(cvode_mem, LS, A);\r
703                 if (flag != CVDLS_SUCCESS) {\r
704                         return;\r
705                 }\r
706                 \r
707                 // Set the user-supplied Jacobian routine Jac\r
708                 flag = CVDlsSetJacFn(cvode_mem, Jac);\r
709                 if (flag != CVDLS_SUCCESS) {\r
710                         return;\r
711                 }\r
712                 \r
713                 // In loop, call CVode, print results, and test for error.\r
714                 // Break out of loop when NOUT preset output times have been reached.\r
715                 System.out.println(" \n3-species kinetics problem\n");\r
716                 \r
717                 iout = 0;\r
718                 tout = T1;\r
719                 \r
720                 while (true) {\r
721                         flag = CVode(cvode_mem, tout, y, t, CV_NORMAL);\r
722                         System.out.println("At t = " + t + " y[0] = " + y.data[0] + " y[1] = " + y.data[1] + " y[2] = " + y.data[2]);\r
723                         \r
724                         if (flag == CV_ROOT_RETURN) {\r
725                             flagr = CVodeGetRootInfo(cvode_mem, rootsfound)     ;\r
726                             if (flagr == CV_MEM_NULL) {\r
727                                 return;\r
728                             }\r
729                             System.out.println("Roots found are " + rootsfound[0] + " and " + rootsfound[1]);\r
730                         }\r
731                         \r
732                         if (flag < 0) {\r
733                                 System.out.println("CVode failed with flag = " + flag);\r
734                                 break;\r
735                         }\r
736                         \r
737                         if (flag == CV_SUCCESS) {\r
738                                 iout++;\r
739                                 tout *= TMULT;\r
740                         }\r
741                         \r
742                         if (iout == NOUT) {\r
743                                 break;\r
744                         }\r
745                 } // while (true)\r
746                 \r
747                 // Print some final statistics\r
748                 PrintFinalStats(cvode_mem);\r
749         }\r
750         \r
751         private void PrintFinalStats(CVodeMemRec cv_mem) {\r
752             System.out.println("Number of integrations steps = " + cv_mem.cv_nst);\r
753             System.out.println("Number of calls to f = " + cv_mem.cv_nfe);\r
754             System.out.println("Number of calls to the linear solver setup routine = " + cv_mem.cv_nsetups);\r
755         }\r
756         \r
757         private void N_VNew_Serial(NVector y, int length) {\r
758             if ((y != null) && (length > 0)) {\r
759                 y.data = new double[length];\r
760                 y.own_data = true;\r
761             }\r
762         }\r
763         \r
764         /**\r
765          * f routine.  Compte function f(t,y).\r
766          * @param t\r
767          * @param yv\r
768          * @param ydotv\r
769          * @param user_data\r
770          * @return\r
771          */\r
772         private int f(double t, NVector yv, NVector ydotv, CVodeMemRec user_data) {\r
773                 double y[] = yv.getData();\r
774                 double ydot[] = ydotv.getData();\r
775                 if (problem == 1) {\r
776                     ydot[0] = -0.04*y[0] + 1.0E4*y[1]*y[2];\r
777                     ydot[2] = 3.0E7*y[1]*y[1];\r
778                     ydot[1] = -ydot[0] - ydot[2];\r
779                 }\r
780                 return 0;\r
781         }\r
782         \r
783         private int fQ(double t, NVector x, NVector y, CVodeMemRec user_data) {\r
784                 return 0;\r
785         }\r
786         \r
787         /**\r
788          * g routine.  Compute functions g_i(t,y) for i = 0,1.\r
789          * @param t\r
790          * @param yv\r
791          * @param gout\r
792          * @param user_data\r
793          * @return\r
794          */\r
795         private int g(double t, NVector yv, double gout[], CVodeMemRec user_data) {\r
796                 double y[] = yv.getData();\r
797                 if (problem == 1) {\r
798                     gout[0] = y[0] - 0.0001;\r
799                     gout[1] = y[2] - 0.01;\r
800                 }\r
801                 return 0;\r
802         }\r
803         \r
804         /**\r
805          * Jacobian routine.  Compute J(t,y) = df/dy.\r
806          * @param t\r
807          * @param y\r
808          * @param fy\r
809          * @param J\r
810          * @param tmp1\r
811          * @param tmp2\r
812          * @return\r
813          */\r
814         private int Jac(double t, NVector yv, NVector fy, double J[][], CVodeMemRec data, NVector tmp1, NVector tmp2, NVector tmp3) {\r
815                 double y[] = yv.getData();\r
816             J[0][0] = -0.04;\r
817             J[0][1] = 1.0E4 * y[2];\r
818             J[0][2] = 1.0E4 * y[1];\r
819             \r
820             J[1][0] = 0.04;\r
821             J[1][1] = -1.0E4 * y[2] - 6.0E7*y[1];\r
822             J[1][2] = -1.0E4 * y[1];\r
823             \r
824             J[2][0] = ZERO;\r
825             J[2][1] = 6.0E7 * y[1];\r
826             J[2][2] = ZERO;\r
827             \r
828             return 0;\r
829         }\r
830         \r
831         \r
832         \r
833         \r
834         \r
835         // Types: struct CVodeMemRec, CVodeMem\r
836         // -----------------------------------------------------------------\r
837         // The type CVodeMem is type pointer to struct CVodeMemRec.\r
838         // This structure contains fields to keep track of problem state.\r
839         \r
840         private class NVector {\r
841                 double data[];\r
842                 boolean own_data;\r
843                 \r
844                 void setData(double data[]) {\r
845                         this.data = data;\r
846                 }\r
847                 double[] getData() {\r
848                     return data;        \r
849                 }\r
850         }\r
851         \r
852           \r
853         private class CVodeMemRec {\r
854             \r
855           double cv_uround;         /* machine unit roundoff                         */   \r
856 \r
857           /*-------------------------- \r
858             Problem Specification Data \r
859             --------------------------*/\r
860 \r
861           //CVRhsFn cv_f;               /* y' = f(t,y(t))                                */\r
862           int cv_f;\r
863           //void *cv_user_data;         /* user pointer passed to f                      */\r
864           CVodeMemRec cv_user_data; \r
865 \r
866           int cv_lmm;                 /* lmm = ADAMS or BDF                            */\r
867           int cv_iter;                /* iter = FUNCTIONAL or NEWTON                   */\r
868 \r
869           int cv_itol;                /* itol = CV_SS, CV_SV, or CV_WF, or CV_NN       */\r
870           double cv_reltol;         /* relative tolerance                            */\r
871           double cv_Sabstol;        /* scalar absolute tolerance                     */\r
872           NVector cv_Vabstol = new NVector();        /* vector absolute tolerance                     */\r
873           boolean cv_user_efun;   /* SUNTRUE if user sets efun                     */\r
874           //CVEwtFn cv_efun;            /* function to set ewt                           */\r
875           int cv_efun_select;\r
876           //void *cv_e_data;            /* user pointer passed to efun                   */\r
877           CVodeMemRec cv_e_data;\r
878 \r
879           /*-----------------------\r
880             Quadrature Related Data \r
881             -----------------------*/\r
882 \r
883           boolean cv_quadr;       /* SUNTRUE if integrating quadratures            */\r
884 \r
885           //CVQuadRhsFn cv_fQ;          /* q' = fQ(t, y(t))                              */\r
886 \r
887           boolean cv_errconQ;     /* SUNTRUE if quadrs. are included in error test */\r
888 \r
889           int cv_itolQ;               /* itolQ = CV_SS or CV_SV                        */\r
890           double cv_reltolQ;        /* relative tolerance for quadratures            */\r
891           double cv_SabstolQ;       /* scalar absolute tolerance for quadratures     */\r
892           NVector cv_VabstolQ = new NVector();       /* vector absolute tolerance for quadratures     */\r
893 \r
894           /*------------------------\r
895             Sensitivity Related Data \r
896             ------------------------*/\r
897 \r
898           boolean cv_sensi;       /* SUNTRUE if computing sensitivities           */\r
899 \r
900           int cv_Ns;                  /* Number of sensitivities                      */\r
901 \r
902           int cv_ism;                 /* ism = SIMULTANEOUS or STAGGERED              */\r
903 \r
904           //CVSensRhsFn cv_fS;          /* fS = (df/dy)*yS + (df/dp)                    */\r
905           //CVSensRhs1Fn cv_fS1;        /* fS1 = (df/dy)*yS_i + (df/dp)                 */\r
906           //void *cv_fS_data;           /* data pointer passed to fS                    */\r
907           CVodeMemRec cv_fS_data;\r
908           boolean cv_fSDQ;        /* SUNTRUE if using internal DQ functions       */\r
909           int cv_ifS;                 /* ifS = ALLSENS or ONESENS                     */\r
910 \r
911           double cv_p[];             /* parameters in f(t,y,p)                       */\r
912           double cv_pbar[];          /* scale factors for parameters                 */\r
913           int cv_plist[];              /* list of sensitivities                        */\r
914           int cv_DQtype;              /* central/forward finite differences           */\r
915           double cv_DQrhomax;       /* cut-off value for separate/simultaneous FD   */\r
916 \r
917           boolean cv_errconS;     /* SUNTRUE if yS are considered in err. control */\r
918 \r
919           int cv_itolS;\r
920           double cv_reltolS;        /* relative tolerance for sensitivities         */\r
921           double cv_SabstolS[];      /* scalar absolute tolerances for sensi.        */\r
922           NVector cv_VabstolS[];      /* vector absolute tolerances for sensi.        */\r
923 \r
924           /*-----------------------------------\r
925             Quadrature Sensitivity Related Data \r
926             -----------------------------------*/\r
927 \r
928           boolean cv_quadr_sensi; /* SUNTRUE if computing sensitivties of quadrs. */\r
929 \r
930           //CVQuadSensRhsFn cv_fQS;     /* fQS = (dfQ/dy)*yS + (dfQ/dp)                 */\r
931           //void *cv_fQS_data;          /* data pointer passed to fQS                   */\r
932           CVodeMemRec cv_fQS_data;\r
933           boolean cv_fQSDQ;       /* SUNTRUE if using internal DQ functions       */\r
934 \r
935           boolean cv_errconQS;    /* SUNTRUE if yQS are considered in err. con.   */\r
936 \r
937           int cv_itolQS;\r
938           double cv_reltolQS;       /* relative tolerance for yQS                   */\r
939           double cv_SabstolQS[];     /* scalar absolute tolerances for yQS           */\r
940           NVector cv_VabstolQS[];     /* vector absolute tolerances for yQS           */\r
941 \r
942           /*-----------------------\r
943             Nordsieck History Array \r
944             -----------------------*/\r
945 \r
946           NVector cv_zn[] = new NVector[L_MAX];   /* Nordsieck array, of size N x (q+1).\r
947                                          zn[j] is a vector of length N (j=0,...,q)\r
948                                          zn[j] = [1/factorial(j)] * h^j * \r
949                                          (jth derivative of the interpolating poly.)  */\r
950 \r
951           /*-------------------\r
952             Vectors of length N \r
953             -------------------*/\r
954 \r
955           NVector cv_ewt = new NVector();            /* error weight vector                          */\r
956           NVector cv_y = new NVector();              // y is used as temporary storage by the solver.\r
957                                       /*   The memory is provided by the user to CVode \r
958                                          where the vector is named yout.              */\r
959           NVector cv_acor = new NVector();           // In the context of the solution of the\r
960                                        /*  nonlinear equation, acor = y_n(m) - y_n(0).\r
961                                          On return, this vector is scaled to give\r
962                                          the estimated local error in y.              */\r
963           NVector cv_tempv = new NVector();          /* temporary storage vector                     */\r
964           NVector cv_ftemp = new NVector();          /* temporary storage vector                     */\r
965           \r
966           /*--------------------------\r
967             Quadrature Related Vectors \r
968             --------------------------*/\r
969 \r
970           NVector cv_znQ[] = new NVector[L_MAX];     /* Nordsieck arrays for quadratures             */\r
971           NVector cv_ewtQ = new NVector();           /* error weight vector for quadratures          */\r
972           NVector cv_yQ = new NVector();             /* Unlike y, yQ is not allocated by the user    */\r
973           NVector cv_acorQ = new NVector();          /* acorQ = yQ_n(m) - yQ_n(0)                    */\r
974           NVector cv_tempvQ = new NVector();         /* temporary storage vector (~ tempv)           */\r
975 \r
976           /*---------------------------\r
977             Sensitivity Related Vectors \r
978             ---------------------------*/\r
979 \r
980           NVector cv_znS[][] = new NVector[L_MAX][];    /* Nordsieck arrays for sensitivities           */\r
981           NVector cv_ewtS[];          /* error weight vectors for sensitivities       */\r
982           NVector cv_yS[];            /* yS=yS0 (allocated by the user)               */\r
983           NVector cv_acorS[];         /* acorS = yS_n(m) - yS_n(0)                    */\r
984           NVector cv_tempvS[];        /* temporary storage vector (~ tempv)           */\r
985           NVector cv_ftempS[];        /* temporary storage vector (~ ftemp)           */\r
986 \r
987           boolean cv_stgr1alloc;  /* Did we allocate ncfS1, ncfnS1, and nniS1?    */\r
988 \r
989           /*--------------------------------------\r
990             Quadrature Sensitivity Related Vectors \r
991             --------------------------------------*/\r
992 \r
993           NVector cv_znQS[][] = new NVector[L_MAX][];   /* Nordsieck arrays for quadr. sensitivities    */\r
994           NVector cv_ewtQS[];         /* error weight vectors for sensitivities       */\r
995           NVector cv_yQS[];           /* Unlike yS, yQS is not allocated by the user  */\r
996           NVector cv_acorQS[];        /* acorQS = yQS_n(m) - yQS_n(0)                 */\r
997           NVector cv_tempvQS[];       /* temporary storage vector (~ tempv)           */\r
998           NVector cv_ftempQ = new NVector();         /* temporary storage vector (~ ftemp)           */\r
999           \r
1000           /*-----------------\r
1001             Tstop information\r
1002             -----------------*/\r
1003 \r
1004           boolean cv_tstopset;\r
1005           double cv_tstop;\r
1006 \r
1007           /*---------\r
1008             Step Data \r
1009             ---------*/\r
1010 \r
1011           int cv_q;                    /* current order                               */\r
1012           int cv_qprime;               /* order to be used on the next step\r
1013                                         * qprime = q-1, q, or q+1                     */\r
1014           int cv_next_q;               /* order to be used on the next step           */\r
1015           int cv_qwait;                /* number of internal steps to wait before\r
1016                                         * considering a change in q                   */\r
1017           int cv_L;                    /* L = q + 1                                   */\r
1018 \r
1019           double cv_hin;\r
1020           double cv_h;               /* current step size                           */\r
1021           double cv_hprime;          /* step size to be used on the next step       */ \r
1022           double cv_next_h;          /* step size to be used on the next step       */ \r
1023           double cv_eta;             /* eta = hprime / h                            */\r
1024           double cv_hscale;          /* value of h used in zn                       */\r
1025           double cv_tn;              /* current internal value of t                 */\r
1026           double cv_tretlast;        /* last value of t returned                    */\r
1027 \r
1028           double cv_tau[] = new double[L_MAX+1];    /* array of previous q+1 successful step\r
1029                                         * sizes indexed from 1 to q+1                 */\r
1030           double cv_tq[] =new double[NUM_TESTS+1]; /* array of test quantities indexed from\r
1031                                         * 1 to NUM_TESTS(=5)                          */\r
1032           double cv_l[] = new double[L_MAX];        /* coefficients of l(x) (degree q poly)        */\r
1033 \r
1034           double cv_rl1;             /* the scalar 1/l[1]                           */\r
1035           double cv_gamma;           /* gamma = h * rl1                             */\r
1036           double cv_gammap;          /* gamma at the last setup call                */\r
1037           double cv_gamrat;          /* gamma / gammap                              */\r
1038 \r
1039           double cv_crate;           /* est. corrector conv. rate in Nls            */\r
1040           double cv_crateS;          /* est. corrector conv. rate in NlsStgr        */\r
1041           double cv_acnrm;           /* | acor |                                    */\r
1042           double cv_acnrmQ;          /* | acorQ |                                   */\r
1043           double cv_acnrmS;          /* | acorS |                                   */\r
1044           double cv_acnrmQS;         /* | acorQS |                                  */\r
1045           double cv_nlscoef;         /* coeficient in nonlinear convergence test    */\r
1046           int  cv_mnewt;               /* Newton iteration counter                    */\r
1047           int  cv_ncfS1[][];              /* Array of Ns local counters for conv.  \r
1048                                         * failures (used in CVStep for STAGGERED1)    */\r
1049                                       /* int[][1]\r
1050 \r
1051           /*------\r
1052             Limits \r
1053             ------*/\r
1054 \r
1055           int cv_qmax;             /* q <= qmax                                       */\r
1056           long cv_mxstep;      /* maximum number of internal steps for one \r
1057                                       user call                                       */\r
1058           int cv_maxcor;           /* maximum number of corrector iterations for \r
1059                                       the solution of the nonlinear equation          */\r
1060           int cv_maxcorS;\r
1061           int cv_mxhnil;           /* max. number of warning messages issued to the\r
1062                                       user that t + h == t for the next internal step */\r
1063           int cv_maxnef;           /* maximum number of error test failures           */\r
1064           int cv_maxncf;           /* maximum number of nonlinear conv. failures      */\r
1065           \r
1066           double cv_hmin;        /* |h| >= hmin                                     */\r
1067           double cv_hmax_inv;    /* |h| <= 1/hmax_inv                               */\r
1068           double cv_etamax;      /* eta <= etamax                                   */\r
1069 \r
1070           /*----------\r
1071             Counters \r
1072             ----------*/\r
1073 \r
1074           long cv_nst;         /* number of internal steps taken                  */\r
1075 \r
1076           long cv_nfe;         /* number of f calls                               */\r
1077           long cv_nfQe;        /* number of fQ calls                              */\r
1078           long cv_nfSe;        /* number of fS calls                              */\r
1079           long cv_nfeS;        /* number of f calls from sensi DQ                 */\r
1080           long cv_nfQSe;       /* number of fQS calls                             */\r
1081           long cv_nfQeS;       /* number of fQ calls from sensi DQ                */\r
1082 \r
1083 \r
1084           long cv_ncfn[] = new long[1];        /* number of corrector convergence failures        */\r
1085           long cv_ncfnS[] = new long[1];       /* number of total sensi. corr. conv. failures     */\r
1086           long cv_ncfnS1[][];     /* number of sensi. corrector conv. failures       */\r
1087                                   // long[][1]\r
1088 \r
1089           long cv_nni;         /* number of nonlinear iterations performed        */\r
1090           long cv_nniS;        /* number of total sensi. nonlinear iterations     */\r
1091           long cv_nniS1[];      /* number of sensi. nonlinear iterations           */\r
1092 \r
1093           long cv_netf[] = new long[1];        /* number of error test failures                   */\r
1094           long cv_netfQ[] = new long[1];       /* number of quadr. error test failures            */\r
1095           long cv_netfS[] = new long[1];       /* number of sensi. error test failures            */\r
1096           long cv_netfQS[] = new long[1];      /* number of quadr. sensi. error test failures     */\r
1097 \r
1098           long cv_nsetups;     /* number of setup calls                           */\r
1099           long cv_nsetupsS;    /* number of setup calls due to sensitivities      */\r
1100 \r
1101           int cv_nhnil;            /* number of messages issued to the user that\r
1102                                       t + h == t for the next iternal step            */\r
1103 \r
1104           /*-----------------------------\r
1105             Space requirements for CVODES \r
1106             -----------------------------*/\r
1107 \r
1108           long cv_lrw1;        /* no. of realtype words in 1 N_Vector y           */ \r
1109           long cv_liw1;        /* no. of integer words in 1 N_Vector y            */ \r
1110           long cv_lrw1Q;       /* no. of realtype words in 1 N_Vector yQ          */ \r
1111           long cv_liw1Q;       /* no. of integer words in 1 N_Vector yQ           */ \r
1112           long cv_lrw;             /* no. of realtype words in CVODES work vectors    */\r
1113           long cv_liw;             /* no. of integer words in CVODES work vectors     */\r
1114 \r
1115           /*----------------\r
1116             Step size ratios\r
1117             ----------------*/\r
1118 \r
1119           double cv_etaqm1;      /* ratio of new to old h for order q-1             */\r
1120           double cv_etaq;        /* ratio of new to old h for order q               */\r
1121           double cv_etaqp1;      /* ratio of new to old h for order q+1             */\r
1122 \r
1123           /*------------------\r
1124             Linear Solver Data \r
1125             ------------------*/\r
1126 \r
1127           /* Linear Solver functions to be called */\r
1128 \r
1129           //int (*cv_linit)(struct CVodeMemRec *cv_mem);\r
1130           int cv_linit_select;\r
1131 \r
1132           //int (*cv_lsetup)(struct CVodeMemRec *cv_mem, int convfail, \r
1133                            //N_Vector ypred, N_Vector fpred, booleantype *jcurPtr, \r
1134                            //N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3); \r
1135           int cv_lsetup_select;\r
1136 \r
1137           //int (*cv_lsolve)(struct CVodeMemRec *cv_mem, N_Vector b, N_Vector weight,\r
1138                            //N_Vector ycur, N_Vector fcur);\r
1139           int cv_lsolve_select;\r
1140 \r
1141           //int (*cv_lfree)(struct CVodeMemRec *cv_mem);\r
1142           int cv_lfree_select;\r
1143 \r
1144           /* Linear Solver specific memory */\r
1145 \r
1146           //void *cv_lmem; \r
1147           CVDlsMemRec cv_lmem;\r
1148 \r
1149           /* Flag to request a call to the setup routine */\r
1150 \r
1151           boolean cv_forceSetup;\r
1152 \r
1153           /*------------\r
1154             Saved Values\r
1155             ------------*/\r
1156 \r
1157           int cv_qu;                   /* last successful q value used                */\r
1158           long cv_nstlp;           /* step number of last setup call              */\r
1159           double cv_h0u;             /* actual initial stepsize                     */\r
1160           double cv_hu;              /* last successful h value used                */\r
1161           double cv_saved_tq5;       /* saved value of tq[5]                        */\r
1162           boolean cv_jcur[] = new boolean[1];         /* is Jacobian info for linear solver current? */\r
1163           int cv_convfail;             /* flag storing previous solver failure mode   */\r
1164           double cv_tolsf;           /* tolerance scale factor                      */\r
1165           int cv_qmax_alloc;           /* qmax used when allocating mem               */\r
1166           int cv_qmax_allocQ;          /* qmax used when allocating quad. mem         */\r
1167           int cv_qmax_allocS;          /* qmax used when allocating sensi. mem        */\r
1168           int cv_qmax_allocQS;         /* qmax used when allocating quad. sensi. mem  */\r
1169           int cv_indx_acor;            /* index of zn vector in which acor is saved   */\r
1170 \r
1171           /*--------------------------------------------------------------------\r
1172             Flags turned ON by CVodeInit, CVodeSensMalloc, and CVodeQuadMalloc \r
1173             and read by CVodeReInit, CVodeSensReInit, and CVodeQuadReInit\r
1174             --------------------------------------------------------------------*/\r
1175 \r
1176           boolean cv_VabstolMallocDone;\r
1177           boolean cv_MallocDone;\r
1178 \r
1179           boolean cv_VabstolQMallocDone;\r
1180           boolean cv_QuadMallocDone;\r
1181 \r
1182           boolean cv_VabstolSMallocDone;\r
1183           boolean cv_SabstolSMallocDone;\r
1184           boolean cv_SensMallocDone;\r
1185 \r
1186           boolean cv_VabstolQSMallocDone;\r
1187           boolean cv_SabstolQSMallocDone;\r
1188           boolean cv_QuadSensMallocDone;\r
1189 \r
1190           /*-------------------------------------------\r
1191             Error handler function and error ouput file \r
1192             -------------------------------------------*/\r
1193 \r
1194           //CVErrHandlerFn cv_ehfun;    /* Error messages are handled by ehfun          */\r
1195           CVodeMemRec cv_eh_data;           /* dats pointer passed to ehfun                 */\r
1196           RandomAccessFile cv_errfp;             /* CVODES error messages are sent to errfp      */    \r
1197 \r
1198           /*-------------------------\r
1199             Stability Limit Detection\r
1200             -------------------------*/\r
1201 \r
1202           boolean cv_sldeton;     /* Is Stability Limit Detection on?             */\r
1203           double cv_ssdat[][] = new double[6][4];    /* scaled data array for STALD                  */\r
1204           int cv_nscon;               /* counter for STALD method                     */\r
1205           long cv_nor;            /* counter for number of order reductions       */\r
1206 \r
1207           /*----------------\r
1208             Rootfinding Data\r
1209             ----------------*/\r
1210 \r
1211           //CVRootFn cv_gfun;        /* Function g for roots sought                     */\r
1212           int cv_gfun;\r
1213           int cv_nrtfn;            /* number of components of g                       */\r
1214           int cv_iroots[];          /* array for root information                      */\r
1215           int cv_rootdir[];         /* array specifying direction of zero-crossing     */\r
1216           double cv_tlo;         /* nearest endpoint of interval in root search     */\r
1217           double cv_thi;         /* farthest endpoint of interval in root search    */\r
1218           double cv_trout;       /* t value returned by rootfinding routine         */\r
1219           double cv_glo[];        /* saved array of g values at t = tlo              */\r
1220           double cv_ghi[];        /* saved array of g values at t = thi              */\r
1221           double cv_grout[];      /* array of g values at t = trout                  */\r
1222           double cv_toutc;       /* copy of tout (if NORMAL mode)                   */\r
1223           double cv_ttol;        /* tolerance on root location trout                */\r
1224           int cv_taskc;            /* copy of parameter itask                         */\r
1225           int cv_irfnd;            /* flag showing whether last step had a root       */\r
1226           long cv_nge;         /* counter for g evaluations                       */\r
1227           boolean cv_gactive[]; /* array with active/inactive event functions      */\r
1228           int cv_mxgnull;          /* number of warning messages about possible g==0  */\r
1229 \r
1230           /*------------------------\r
1231             Adjoint sensitivity data\r
1232             ------------------------*/\r
1233 \r
1234           boolean cv_adj;             /* SUNTRUE if performing ASA                */\r
1235 \r
1236           CVodeMemRec cv_adj_mem; /* Pointer to adjoint memory structure      */\r
1237 \r
1238           boolean cv_adjMallocDone;\r
1239 \r
1240         } ;\r
1241 \r
1242 \r
1243         \r
1244         // Call CVodeCreate to create the solver memory and specify the\r
1245         // Backward Differentiation Formula and the use of a Newton \r
1246         // iteration\r
1247         \r
1248          // CVodeCreate creates an internal memory block for a problem to\r
1249          // be solved by CVODES.\r
1250          //\r
1251          // lmm  is the type of linear multistep method to be used.\r
1252          //      The legal values are CV_ADAMS and CV_BDF (see previous\r
1253          //      description).\r
1254          //\r
1255          // iter  is the type of iteration used to solve the nonlinear\r
1256          //       system that arises during each internal time step.\r
1257          //       The legal values are CV_FUNCTIONAL and CV_NEWTON.\r
1258          //\r
1259          //If successful, CVodeCreate returns a pointer to initialized\r
1260          // problem memory. This pointer should be passed to CVodeInit.\r
1261          // If an initialization error occurs, CVodeCreate prints an error\r
1262          // message to standard err and returns NULL.\r
1263 \r
1264      private CVodeMemRec CVodeCreate(int lmm, int iter) {\r
1265          int maxord;\r
1266          CVodeMemRec cv_mem;\r
1267 \r
1268           /* Test inputs */\r
1269 \r
1270           if ((lmm != CV_ADAMS) && (lmm != CV_BDF)) {\r
1271             cvProcessError(null, 0, "CVODES", "CVodeCreate", MSGCV_BAD_LMM);\r
1272             return null;\r
1273           }\r
1274           \r
1275           if ((iter != CV_FUNCTIONAL) && (iter != CV_NEWTON)) {\r
1276             cvProcessError(null, 0, "CVODES", "CVodeCreate", MSGCV_BAD_ITER);\r
1277             return null;\r
1278           }\r
1279 \r
1280           cv_mem = null;\r
1281           cv_mem = new CVodeMemRec();\r
1282 \r
1283           /* Zero out cv_mem */\r
1284           //memset(cv_mem, 0, sizeof(struct CVodeMemRec));\r
1285 \r
1286           maxord = (lmm == CV_ADAMS) ? ADAMS_Q_MAX : BDF_Q_MAX;\r
1287 \r
1288           /* copy input parameters into cv_mem */\r
1289 \r
1290           cv_mem.cv_lmm  = lmm;\r
1291           cv_mem.cv_iter = iter;\r
1292 \r
1293           /* Set uround */\r
1294 \r
1295           cv_mem.cv_uround = UNIT_ROUNDOFF;\r
1296 \r
1297           /* Set default values for integrator optional inputs */\r
1298 \r
1299           //cv_mem.cv_f          = null;\r
1300           cv_mem.cv_f = -1;\r
1301           //cv_mem.cv_user_data  = null;\r
1302           cv_mem.cv_itol       = CV_NN;\r
1303           cv_mem.cv_user_efun  = false;\r
1304           //cv_mem.cv_efun       = null;\r
1305           cv_mem.cv_efun_select = -1;\r
1306           cv_mem.cv_e_data     = null;\r
1307           //cv_mem.cv_ehfun      = cvErrHandler;\r
1308           cv_mem.cv_eh_data    = cv_mem;\r
1309           //cv_mem.cv_errfp      = stderr;\r
1310           cv_mem.cv_qmax       = maxord;\r
1311           cv_mem.cv_mxstep     = MXSTEP_DEFAULT;\r
1312           cv_mem.cv_mxhnil     = MXHNIL_DEFAULT;\r
1313           cv_mem.cv_sldeton    = false;\r
1314           cv_mem.cv_hin        = ZERO;\r
1315           cv_mem.cv_hmin       = HMIN_DEFAULT;\r
1316           cv_mem.cv_hmax_inv   = HMAX_INV_DEFAULT;\r
1317           cv_mem.cv_tstopset   = false;\r
1318           cv_mem.cv_maxcor     = NLS_MAXCOR;\r
1319           cv_mem.cv_maxnef     = MXNEF;\r
1320           cv_mem.cv_maxncf     = MXNCF;\r
1321           cv_mem.cv_nlscoef    = CORTES;\r
1322 \r
1323           /* Initialize root finding variables */\r
1324 \r
1325           cv_mem.cv_glo        = null;\r
1326           cv_mem.cv_ghi        = null;\r
1327           cv_mem.cv_grout      = null;\r
1328           cv_mem.cv_iroots     = null;\r
1329           cv_mem.cv_rootdir    = null;\r
1330           //cv_mem.cv_gfun       = null;\r
1331           cv_mem.cv_gfun = -1;\r
1332           cv_mem.cv_nrtfn      = 0;  \r
1333           cv_mem.cv_gactive    = null;\r
1334           cv_mem.cv_mxgnull    = 1;\r
1335 \r
1336           /* Set default values for quad. optional inputs */\r
1337 \r
1338           cv_mem.cv_quadr      = false;\r
1339           //cv_mem.cv_fQ         = null;\r
1340           cv_mem.cv_errconQ    = false;\r
1341           cv_mem.cv_itolQ      = CV_NN;\r
1342 \r
1343           /* Set default values for sensi. optional inputs */\r
1344 \r
1345           cv_mem.cv_sensi      = false;\r
1346           //cv_mem.cv_fS_data    = mull;\r
1347           //cv_mem.cv_fS         = cvSensRhsInternalDQ;\r
1348           //cv_mem.cv_fS1        = cvSensRhs1InternalDQ;\r
1349           cv_mem.cv_fSDQ       = true;\r
1350           cv_mem.cv_ifS        = CV_ONESENS;\r
1351           cv_mem.cv_DQtype     = CV_CENTERED;\r
1352           cv_mem.cv_DQrhomax   = ZERO;\r
1353           cv_mem.cv_p          = null;\r
1354           cv_mem.cv_pbar       = null;\r
1355           cv_mem.cv_plist      = null;\r
1356           cv_mem.cv_errconS    = false;\r
1357           cv_mem.cv_maxcorS    = NLS_MAXCOR;\r
1358           cv_mem.cv_ncfS1      = null;\r
1359           cv_mem.cv_ncfnS1     = null;\r
1360           cv_mem.cv_nniS1      = null;\r
1361           cv_mem.cv_itolS      = CV_NN;\r
1362 \r
1363           /* Set default values for quad. sensi. optional inputs */\r
1364 \r
1365           cv_mem.cv_quadr_sensi = false;\r
1366           //cv_mem.cv_fQS         = null;\r
1367           //cv_mem.cv_fQS_data    = null;\r
1368           cv_mem.cv_fQSDQ       = true;\r
1369           cv_mem.cv_errconQS    = false;\r
1370           cv_mem.cv_itolQS      = CV_NN;\r
1371 \r
1372           /* Set default for ASA */\r
1373 \r
1374           cv_mem.cv_adj         = false;\r
1375           cv_mem.cv_adj_mem     = null;\r
1376 \r
1377           /* Set the saved values for qmax_alloc */\r
1378 \r
1379           cv_mem.cv_qmax_alloc  = maxord;\r
1380           cv_mem.cv_qmax_allocQ = maxord;\r
1381           cv_mem.cv_qmax_allocS = maxord;\r
1382 \r
1383           /* Initialize lrw and liw */\r
1384 \r
1385           cv_mem.cv_lrw = 65 + 2*L_MAX + NUM_TESTS;\r
1386           cv_mem.cv_liw = 52;\r
1387 \r
1388           /* No mallocs have been done yet */\r
1389 \r
1390           cv_mem.cv_VabstolMallocDone   = false;\r
1391           cv_mem.cv_MallocDone          = false;\r
1392 \r
1393           cv_mem.cv_VabstolQMallocDone  = false;\r
1394           cv_mem.cv_QuadMallocDone      = false;\r
1395 \r
1396           cv_mem.cv_VabstolSMallocDone  = false;\r
1397           cv_mem.cv_SabstolSMallocDone  = false;\r
1398           cv_mem.cv_SensMallocDone      = false;\r
1399 \r
1400           cv_mem.cv_VabstolQSMallocDone = false;\r
1401           cv_mem.cv_SabstolQSMallocDone = false;\r
1402           cv_mem.cv_QuadSensMallocDone  = false;\r
1403 \r
1404           cv_mem.cv_adjMallocDone       = false;\r
1405 \r
1406           /* Return pointer to CVODES memory block */\r
1407 \r
1408           return cv_mem;\r
1409 \r
1410      }\r
1411      \r
1412      /*\r
1413       * cvProcessError is a high level error handling function.\r
1414       * - If cv_mem==NULL it prints the error message to stderr.\r
1415       * - Otherwise, it sets up and calls the error handling function \r
1416       *   pointed to by cv_ehfun.\r
1417       */\r
1418 \r
1419      void cvProcessError(CVodeMemRec cv_mem, \r
1420                          int error_code,  String module, String fname, \r
1421                          String msgfmt, double... numbers)\r
1422      {\r
1423          MipavUtil.displayError("In module " + module + " in function " + fname + " msgfmt");\r
1424        //va_list ap;\r
1425        //char msg[256];\r
1426 \r
1427        /* Initialize the argument pointer variable \r
1428           (msgfmt is the last required argument to cvProcessError) */\r
1429 \r
1430        //va_start(ap, msgfmt);\r
1431 \r
1432        /* Compose the message */\r
1433 \r
1434        //vsprintf(msg, msgfmt, ap);\r
1435 \r
1436        //if (cv_mem == NULL) {    /* We write to stderr */\r
1437      //#ifndef NO_FPRINTF_OUTPUT\r
1438          //fprintf(stderr, "\n[%s ERROR]  %s\n  ", module, fname);\r
1439          //fprintf(stderr, "%s\n\n", msg);\r
1440      //#endif\r
1441 \r
1442        //} else {                 /* We can call ehfun */\r
1443          //cv_mem->cv_ehfun(error_code, module, fname, msg, cv_mem->cv_eh_data);\r
1444        //}\r
1445 \r
1446        /* Finalize argument processing */\r
1447        //va_end(ap);\r
1448 \r
1449        return;\r
1450      }\r
1451      \r
1452      /*-----------------------------------------------------------------*/\r
1453 \r
1454      /*\r
1455       * CVodeInit\r
1456       * \r
1457       * CVodeInit allocates and initializes memory for a problem. All \r
1458       * problem inputs are checked for errors. If any error occurs during \r
1459       * initialization, it is reported to the file whose file pointer is \r
1460       * errfp and an error flag is returned. Otherwise, it returns CV_SUCCESS\r
1461       */\r
1462 \r
1463      int CVodeInit(CVodeMemRec cv_mem, int f, double t0, NVector y0)\r
1464      {\r
1465        boolean nvectorOK, allocOK;\r
1466        long lrw1[] = new long[1];\r
1467        long liw1[] = new long[1];\r
1468        int i,k;\r
1469 \r
1470        /* Check cvode_mem */\r
1471 \r
1472        if (cv_mem==null) {\r
1473          cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeInit",\r
1474                         MSGCV_NO_MEM);\r
1475          return(CV_MEM_NULL);\r
1476        }\r
1477 \r
1478        /* Check for legal input parameters */\r
1479 \r
1480        if (y0==null) {\r
1481          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
1482                         MSGCV_NULL_Y0);\r
1483          return(CV_ILL_INPUT);\r
1484        }\r
1485 \r
1486        if (f < 0) {\r
1487          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
1488                         MSGCV_NULL_F);\r
1489          return(CV_ILL_INPUT);\r
1490        }\r
1491 \r
1492        /* Test if all required vector operations are implemented */\r
1493 \r
1494        /*nvectorOK = cvCheckNvector(y0);\r
1495        if(!nvectorOK) {\r
1496          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeInit",\r
1497                         MSGCV_BAD_NVECTOR);\r
1498          return(CV_ILL_INPUT);\r
1499        }*/\r
1500 \r
1501        /* Set space requirements for one N_Vector */\r
1502 \r
1503        if (y0.getData() != null) {\r
1504          N_VSpace_Serial(y0, lrw1, liw1);\r
1505        } else {\r
1506          lrw1[0] = 0;\r
1507          liw1[0] = 0;\r
1508        }\r
1509        cv_mem.cv_lrw1 = lrw1[0];\r
1510        cv_mem.cv_liw1 = liw1[0];\r
1511 \r
1512        /* Allocate the vectors (using y0 as a template) */\r
1513 \r
1514       allocOK = cvAllocVectors(cv_mem, y0);\r
1515        if (!allocOK) {\r
1516          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeInit",\r
1517                         MSGCV_MEM_FAIL);\r
1518          return(CV_MEM_FAIL);\r
1519        }\r
1520 \r
1521        /* All error checking is complete at this point */\r
1522 \r
1523        /* Copy the input parameters into CVODES state */\r
1524 \r
1525        //cv_mem.cv_f  = f;\r
1526        cv_mem.cv_tn = t0;\r
1527 \r
1528        /* Set step parameters */\r
1529 \r
1530        cv_mem.cv_q      = 1;\r
1531        cv_mem.cv_L      = 2;\r
1532        cv_mem.cv_qwait  = cv_mem.cv_L;\r
1533        cv_mem.cv_etamax = ETAMX1;\r
1534 \r
1535        cv_mem.cv_qu     = 0;\r
1536        cv_mem.cv_hu     = ZERO;\r
1537        cv_mem.cv_tolsf  = ONE;\r
1538 \r
1539        /* Set the linear solver addresses to NULL.\r
1540           (We check != NULL later, in CVode, if using CV_NEWTON.) */\r
1541 \r
1542        //cv_mem.cv_linit  = null;\r
1543        cv_mem.cv_linit_select = -1;\r
1544        //cv_mem.cv_lsetup = null;\r
1545        cv_mem.cv_lsetup_select = -1;\r
1546        //cv_mem.cv_lsolve = null;\r
1547        cv_mem.cv_lsolve_select = -1;\r
1548        //cv_mem.cv_lfree  = null;\r
1549        cv_mem.cv_lfree_select = -1;\r
1550        //cv_mem.cv_lmem   = null;\r
1551 \r
1552        /* Set forceSetup to SUNFALSE */\r
1553 \r
1554        cv_mem.cv_forceSetup = false;\r
1555 \r
1556        /* Initialize zn[0] in the history array */\r
1557 \r
1558        N_VScale_Serial(ONE, y0, cv_mem.cv_zn[0]);\r
1559 \r
1560        /* Initialize all the counters */\r
1561 \r
1562        cv_mem.cv_nst     = 0;\r
1563        cv_mem.cv_nfe     = 0;\r
1564        cv_mem.cv_ncfn[0]    = 0;\r
1565        cv_mem.cv_netf[0]    = 0;\r
1566        cv_mem.cv_nni     = 0;\r
1567        cv_mem.cv_nsetups = 0;\r
1568        cv_mem.cv_nhnil   = 0;\r
1569        cv_mem.cv_nstlp   = 0;\r
1570        cv_mem.cv_nscon   = 0;\r
1571        cv_mem.cv_nge     = 0;\r
1572 \r
1573        cv_mem.cv_irfnd   = 0;\r
1574 \r
1575        /* Initialize other integrator optional outputs */\r
1576 \r
1577        cv_mem.cv_h0u      = ZERO;\r
1578        cv_mem.cv_next_h   = ZERO;\r
1579        cv_mem.cv_next_q   = 0;\r
1580 \r
1581        /* Initialize Stablilty Limit Detection data */\r
1582        /* NOTE: We do this even if stab lim det was not\r
1583           turned on yet. This way, the user can turn it\r
1584           on at any time */\r
1585 \r
1586        cv_mem.cv_nor = 0;\r
1587        for (i = 1; i <= 5; i++)\r
1588          for (k = 1; k <= 3; k++) \r
1589            cv_mem.cv_ssdat[i-1][k-1] = ZERO;\r
1590 \r
1591        /* Problem has been successfully initialized */\r
1592 \r
1593        cv_mem.cv_MallocDone = true;\r
1594 \r
1595        return(CV_SUCCESS);\r
1596      }\r
1597      \r
1598      void N_VSpace_Serial(NVector v, long lrw[], long liw[])\r
1599      {\r
1600        lrw[0] = v.getData().length;\r
1601        liw[0] = 1;\r
1602 \r
1603        return;\r
1604      }\r
1605 \r
1606      \r
1607      /*\r
1608       * cvCheckNvector\r
1609       * This routine checks if all required vector operations are present.\r
1610       * If any of them is missing it returns SUNFALSE.\r
1611       */\r
1612 \r
1613      /*static booleantype cvCheckNvector(N_Vector tmpl)\r
1614      {\r
1615        if((tmpl->ops->nvclone     == NULL) || // clone\r
1616           (tmpl->ops->nvdestroy   == NULL) || // set null\r
1617           (tmpl->ops->nvlinearsum == NULL) || // (double a, double x[], double b, double y[], double z[])\r
1618                                               // z[i] = a*x[i] + b*y[i];\r
1619           (tmpl->ops->nvconst     == NULL) || // set all to constant\r
1620           (tmpl->ops->nvprod      == NULL) || // (double x[], double y[], double z[])\r
1621                                               // z[i] = x[i] * y[i]\r
1622           (tmpl->ops->nvdiv       == NULL) || // double x[], double y[], double z[])\r
1623                                               // z[i] = x[i]/y[i]\r
1624           (tmpl->ops->nvscale     == NULL) || // (double c, double x[], double z[])\r
1625                                               // z[i] = c * x[i]\r
1626           (tmpl->ops->nvabs       == NULL) || // (double x[], double z[])\r
1627                                               // z[i] = abs(x[i])\r
1628           (tmpl->ops->nvinv       == NULL) || // (double z[], double z[])\r
1629                                               // z[i] = 1.0/x[i]\r
1630           (tmpl->ops->nvaddconst  == NULL) || // (double x[], double b, double z[])\r
1631                                               // z[i] = x[i] + b\r
1632           (tmpl->ops->nvmaxnorm   == NULL) || // (double x[])\r
1633                                               // max(abs(x[i]))\r
1634           (tmpl->ops->nvwrmsnorm  == NULL) || // (double x[], double w[])\r
1635                                               // prod = x[i] * w[i]\r
1636                                               // sum = sum over all i of prod[i]*prod[i]\r
1637                                               // answer = sqrt(sum/N)\r
1638           (tmpl->ops->nvmin       == NULL))   // (double x[])\r
1639                                               // min(x[i])\r
1640          return(SUNFALSE);\r
1641        else\r
1642          return(SUNTRUE);\r
1643      }*/\r
1644      \r
1645      /*\r
1646       * cvAllocVectors\r
1647       *\r
1648       * This routine allocates the CVODES vectors ewt, acor, tempv, ftemp, and\r
1649       * zn[0], ..., zn[maxord].\r
1650       * If all memory allocations are successful, cvAllocVectors returns SUNTRUE. \r
1651       * Otherwise all allocated memory is freed and cvAllocVectors returns SUNFALSE.\r
1652       * This routine also sets the optional outputs lrw and liw, which are\r
1653       * (respectively) the lengths of the real and integer work spaces\r
1654       * allocated here.\r
1655       */\r
1656 \r
1657      private boolean cvAllocVectors(CVodeMemRec cv_mem, NVector tmpl)\r
1658      {\r
1659        int i, j;\r
1660 \r
1661        /* Allocate ewt, acor, tempv, ftemp */\r
1662        \r
1663        cv_mem.cv_ewt = N_VClone(tmpl);\r
1664        if (cv_mem.cv_ewt == null) return(false);\r
1665 \r
1666        cv_mem.cv_acor = N_VClone(tmpl);\r
1667        if (cv_mem.cv_acor == null) {\r
1668          N_VDestroy(cv_mem.cv_ewt);\r
1669          return false;\r
1670        }\r
1671 \r
1672        cv_mem.cv_tempv = N_VClone(tmpl);\r
1673        if (cv_mem.cv_tempv == null) {\r
1674          N_VDestroy(cv_mem.cv_ewt);\r
1675          N_VDestroy(cv_mem.cv_acor);\r
1676          return false;\r
1677        }\r
1678 \r
1679        cv_mem.cv_ftemp = N_VClone(tmpl);\r
1680        if (cv_mem.cv_ftemp == null) {\r
1681          N_VDestroy(cv_mem.cv_tempv);\r
1682          N_VDestroy(cv_mem.cv_ewt);\r
1683          N_VDestroy(cv_mem.cv_acor);\r
1684          return false;\r
1685        }\r
1686 \r
1687        /* Allocate zn[0] ... zn[qmax] */\r
1688 \r
1689        for (j=0; j <= cv_mem.cv_qmax; j++) {\r
1690          cv_mem.cv_zn[j] = N_VClone(tmpl);\r
1691          if (cv_mem.cv_zn[j] == null) {\r
1692            N_VDestroy(cv_mem.cv_ewt);\r
1693            N_VDestroy(cv_mem.cv_acor);\r
1694            N_VDestroy(cv_mem.cv_tempv);\r
1695            N_VDestroy(cv_mem.cv_ftemp);\r
1696            for (i=0; i < j; i++) N_VDestroy(cv_mem.cv_zn[i]);\r
1697            return false;\r
1698          }\r
1699        }\r
1700 \r
1701        /* Update solver workspace lengths  */\r
1702        cv_mem.cv_lrw += (cv_mem.cv_qmax + 5)*cv_mem.cv_lrw1;\r
1703        cv_mem.cv_liw += (cv_mem.cv_qmax + 5)*cv_mem.cv_liw1;\r
1704 \r
1705        /* Store the value of qmax used here */\r
1706        cv_mem.cv_qmax_alloc = cv_mem.cv_qmax;\r
1707 \r
1708        return true;\r
1709      }\r
1710      \r
1711      private NVector N_VClone(NVector y) {\r
1712          NVector x = new NVector();\r
1713          x.data = y.data.clone();\r
1714          x.own_data = y.own_data;\r
1715          return x;\r
1716      }\r
1717 \r
1718      private void N_VDestroy(NVector x) {\r
1719          x.data = null;\r
1720          x = null;\r
1721      }\r
1722      \r
1723      int CVodeSVtolerances(CVodeMemRec cv_mem, double reltol, NVector abstol)\r
1724      {\r
1725 \r
1726        if (cv_mem==null) {\r
1727          cvProcessError(null, CV_MEM_NULL, "CVODES",\r
1728                         "CVodeSVtolerances", MSGCV_NO_MEM);\r
1729          return(CV_MEM_NULL);\r
1730        }\r
1731 \r
1732        if (cv_mem.cv_MallocDone == false) {\r
1733          cvProcessError(cv_mem, CV_NO_MALLOC, "CVODES",\r
1734                         "CVodeSVtolerances", MSGCV_NO_MALLOC);\r
1735          return(CV_NO_MALLOC);\r
1736        }\r
1737 \r
1738        /* Check inputs */\r
1739 \r
1740        if (reltol < ZERO) {\r
1741          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
1742                         "CVodeSVtolerances", MSGCV_BAD_RELTOL);\r
1743          return(CV_ILL_INPUT);\r
1744        }\r
1745 \r
1746        if (N_VMin(abstol) < ZERO) {\r
1747          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES",\r
1748                         "CVodeSVtolerances", MSGCV_BAD_ABSTOL);\r
1749          return(CV_ILL_INPUT);\r
1750        }\r
1751 \r
1752        /* Copy tolerances into memory */\r
1753        \r
1754        if ( !(cv_mem.cv_VabstolMallocDone) ) {\r
1755          cv_mem.cv_Vabstol = N_VClone(cv_mem.cv_ewt);\r
1756          cv_mem.cv_lrw += cv_mem.cv_lrw1;\r
1757          cv_mem.cv_liw += cv_mem.cv_liw1;\r
1758          cv_mem.cv_VabstolMallocDone = true;\r
1759        }\r
1760 \r
1761        cv_mem.cv_reltol = reltol;\r
1762        N_VScale_Serial(ONE, abstol, cv_mem.cv_Vabstol);\r
1763 \r
1764        cv_mem.cv_itol = CV_SV;\r
1765 \r
1766        cv_mem.cv_user_efun = false;\r
1767        //cv_mem.cv_efun = cvEwtSet;\r
1768        cv_mem.cv_efun_select = cvEwtSet_select;\r
1769        cv_mem.cv_e_data = null; /* will be set to cvode_mem in InitialSetup */\r
1770 \r
1771        return(CV_SUCCESS);\r
1772      }\r
1773 \r
1774      private double N_VMin(NVector x) {\r
1775          int i;\r
1776          if ((x == null) || (x.data == null) || (x.data.length == 0)) {\r
1777                  return Double.NaN;\r
1778          }\r
1779          double data[] = x.data;\r
1780          int n = data.length;\r
1781          double min = data[0];\r
1782          for (i = 1; i < n; i++) {\r
1783                  if (data[i] < min) {\r
1784                          min = data[i];\r
1785                  }\r
1786          }\r
1787          return min;\r
1788      }\r
1789      \r
1790      /*\r
1791       * CVodeRootInit\r
1792       *\r
1793       * CVodeRootInit initializes a rootfinding problem to be solved\r
1794       * during the integration of the ODE system.  It loads the root\r
1795       * function pointer and the number of root functions, and allocates\r
1796       * workspace memory.  The return value is CV_SUCCESS = 0 if no errors\r
1797       * occurred, or a negative value otherwise.\r
1798       */\r
1799 \r
1800      int CVodeRootInit(CVodeMemRec cv_mem, int nrtfn, int g)\r
1801      {\r
1802        int i, nrt;\r
1803 \r
1804        /* Check cvode_mem */\r
1805        if (cv_mem==null) {\r
1806          cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeRootInit",\r
1807                         MSGCV_NO_MEM);\r
1808          return(CV_MEM_NULL);\r
1809        }\r
1810 \r
1811        nrt = (nrtfn < 0) ? 0 : nrtfn;\r
1812 \r
1813        /* If rerunning CVodeRootInit() with a different number of root\r
1814           functions (changing number of gfun components), then free\r
1815           currently held memory resources */\r
1816        if ((nrt != cv_mem.cv_nrtfn) && (cv_mem.cv_nrtfn > 0)) {\r
1817          cv_mem.cv_glo = null;\r
1818          cv_mem.cv_ghi = null;\r
1819          cv_mem.cv_grout = null;\r
1820          cv_mem.cv_iroots = null;\r
1821          cv_mem.cv_rootdir = null;\r
1822          cv_mem.cv_gactive = null;\r
1823 \r
1824          cv_mem.cv_lrw -= 3 * (cv_mem.cv_nrtfn);\r
1825          cv_mem.cv_liw -= 3 * (cv_mem.cv_nrtfn);\r
1826        }\r
1827 \r
1828        /* If CVodeRootInit() was called with nrtfn == 0, then set cv_nrtfn to\r
1829           zero and cv_gfun to NULL before returning */\r
1830        if (nrt == 0) {\r
1831          cv_mem.cv_nrtfn = nrt;\r
1832          cv_mem.cv_gfun = -1;\r
1833          return(CV_SUCCESS);\r
1834        }\r
1835 \r
1836        /* If rerunning CVodeRootInit() with the same number of root functions\r
1837           (not changing number of gfun components), then check if the root\r
1838           function argument has changed */\r
1839        /* If g != NULL then return as currently reserved memory resources\r
1840           will suffice */\r
1841        if (nrt == cv_mem.cv_nrtfn) {\r
1842          if (g != cv_mem.cv_gfun) {\r
1843            if (g < 0) {\r
1844              cv_mem.cv_glo = null;\r
1845              cv_mem.cv_ghi = null;\r
1846              cv_mem.cv_grout = null;\r
1847              cv_mem.cv_iroots = null;\r
1848              cv_mem.cv_rootdir = null;\r
1849              cv_mem.cv_gactive = null;\r
1850 \r
1851              cv_mem.cv_lrw -= 3*nrt;\r
1852              cv_mem.cv_liw -= 3*nrt;\r
1853 \r
1854              cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeRootInit",\r
1855                             MSGCV_NULL_G);\r
1856              return(CV_ILL_INPUT);\r
1857            }\r
1858            else {\r
1859              cv_mem.cv_gfun = g;\r
1860              return(CV_SUCCESS);\r
1861            }\r
1862          }\r
1863          else return(CV_SUCCESS);\r
1864        }\r
1865 \r
1866        /* Set variable values in CVode memory block */\r
1867        cv_mem.cv_nrtfn = nrt;\r
1868        if (g < 0) {\r
1869          cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeRootInit",\r
1870                         MSGCV_NULL_G);\r
1871          return(CV_ILL_INPUT);\r
1872        }\r
1873        else cv_mem.cv_gfun = g;\r
1874 \r
1875        /* Allocate necessary memory and return */\r
1876        cv_mem.cv_glo = null;\r
1877        cv_mem.cv_glo = new double[nrt];\r
1878        if (cv_mem.cv_glo == null) {\r
1879          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1880                         MSGCV_MEM_FAIL);\r
1881          return(CV_MEM_FAIL);\r
1882        }\r
1883          \r
1884        cv_mem.cv_ghi = null;\r
1885        cv_mem.cv_ghi = new double[nrt];\r
1886        if (cv_mem.cv_ghi == null) {\r
1887          cv_mem.cv_glo = null;\r
1888          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1889                         MSGCV_MEM_FAIL);\r
1890          return(CV_MEM_FAIL);\r
1891        }\r
1892          \r
1893        cv_mem.cv_grout = null;\r
1894        cv_mem.cv_grout = new double[nrt];\r
1895        if (cv_mem.cv_grout == null) {\r
1896          cv_mem.cv_glo = null;\r
1897          cv_mem.cv_ghi = null;\r
1898          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1899                         MSGCV_MEM_FAIL);\r
1900          return(CV_MEM_FAIL);\r
1901        }\r
1902 \r
1903        cv_mem.cv_iroots = null;\r
1904        cv_mem.cv_iroots = new int[nrt];\r
1905        if (cv_mem.cv_iroots == null) {\r
1906          cv_mem.cv_glo = null;\r
1907          cv_mem.cv_ghi = null;\r
1908          cv_mem.cv_grout = null;\r
1909          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1910                         MSGCV_MEM_FAIL);\r
1911          return(CV_MEM_FAIL);\r
1912        }\r
1913 \r
1914        cv_mem.cv_rootdir = null;\r
1915        cv_mem.cv_rootdir = new int[nrt];\r
1916        if (cv_mem.cv_rootdir == null) {\r
1917          cv_mem.cv_glo = null; \r
1918          cv_mem.cv_ghi = null;\r
1919          cv_mem.cv_grout = null;\r
1920          cv_mem.cv_iroots = null;\r
1921          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1922                         MSGCV_MEM_FAIL);\r
1923          return(CV_MEM_FAIL);\r
1924        }\r
1925 \r
1926 \r
1927        cv_mem.cv_gactive = null;\r
1928        cv_mem.cv_gactive = new boolean[nrt];\r
1929        if (cv_mem.cv_gactive == null) {\r
1930          cv_mem.cv_glo = null; \r
1931          cv_mem.cv_ghi = null;\r
1932          cv_mem.cv_grout = null;\r
1933          cv_mem.cv_iroots = null;\r
1934          cv_mem.cv_rootdir = null;\r
1935          cvProcessError(cv_mem, CV_MEM_FAIL, "CVODES", "CVodeRootInit",\r
1936                         MSGCV_MEM_FAIL);\r
1937          return(CV_MEM_FAIL);\r
1938        }\r
1939 \r
1940 \r
1941        /* Set default values for rootdir (both directions) */\r
1942        for(i=0; i<nrt; i++) cv_mem.cv_rootdir[i] = 0;\r
1943 \r
1944        /* Set default values for gactive (all active) */\r
1945        for(i=0; i<nrt; i++) cv_mem.cv_gactive[i] = true;\r
1946 \r
1947        cv_mem.cv_lrw += 3*nrt;\r
1948        cv_mem.cv_liw += 3*nrt;\r
1949 \r
1950        return(CV_SUCCESS);\r
1951      }\r
1952      \r
1953      class SUNLinearSolver {\r
1954          long N; // Size of the linear system, the number of matrix rows\r
1955          int pivots[]; // Array of size N, index array for partial pivoting in LU factorization\r
1956          int last_flag; // last error return flag from internal setup/solve\r
1957          SUNLinearSolver_Type type;\r
1958      }\r
1959      \r
1960      private SUNLinearSolver SUNDenseLinearSolver(NVector y, double A[][]) {\r
1961          int matrixRows = A.length;\r
1962          int matrixColumns = A[0].length;\r
1963          if (matrixRows != matrixColumns) {\r
1964                  MipavUtil.displayError("Did not have the matrix rows = matrix columns required for a LinearSolver");\r
1965                  return null;\r
1966          }\r
1967          int vecLength = y.getData().length;\r
1968          if (matrixRows != vecLength) {\r
1969                  MipavUtil.displayError("Did not have the matrix rows equal to vector length required for a LinearSolver");\r
1970                  return null;\r
1971          }\r
1972          SUNLinearSolver S = new SUNLinearSolver();\r
1973          S.N = matrixRows;\r
1974          S.last_flag = 0;\r
1975          S.pivots = new int[matrixRows];\r
1976          S.type = SUNLinearSolver_Type.SUNLINEARSOLVER_DIRECT;\r
1977          return S;\r
1978      }\r
1979      \r
1980      /*---------------------------------------------------------------\r
1981      CVDlsSetLinearSolver specifies the direct linear solver.\r
1982     ---------------------------------------------------------------*/\r
1983     int CVDlsSetLinearSolver(CVodeMemRec cv_mem, SUNLinearSolver LS,\r
1984                            double A[][])\r
1985     {\r
1986       CVDlsMemRec cvdls_mem;\r
1987 \r
1988       /* Return immediately if any input is NULL */\r
1989       if (cv_mem == null) {\r
1990         cvProcessError(null, CVDLS_MEM_NULL, "CVSDLS", \r
1991                        "CVDlsSetLinearSolver", MSGD_CVMEM_NULL);\r
1992         return(CVDLS_MEM_NULL);\r
1993       }\r
1994       if ( (LS == null)  || (A == null) ) {\r
1995         cvProcessError(null, CVDLS_ILL_INPUT, "CVSDLS", \r
1996                        "CVDlsSetLinearSolver",\r
1997                         "Both LS and A must be non-NULL");\r
1998         return(CVDLS_ILL_INPUT);\r
1999       }\r
2000 \r
2001       /* Test if solver and vector are compatible with DLS */\r
2002       if (LS.type != SUNLinearSolver_Type.SUNLINEARSOLVER_DIRECT) {\r
2003         cvProcessError(cv_mem, CVDLS_ILL_INPUT, "CVSDLS", \r
2004                        "CVDlsSetLinearSolver", \r
2005                        "Non-direct LS supplied to CVDls interface");\r
2006         return(CVDLS_ILL_INPUT);\r
2007       }\r
2008       //if (cv_mem.cv_tempv.ops.nvgetarraypointer == null ||\r
2009           //cv_mem.cv_tempv.ops.nvsetarraypointer ==  null) {\r
2010         //cvProcessError(cv_mem, CVDLS_ILL_INPUT, "CVSDLS", \r
2011                        //"CVDlsSetLinearSolver", MSGD_BAD_NVECTOR);\r
2012         //return(CVDLS_ILL_INPUT);\r
2013       //}\r
2014 \r
2015       /* free any existing system solver attached to CVode */\r
2016       if (cv_mem.cv_lfree_select > 0)  cv_lfree(cv_mem, cv_mem.cv_lfree_select);\r
2017 \r
2018       /* Set four main system linear solver function fields in cv_mem */\r
2019       //cv_mem.cv_linit  = cvDlsInitialize;\r
2020       cv_mem.cv_linit_select = cvDlsInitialize_select;\r
2021       //cv_mem.cv_lsetup = cvDlsSetup;\r
2022       cv_mem.cv_lsetup_select = cvDlsSetup_select;\r
2023       //cv_mem.cv_lsolve = cvDlsSolve;\r
2024       cv_mem.cv_lsolve_select = cvDlsSolve_select;\r
2025       //cv_mem.cv_lfree  = cvDlsFree;\r
2026       cv_mem.cv_lfree_select = cvDlsFree_select;\r
2027       \r
2028       /* Get memory for CVDlsMemRec */\r
2029       cvdls_mem = null;\r
2030       cvdls_mem = new CVDlsMemRec();\r
2031 \r
2032       /* set SUNLinearSolver pointer */\r
2033       cvdls_mem.LS = LS;\r
2034       \r
2035       /* Initialize Jacobian-related data */\r
2036       cvdls_mem.jacDQ = true;\r
2037       cvdls_mem.jac = cvDlsDQJac;\r
2038       cvdls_mem.J_data = cv_mem;\r
2039       cvdls_mem.last_flag = CVDLS_SUCCESS;\r
2040 \r
2041       /* Initialize counters */\r
2042       cvDlsInitializeCounters(cvdls_mem);\r
2043 \r
2044       /* Store pointer to A and create saved_J */\r
2045       cvdls_mem.A = A;\r
2046       cvdls_mem.savedJ = A.clone();\r
2047       if (cvdls_mem.savedJ == null) {\r
2048         cvProcessError(cv_mem, CVDLS_MEM_FAIL, "CVSDLS", \r
2049                         "CVDlsSetLinearSolver", MSGD_MEM_FAIL);\r
2050         cvdls_mem = null;\r
2051         return(CVDLS_MEM_FAIL);\r
2052       }\r
2053 \r
2054       /* Allocate memory for x */\r
2055       cvdls_mem.x = N_VClone(cv_mem.cv_tempv);\r
2056       if (cvdls_mem.x == null) {\r
2057         cvProcessError(cv_mem, CVDLS_MEM_FAIL, "CVSDLS", \r
2058                         "CVDlsSetLinearSolver", MSGD_MEM_FAIL);\r
2059         cvdls_mem.savedJ = null;\r
2060         cvdls_mem = null;\r
2061         return(CVDLS_MEM_FAIL);\r
2062       }\r
2063       /* Attach linear solver memory to integrator memory */\r
2064       cv_mem.cv_lmem = cvdls_mem;\r
2065 \r
2066       return(CVDLS_SUCCESS);\r
2067     }\r
2068     \r
2069     //-----------------------------------------------------------------\r
2070     //cvDlsInitializeCounters\r
2071     //-----------------------------------------------------------------\r
2072     //This routine resets the counters inside the CVDlsMem object.\r
2073     //-----------------------------------------------------------------*/\r
2074   int cvDlsInitializeCounters(CVDlsMemRec cvdls_mem)\r
2075   {\r
2076     cvdls_mem.nje   = 0;\r
2077     cvdls_mem.nfeDQ = 0;\r
2078     cvdls_mem.nstlj = 0;\r
2079     return(0);\r
2080   }\r
2081 \r
2082     \r
2083     \r
2084     \r
2085     //-----------------------------------------------------------------\r
2086     //CVDlsMem is pointer to a CVDlsMemRec structure.\r
2087     //-----------------------------------------------------------------*/\r
2088 \r
2089   class CVDlsMemRec {\r
2090 \r
2091     boolean jacDQ;    /* true if using internal DQ Jacobian approx. */\r
2092     //CVDlsJacFn jac;       /* Jacobian routine to be called                 */\r
2093     int jac;\r
2094     //void *J_data;         /* data pointer passed to jac                    */\r
2095     CVodeMemRec J_data;\r
2096 \r
2097     double A[][];          /* A = I - gamma * df/dy                         */\r
2098     double savedJ[][];     /* savedJ = old Jacobian                         */\r
2099 \r
2100     SUNLinearSolver LS;   /* generic direct linear solver object           */\r
2101 \r
2102     NVector x;           /* solution vector used by SUNLinearSolver       */\r
2103     \r
2104     long nstlj;       /* nstlj = nst at last Jacobian eval.            */\r
2105 \r
2106     long nje;         /* nje = no. of calls to jac                     */\r
2107 \r
2108     long nfeDQ;       /* no. of calls to f due to DQ Jacobian approx.  */\r
2109 \r
2110     int last_flag;   /* last error return flag                        */\r
2111 \r
2112   };\r
2113 \r
2114 \r
2115   /* CVDlsSetJacFn specifies the Jacobian function. */\r
2116   int CVDlsSetJacFn(CVodeMemRec cv_mem, int jac)\r
2117   {\r
2118     CVDlsMemRec cvdls_mem;\r
2119 \r
2120     /* Return immediately if cvode_mem or cv_mem->cv_lmem are NULL */\r
2121     if (cv_mem == null) {\r
2122       cvProcessError(null, CVDLS_MEM_NULL, "CVSDLS",\r
2123                      "CVDlsSetJacFn", MSGD_CVMEM_NULL);\r
2124       return(CVDLS_MEM_NULL);\r
2125     }\r
2126   \r
2127     if (cv_mem.cv_lmem == null) {\r
2128       cvProcessError(cv_mem, CVDLS_LMEM_NULL, "CVSDLS",\r
2129                      "CVDlsSetJacFn", MSGD_LMEM_NULL);\r
2130       return(CVDLS_LMEM_NULL);\r
2131     }\r
2132     cvdls_mem = cv_mem.cv_lmem;\r
2133 \r
2134     if (jac >= 0) {\r
2135       cvdls_mem.jacDQ  = false;\r
2136       cvdls_mem.jac    = jac;\r
2137       cvdls_mem.J_data = cv_mem.cv_user_data;\r
2138     } else {\r
2139       cvdls_mem.jacDQ  = true;\r
2140       cvdls_mem.jac    = cvDlsDQJac;\r
2141       cvdls_mem.J_data = cv_mem;\r
2142     }\r
2143 \r
2144     return(CVDLS_SUCCESS);\r
2145   }\r
2146   \r
2147   /*\r
2148    * CVode\r
2149    *\r
2150    * This routine is the main driver of the CVODES package. \r
2151    *\r
2152    * It integrates over a time interval defined by the user, by calling\r
2153    * cvStep to do internal time steps.\r
2154    *\r
2155    * The first time that CVode is called for a successfully initialized\r
2156    * problem, it computes a tentative initial step size h.\r
2157    *\r
2158    * CVode supports two modes, specified by itask: CV_NORMAL, CV_ONE_STEP.\r
2159    * In the CV_NORMAL mode, the solver steps until it reaches or passes tout\r
2160    * and then interpolates to obtain y(tout).\r
2161    * In the CV_ONE_STEP mode, it takes one internal step and returns.\r
2162    */\r
2163 \r
2164   int CVode(CVodeMemRec cv_mem, double tout, NVector yout, \r
2165             double tret[], int itask)\r
2166   {\r
2167     long nstloc; \r
2168     int retval, hflag, kflag, istate, is, ir, ier, irfndp;\r
2169     double troundoff, tout_hin, rh, nrm;\r
2170     boolean inactive_roots;\r
2171 \r
2172     /*\r
2173      * -------------------------------------\r
2174      * 1. Check and process inputs\r
2175      * -------------------------------------\r
2176      */\r
2177 \r
2178     /* Check if cvode_mem exists */\r
2179     if (cv_mem == null) {\r
2180       cvProcessError(null, CV_MEM_NULL, "CVODES", "CVode",\r
2181                      MSGCV_NO_MEM);\r
2182       return(CV_MEM_NULL);\r
2183     }\r
2184 \r
2185     /* Check if cvode_mem was allocated */\r
2186     if (cv_mem.cv_MallocDone == false) {\r
2187       cvProcessError(cv_mem, CV_NO_MALLOC, "CVODES", "CVode",\r
2188                      MSGCV_NO_MALLOC);\r
2189       return(CV_NO_MALLOC);\r
2190     }\r
2191     \r
2192     /* Check for yout != NULL */\r
2193     if ((cv_mem.cv_y = yout) == null) {\r
2194       cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2195                      MSGCV_YOUT_NULL);\r
2196       return(CV_ILL_INPUT);\r
2197     }\r
2198     \r
2199     /* Check for tret != NULL */\r
2200     if (tret == null) {\r
2201       cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2202                      MSGCV_TRET_NULL);\r
2203       return(CV_ILL_INPUT);\r
2204     }\r
2205 \r
2206     /* Check for valid itask */\r
2207     if ( (itask != CV_NORMAL) && (itask != CV_ONE_STEP) ) {\r
2208       cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2209                      MSGCV_BAD_ITASK);\r
2210       return(CV_ILL_INPUT);\r
2211     }\r
2212 \r
2213     if (itask == CV_NORMAL) cv_mem.cv_toutc = tout;\r
2214     cv_mem.cv_taskc = itask;\r
2215 \r
2216     /*\r
2217      * ----------------------------------------\r
2218      * 2. Initializations performed only at\r
2219      *    the first step (nst=0):\r
2220      *    - initial setup\r
2221      *    - initialize Nordsieck history array\r
2222      *    - compute initial step size\r
2223      *    - check for approach to tstop\r
2224      *    - check for approach to a root\r
2225      * ----------------------------------------\r
2226      */\r
2227 \r
2228     if (cv_mem.cv_nst == 0) {\r
2229 \r
2230       cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2231 \r
2232       /* Check inputs for corectness */\r
2233 \r
2234       ier = cvInitialSetup(cv_mem);\r
2235       if (ier!= CV_SUCCESS) return(ier);\r
2236 \r
2237       /* \r
2238        * Call f at (t0,y0), set zn[1] = y'(t0). \r
2239        * If computing any quadratures, call fQ at (t0,y0), set znQ[1] = yQ'(t0)\r
2240        * If computing sensitivities, call fS at (t0,y0,yS0), set znS[1][is] = yS'(t0), is=1,...,Ns.\r
2241        * If computing quadr. sensi., call fQS at (t0,y0,yS0), set znQS[1][is] = yQS'(t0), is=1,...,Ns.\r
2242        */\r
2243 \r
2244       retval = f(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
2245                             cv_mem.cv_zn[1], cv_mem.cv_user_data); \r
2246       cv_mem.cv_nfe++;\r
2247       if (retval < 0) {\r
2248         cvProcessError(cv_mem, CV_RHSFUNC_FAIL, "CVODES", "CVode",\r
2249                        MSGCV_RHSFUNC_FAILED, cv_mem.cv_tn);\r
2250         return(CV_RHSFUNC_FAIL);\r
2251       }\r
2252       if (retval > 0) {\r
2253         cvProcessError(cv_mem, CV_FIRST_RHSFUNC_ERR, "CVODES", "CVode",\r
2254                        MSGCV_RHSFUNC_FIRST);\r
2255         return(CV_FIRST_RHSFUNC_ERR);\r
2256       }\r
2257 \r
2258       if (cv_mem.cv_quadr) {\r
2259         retval = fQ(cv_mem.cv_tn, cv_mem.cv_zn[0],\r
2260                                cv_mem.cv_znQ[1], cv_mem.cv_user_data);\r
2261         cv_mem.cv_nfQe++;\r
2262         if (retval < 0) {\r
2263           cvProcessError(cv_mem, CV_QRHSFUNC_FAIL, "CVODES", "CVode",\r
2264                          MSGCV_QRHSFUNC_FAILED, cv_mem.cv_tn);\r
2265           return(CV_QRHSFUNC_FAIL);\r
2266         }\r
2267         if (retval > 0) {\r
2268           cvProcessError(cv_mem, CV_FIRST_QRHSFUNC_ERR, "CVODES",\r
2269                          "CVode", MSGCV_QRHSFUNC_FIRST);\r
2270           return(CV_FIRST_QRHSFUNC_ERR);\r
2271         }\r
2272       }\r
2273 \r
2274       if (cv_mem.cv_sensi) {\r
2275         retval = cvSensRhsWrapper(cv_mem, cv_mem.cv_tn, cv_mem.cv_zn[0],\r
2276                                   cv_mem.cv_zn[1], cv_mem.cv_znS[0],\r
2277                                   cv_mem.cv_znS[1], cv_mem.cv_tempv,\r
2278                                   cv_mem.cv_ftemp);\r
2279         if (retval < 0) {\r
2280           cvProcessError(cv_mem, CV_SRHSFUNC_FAIL, "CVODES", "CVode",\r
2281                          MSGCV_SRHSFUNC_FAILED, cv_mem.cv_tn);\r
2282           return(CV_SRHSFUNC_FAIL);\r
2283         } \r
2284         if (retval > 0) {\r
2285           cvProcessError(cv_mem, CV_FIRST_SRHSFUNC_ERR, "CVODES",\r
2286                          "CVode", MSGCV_SRHSFUNC_FIRST);\r
2287           return(CV_FIRST_SRHSFUNC_ERR);\r
2288         }\r
2289       }\r
2290 \r
2291       if (cv_mem.cv_quadr_sensi) {\r
2292         retval = cvQuadSensRhsInternalDQ(cv_mem.cv_Ns, cv_mem.cv_tn, cv_mem.cv_zn[0],\r
2293                                 cv_mem.cv_znS[0], cv_mem.cv_znQ[1],\r
2294                                 cv_mem.cv_znQS[1], cv_mem.cv_fQS_data,\r
2295                                 cv_mem.cv_tempv, cv_mem.cv_tempvQ); \r
2296         cv_mem.cv_nfQSe++;\r
2297         if (retval < 0) {\r
2298           cvProcessError(cv_mem, CV_QSRHSFUNC_FAIL, "CVODES", "CVode",\r
2299                          MSGCV_QSRHSFUNC_FAILED, cv_mem.cv_tn);\r
2300           return(CV_QSRHSFUNC_FAIL);\r
2301         } \r
2302         if (retval > 0) {\r
2303           cvProcessError(cv_mem, CV_FIRST_QSRHSFUNC_ERR, "CVODES",\r
2304                          "CVode", MSGCV_QSRHSFUNC_FIRST);\r
2305           return(CV_FIRST_QSRHSFUNC_ERR);\r
2306         }\r
2307       } \r
2308 \r
2309       /* Test input tstop for legality. */\r
2310 \r
2311       if (cv_mem.cv_tstopset) {\r
2312         if ( (cv_mem.cv_tstop - cv_mem.cv_tn)*(tout - cv_mem.cv_tn) <= ZERO ) {\r
2313           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2314                          MSGCV_BAD_TSTOP, cv_mem.cv_tstop, cv_mem.cv_tn);\r
2315           return(CV_ILL_INPUT);\r
2316         }\r
2317       }\r
2318 \r
2319       /* Set initial h (from H0 or cvHin). */\r
2320       \r
2321       cv_mem.cv_h = cv_mem.cv_hin;\r
2322       if ( (cv_mem.cv_h != ZERO) && ((tout-cv_mem.cv_tn)*cv_mem.cv_h < ZERO) ) {\r
2323         cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode", MSGCV_BAD_H0);\r
2324         return(CV_ILL_INPUT);\r
2325       }\r
2326       if (cv_mem.cv_h == ZERO) {\r
2327         tout_hin = tout;\r
2328         if ( cv_mem.cv_tstopset &&\r
2329              (tout-cv_mem.cv_tn)*(tout-cv_mem.cv_tstop) > ZERO )\r
2330           tout_hin = cv_mem.cv_tstop; \r
2331         hflag = cvHin(cv_mem, tout_hin);\r
2332         if (hflag != CV_SUCCESS) {\r
2333           istate = cvHandleFailure(cv_mem, hflag);\r
2334           return(istate);\r
2335         }\r
2336       }\r
2337       rh = Math.abs(cv_mem.cv_h)*cv_mem.cv_hmax_inv;\r
2338       if (rh > ONE) cv_mem.cv_h /= rh;\r
2339       if (Math.abs(cv_mem.cv_h) < cv_mem.cv_hmin)\r
2340         cv_mem.cv_h *= cv_mem.cv_hmin/Math.abs(cv_mem.cv_h);\r
2341 \r
2342       /* Check for approach to tstop */\r
2343 \r
2344       if (cv_mem.cv_tstopset) {\r
2345         if ( (cv_mem.cv_tn + cv_mem.cv_h - cv_mem.cv_tstop)*cv_mem.cv_h > ZERO ) \r
2346           cv_mem.cv_h = (cv_mem.cv_tstop - cv_mem.cv_tn)*(ONE-FOUR*cv_mem.cv_uround);\r
2347       }\r
2348 \r
2349       /* \r
2350        * Scale zn[1] by h.\r
2351        * If computing any quadratures, scale znQ[1] by h.\r
2352        * If computing sensitivities,  scale znS[1][is] by h. \r
2353        * If computing quadrature sensitivities,  scale znQS[1][is] by h. \r
2354        */\r
2355 \r
2356       cv_mem.cv_hscale = cv_mem.cv_h;\r
2357       cv_mem.cv_h0u    = cv_mem.cv_h;\r
2358       cv_mem.cv_hprime = cv_mem.cv_h;\r
2359 \r
2360       N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_zn[1], cv_mem.cv_zn[1]);\r
2361       \r
2362       if (cv_mem.cv_quadr)\r
2363         N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_znQ[1], cv_mem.cv_znQ[1]);\r
2364 \r
2365       if (cv_mem.cv_sensi)\r
2366         for (is=0; is<cv_mem.cv_Ns; is++) \r
2367           N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_znS[1][is], cv_mem.cv_znS[1][is]);\r
2368 \r
2369       if (cv_mem.cv_quadr_sensi)\r
2370         for (is=0; is<cv_mem.cv_Ns; is++) \r
2371           N_VScale_Serial(cv_mem.cv_h, cv_mem.cv_znQS[1][is], cv_mem.cv_znQS[1][is]);\r
2372       \r
2373       /* Check for zeros of root function g at and near t0. */\r
2374 \r
2375       if (cv_mem.cv_nrtfn > 0) {\r
2376 \r
2377         retval = cvRcheck1(cv_mem);\r
2378 \r
2379         if (retval == CV_RTFUNC_FAIL) {\r
2380           cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck1",\r
2381                          MSGCV_RTFUNC_FAILED, cv_mem.cv_tn);\r
2382           return(CV_RTFUNC_FAIL);\r
2383         }\r
2384 \r
2385       }\r
2386 \r
2387     } /* end first call block */ // if (cv_mem.cv_nst == 0)\r
2388 \r
2389     /*\r
2390      * ------------------------------------------------------\r
2391      * 3. At following steps, perform stop tests:\r
2392      *    - check for root in last step\r
2393      *    - check if we passed tstop\r
2394      *    - check if we passed tout (NORMAL mode)\r
2395      *    - check if current tn was returned (ONE_STEP mode)\r
2396      *    - check if we are close to tstop\r
2397      *      (adjust step size if needed)\r
2398      * -------------------------------------------------------\r
2399      */\r
2400 \r
2401     if (cv_mem.cv_nst > 0) {\r
2402 \r
2403       /* Estimate an infinitesimal time interval to be used as\r
2404          a roundoff for time quantities (based on current time \r
2405          and step size) */\r
2406       troundoff = FUZZ_FACTOR * cv_mem.cv_uround *\r
2407         (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_h));\r
2408 \r
2409       /* First check for a root in the last step taken, other than the\r
2410          last root found, if any.  If itask = CV_ONE_STEP and y(tn) was not\r
2411          returned because of an intervening root, return y(tn) now.     */\r
2412       if (cv_mem.cv_nrtfn > 0) {\r
2413         \r
2414         irfndp = cv_mem.cv_irfnd;\r
2415         \r
2416         retval = cvRcheck2(cv_mem);\r
2417 \r
2418         if (retval == CLOSERT) {\r
2419           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvRcheck2",\r
2420                          MSGCV_CLOSE_ROOTS, cv_mem.cv_tlo);\r
2421           return(CV_ILL_INPUT);\r
2422         } else if (retval == CV_RTFUNC_FAIL) {\r
2423           cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck2",\r
2424                          MSGCV_RTFUNC_FAILED, cv_mem.cv_tlo);\r
2425           return(CV_RTFUNC_FAIL);\r
2426         } else if (retval == RTFOUND) {\r
2427           cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tlo;\r
2428           return(CV_ROOT_RETURN);\r
2429         }\r
2430         \r
2431         /* If tn is distinct from tretlast (within roundoff),\r
2432            check remaining interval for roots */\r
2433         if ( Math.abs(cv_mem.cv_tn - cv_mem.cv_tretlast) > troundoff ) {\r
2434 \r
2435           retval = cvRcheck3(cv_mem);\r
2436 \r
2437           if (retval == CV_SUCCESS) {     /* no root found */\r
2438             cv_mem.cv_irfnd = 0;\r
2439             if ((irfndp == 1) && (itask == CV_ONE_STEP)) {\r
2440               cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2441               N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
2442               return(CV_SUCCESS);\r
2443             }\r
2444           } else if (retval == RTFOUND) {  /* a new root was found */\r
2445             cv_mem.cv_irfnd = 1;\r
2446             cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tlo;\r
2447             return(CV_ROOT_RETURN);\r
2448           } else if (retval == CV_RTFUNC_FAIL) {  /* g failed */\r
2449             cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck3", \r
2450                            MSGCV_RTFUNC_FAILED, cv_mem.cv_tlo);\r
2451             return(CV_RTFUNC_FAIL);\r
2452           }\r
2453 \r
2454         }\r
2455         \r
2456       } /* end of root stop check */\r
2457       \r
2458       /* In CV_NORMAL mode, test if tout was reached */\r
2459       if ( (itask == CV_NORMAL) && ((cv_mem.cv_tn-tout)*cv_mem.cv_h >= ZERO) ) {\r
2460         cv_mem.cv_tretlast = tret[0] = tout;\r
2461         ier =  CVodeGetDky(cv_mem, tout, 0, yout);\r
2462         if (ier != CV_SUCCESS) {\r
2463           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2464                          MSGCV_BAD_TOUT, tout);\r
2465           return(CV_ILL_INPUT);\r
2466         }\r
2467         return(CV_SUCCESS);\r
2468       }\r
2469       \r
2470       /* In CV_ONE_STEP mode, test if tn was returned */\r
2471       if ( itask == CV_ONE_STEP &&\r
2472            Math.abs(cv_mem.cv_tn - cv_mem.cv_tretlast) > troundoff ) {\r
2473         cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2474         N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
2475         return(CV_SUCCESS);\r
2476       }\r
2477       \r
2478       /* Test for tn at tstop or near tstop */\r
2479       if ( cv_mem.cv_tstopset ) {\r
2480         \r
2481         if ( Math.abs(cv_mem.cv_tn - cv_mem.cv_tstop) <= troundoff ) {\r
2482           ier =  CVodeGetDky(cv_mem, cv_mem.cv_tstop, 0, yout);\r
2483           if (ier != CV_SUCCESS) {\r
2484             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2485                            MSGCV_BAD_TSTOP, cv_mem.cv_tstop, cv_mem.cv_tn);\r
2486             return(CV_ILL_INPUT);\r
2487           }\r
2488           cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tstop;\r
2489           cv_mem.cv_tstopset = false;\r
2490           return(CV_TSTOP_RETURN);\r
2491         }\r
2492         \r
2493         /* If next step would overtake tstop, adjust stepsize */\r
2494         if ( (cv_mem.cv_tn + cv_mem.cv_hprime - cv_mem.cv_tstop)*cv_mem.cv_h > ZERO ) {\r
2495           cv_mem.cv_hprime = (cv_mem.cv_tstop - cv_mem.cv_tn)*(ONE-FOUR*cv_mem.cv_uround);\r
2496           cv_mem.cv_eta = cv_mem.cv_hprime / cv_mem.cv_h;\r
2497         }\r
2498         \r
2499       }\r
2500       \r
2501     } /* end stopping tests block at nst>0 */  \r
2502     \r
2503     /*\r
2504      * --------------------------------------------------\r
2505      * 4. Looping point for internal steps\r
2506      *\r
2507      *    4.1. check for errors (too many steps, too much\r
2508      *         accuracy requested, step size too small)\r
2509      *    4.2. take a new step (call cvStep)\r
2510      *    4.3. stop on error \r
2511      *    4.4. perform stop tests:\r
2512      *         - check for root in last step\r
2513      *         - check if tout was passed\r
2514      *         - check if close to tstop\r
2515      *         - check if in ONE_STEP mode (must return)\r
2516      * --------------------------------------------------\r
2517      */  \r
2518     \r
2519     nstloc = 0;\r
2520     for(;;) {\r
2521      \r
2522       cv_mem.cv_next_h = cv_mem.cv_h;\r
2523       cv_mem.cv_next_q = cv_mem.cv_q;\r
2524       \r
2525       /* Reset and check ewt, ewtQ, ewtS */   \r
2526       if (cv_mem.cv_nst > 0) {\r
2527 \r
2528         ier = cv_efun(cv_mem.cv_zn[0], cv_mem.cv_ewt, cv_mem.cv_e_data, cv_mem.cv_efun_select);\r
2529         if(ier != 0) {\r
2530           if (cv_mem.cv_itol == CV_WF)\r
2531             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2532                            MSGCV_EWT_NOW_FAIL, cv_mem.cv_tn);\r
2533           else\r
2534             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2535                            MSGCV_EWT_NOW_BAD, cv_mem.cv_tn);\r
2536           istate = CV_ILL_INPUT;\r
2537           cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2538           N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
2539           break;\r
2540         }\r
2541 \r
2542         if (cv_mem.cv_quadr && cv_mem.cv_errconQ) {\r
2543           ier = cvQuadEwtSet(cv_mem, cv_mem.cv_znQ[0], cv_mem.cv_ewtQ);\r
2544           if(ier != 0) {\r
2545             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2546                            MSGCV_EWTQ_NOW_BAD, cv_mem.cv_tn);\r
2547             istate = CV_ILL_INPUT;\r
2548             cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2549             N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
2550             break;\r
2551           }\r
2552         }\r
2553 \r
2554         if (cv_mem.cv_sensi) {\r
2555           ier = cvSensEwtSet(cv_mem, cv_mem.cv_znS[0], cv_mem.cv_ewtS);\r
2556           if (ier != 0) {\r
2557             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2558                            MSGCV_EWTS_NOW_BAD, cv_mem.cv_tn);\r
2559             istate = CV_ILL_INPUT;\r
2560             cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2561             N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
2562             break;\r
2563           }\r
2564         }\r
2565 \r
2566         if (cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS) {\r
2567           ier = cvQuadSensEwtSet(cv_mem, cv_mem.cv_znQS[0], cv_mem.cv_ewtQS);\r
2568           if (ier != 0) {\r
2569             cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
2570                            MSGCV_EWTQS_NOW_BAD, cv_mem.cv_tn);\r
2571             istate = CV_ILL_INPUT;\r
2572             cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2573             N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
2574             break;\r
2575           }\r
2576         }\r
2577 \r
2578       }\r
2579 \r
2580       /* Check for too many steps */\r
2581       if ( (cv_mem.cv_mxstep>0) && (nstloc >= cv_mem.cv_mxstep) ) {\r
2582         cvProcessError(cv_mem, CV_TOO_MUCH_WORK, "CVODES", "CVode",\r
2583                        MSGCV_MAX_STEPS, cv_mem.cv_tn);\r
2584         istate = CV_TOO_MUCH_WORK;\r
2585         cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2586         N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
2587         break;\r
2588       }\r
2589 \r
2590       /* Check for too much accuracy requested */\r
2591       nrm = N_VWrmsNorm_Serial(cv_mem.cv_zn[0], cv_mem.cv_ewt);\r
2592       if (cv_mem.cv_quadr && cv_mem.cv_errconQ) {\r
2593         nrm = cvQuadUpdateNorm(cv_mem, nrm, cv_mem.cv_znQ[0], cv_mem.cv_ewtQ); \r
2594       }\r
2595       if (cv_mem.cv_sensi && cv_mem.cv_errconS) {\r
2596         nrm = cvSensUpdateNorm(cv_mem, nrm, cv_mem.cv_znS[0], cv_mem.cv_ewtS);\r
2597       }\r
2598       if (cv_mem.cv_quadr_sensi && cv_mem.cv_errconQS) {\r
2599         nrm = cvQuadSensUpdateNorm(cv_mem, nrm, cv_mem.cv_znQS[0], cv_mem.cv_ewtQS);\r
2600       }\r
2601       cv_mem.cv_tolsf = cv_mem.cv_uround * nrm;\r
2602       if (cv_mem.cv_tolsf > ONE) {\r
2603         cvProcessError(cv_mem, CV_TOO_MUCH_ACC, "CVODES", "CVode",\r
2604                        MSGCV_TOO_MUCH_ACC, cv_mem.cv_tn);\r
2605         istate = CV_TOO_MUCH_ACC;\r
2606         cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2607         N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
2608         cv_mem.cv_tolsf *= TWO;\r
2609         break;\r
2610       } else {\r
2611         cv_mem.cv_tolsf = ONE;\r
2612       }\r
2613       \r
2614       /* Check for h below roundoff level in tn */\r
2615       if (cv_mem.cv_tn + cv_mem.cv_h == cv_mem.cv_tn) {\r
2616         cv_mem.cv_nhnil++;\r
2617         if (cv_mem.cv_nhnil <= cv_mem.cv_mxhnil) \r
2618           cvProcessError(cv_mem, CV_WARNING, "CVODES", "CVode", MSGCV_HNIL,\r
2619                          cv_mem.cv_tn, cv_mem.cv_h);\r
2620         if (cv_mem.cv_nhnil == cv_mem.cv_mxhnil) \r
2621           cvProcessError(cv_mem, CV_WARNING, "CVODES", "CVode", MSGCV_HNIL_DONE);\r
2622       }\r
2623 \r
2624       /* Call cvStep to take a step */\r
2625       kflag = cvStep(cv_mem);\r
2626 \r
2627       /* Process failed step cases, and exit loop */\r
2628       if (kflag != CV_SUCCESS) {\r
2629         istate = cvHandleFailure(cv_mem, kflag);\r
2630         cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2631         N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
2632         break;\r
2633       }\r
2634       \r
2635       nstloc++;\r
2636 \r
2637       /* If tstop is set and was reached, reset tn = tstop */\r
2638       if ( cv_mem.cv_tstopset ) {\r
2639         troundoff = FUZZ_FACTOR * cv_mem.cv_uround *\r
2640           (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_h));\r
2641         if ( Math.abs(cv_mem.cv_tn - cv_mem.cv_tstop) <= troundoff)\r
2642           cv_mem.cv_tn = cv_mem.cv_tstop;\r
2643       }\r
2644 \r
2645       /* Check for root in last step taken. */    \r
2646       if (cv_mem.cv_nrtfn > 0) {\r
2647         \r
2648         retval = cvRcheck3(cv_mem);\r
2649         \r
2650         if (retval == RTFOUND) {  /* A new root was found */\r
2651           cv_mem.cv_irfnd = 1;\r
2652           istate = CV_ROOT_RETURN;\r
2653           cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tlo;\r
2654           break;\r
2655         } else if (retval == CV_RTFUNC_FAIL) { /* g failed */\r
2656           cvProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODES", "cvRcheck3",\r
2657                          MSGCV_RTFUNC_FAILED, cv_mem.cv_tlo);\r
2658           istate = CV_RTFUNC_FAIL;\r
2659           break;\r
2660         }\r
2661 \r
2662         /* If we are at the end of the first step and we still have\r
2663          * some event functions that are inactive, issue a warning\r
2664          * as this may indicate a user error in the implementation\r
2665          * of the root function. */\r
2666 \r
2667         if (cv_mem.cv_nst==1) {\r
2668           inactive_roots = false;\r
2669           for (ir=0; ir<cv_mem.cv_nrtfn; ir++) { \r
2670             if (!cv_mem.cv_gactive[ir]) {\r
2671               inactive_roots = true;\r
2672               break;\r
2673             }\r
2674           }\r
2675           if ((cv_mem.cv_mxgnull > 0) && inactive_roots) {\r
2676             cvProcessError(cv_mem, CV_WARNING, "CVODES", "CVode",\r
2677                            MSGCV_INACTIVE_ROOTS);\r
2678           }\r
2679         }\r
2680 \r
2681       }\r
2682 \r
2683       /* In NORMAL mode, check if tout reached */\r
2684       if ( (itask == CV_NORMAL) &&  (cv_mem.cv_tn-tout)*cv_mem.cv_h >= ZERO ) {\r
2685         istate = CV_SUCCESS;\r
2686         cv_mem.cv_tretlast = tret[0] = tout;\r
2687         CVodeGetDky(cv_mem, tout, 0, yout);\r
2688         cv_mem.cv_next_q = cv_mem.cv_qprime;\r
2689         cv_mem.cv_next_h = cv_mem.cv_hprime;\r
2690         break;\r
2691       }\r
2692 \r
2693       /* Check if tn is at tstop, or about to pass tstop */\r
2694       if ( cv_mem.cv_tstopset ) {\r
2695 \r
2696         troundoff = FUZZ_FACTOR * cv_mem.cv_uround *\r
2697           (Math.abs(cv_mem.cv_tn) + Math.abs(cv_mem.cv_h));\r
2698         if ( Math.abs(cv_mem.cv_tn - cv_mem.cv_tstop) <= troundoff) {\r
2699           CVodeGetDky(cv_mem, cv_mem.cv_tstop, 0, yout);\r
2700           cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tstop;\r
2701           cv_mem.cv_tstopset = false;\r
2702           istate = CV_TSTOP_RETURN;\r
2703           break;\r
2704         }\r
2705 \r
2706         if ( (cv_mem.cv_tn + cv_mem.cv_hprime - cv_mem.cv_tstop)*cv_mem.cv_h > ZERO ) {\r
2707           cv_mem.cv_hprime = (cv_mem.cv_tstop - cv_mem.cv_tn)*(ONE-FOUR*cv_mem.cv_uround);\r
2708           cv_mem.cv_eta = cv_mem.cv_hprime / cv_mem.cv_h;\r
2709         }\r
2710 \r
2711       }\r
2712 \r
2713       /* In ONE_STEP mode, copy y and exit loop */\r
2714       if (itask == CV_ONE_STEP) {\r
2715         istate = CV_SUCCESS;\r
2716         cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
2717         N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
2718         cv_mem.cv_next_q = cv_mem.cv_qprime;\r
2719         cv_mem.cv_next_h = cv_mem.cv_hprime;\r
2720         break;\r
2721       }\r
2722 \r
2723     } /* end looping for internal steps */\r