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