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