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