fix linux compilation with --enable-stop-on-warning argument.
[scilab.git] / scilab / modules / scicos / src / scicos_sundials / src / cvode / cvode.c
1 /*
2  * -----------------------------------------------------------------
3  * $Revision: 1.4 $
4  * $Date: 2006/10/17 21:00:02 $
5  * -----------------------------------------------------------------
6  * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh, Radu Serban,
7  *                and Dan Shumaker @ LLNL
8  * -----------------------------------------------------------------
9  * Copyright (c) 2002, The Regents of the University of California.
10  * Produced at the Lawrence Livermore National Laboratory.
11  * All rights reserved.
12  * For details, see the LICENSE file.
13  * -----------------------------------------------------------------
14  * This is the implementation file for the main CVODE integrator.
15  * It is independent of the CVODE linear solver in use.
16  * -----------------------------------------------------------------
17  */
18
19 /*=================================================================*/
20 /*             Import Header Files                                 */
21 /*=================================================================*/
22
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <stdarg.h>
26
27 #include "cvode_impl.h"
28 #include <sundials/sundials_math.h>
29 #include <sundials/sundials_types.h>
30
31 /* SUNDIALS EXTENSION */
32 #include "sundials_extension.h"
33 #include "nvector/nvector_serial.h"
34
35 /*=================================================================*/
36 /*             CVODE Private Constants                             */
37 /*=================================================================*/
38
39 #define ZERO   RCONST(0.0)     /* real 0.0     */
40 #define TINY   RCONST(1.0e-10) /* small number */
41 #define TENTH  RCONST(0.1)     /* real 0.1     */
42 #define POINT2 RCONST(0.2)     /* real 0.2     */
43 #define FOURTH RCONST(0.25)    /* real 0.25    */
44 #define HALF   RCONST(0.5)     /* real 0.5     */
45 #define ONE    RCONST(1.0)     /* real 1.0     */
46 #define TWO    RCONST(2.0)     /* real 2.0     */
47 #define THREE  RCONST(3.0)     /* real 3.0     */
48 #define FOUR   RCONST(4.0)     /* real 4.0     */
49 #define FIVE   RCONST(5.0)     /* real 5.0     */
50 #define TWELVE RCONST(12.0)    /* real 12.0    */
51 #define HUN    RCONST(100.0)   /* real 100.0   */
52
53 /*=================================================================*/
54 /*             CVODE Routine-Specific Constants                   */
55 /*=================================================================*/
56
57 /* 
58  * Control constants for lower-level functions used by CVStep 
59  * ----------------------------------------------------------
60  *
61  * CVHin return values:
62  *    CV_SUCCESS
63  *    CV_RHSFUNC_FAIL
64  *    CV_TOO_CLOSE
65  *
66  * CVStep control constants:
67  *    DO_ERROR_TEST
68  *    PREDICT_AGAIN
69  *
70  * CVStep return values: 
71  *    CV_SUCCESS,
72  *    CV_LSETUP_FAIL,  CV_LSOLVE_FAIL, 
73  *    CV_RHSFUNC_FAIL, CV_RTFUNC_FAIL
74  *    CV_CONV_FAILURE, CV_ERR_FAILURE,
75  *    CV_FIRST_RHSFUNC_ERR
76  *
77  * CVNls input nflag values:
78  *    FIRST_CALL
79  *    PREV_CONV_FAIL
80  *    PREV_ERR_FAIL
81  *    
82  * CVNls return values: 
83  *    CV_SUCCESS,
84  *    CV_LSETUP_FAIL, CV_LSOLVE_FAIL, CV_RHSFUNC_FAIL,
85  *    CONV_FAIL, RHSFUNC_RECVR
86  * 
87  * CVNewtonIteration return values:
88  *    CV_SUCCESS, 
89  *    CV_LSOLVE_FAIL, CV_RHSFUNC_FAIL
90  *    CONV_FAIL, RHSFUNC_RECVR,
91  *    TRY_AGAIN
92  * 
93  */
94
95 #define DO_ERROR_TEST    +2
96 #define PREDICT_AGAIN    +3
97
98 #define CONV_FAIL        +4 
99 #define TRY_AGAIN        +5
100
101 #define FIRST_CALL       +6
102 #define PREV_CONV_FAIL   +7
103 #define PREV_ERR_FAIL    +8
104
105 #define RHSFUNC_RECVR    +9
106
107 /*
108  * Control constants for lower-level rootfinding functions
109  * -------------------------------------------------------
110  *
111  * CVRcheck1 return values:
112  *    CV_SUCCESS,
113  *    CV_RTFUNC_FAIL,
114  *    INITROOT
115  * CVRcheck2 return values:
116  *    CV_SUCCESS
117  *    CV_RTFUNC_FAIL,
118  *    CLOSERT
119  *    RTFOUND
120  * CVRcheck3 return values:
121  *    CV_SUCCESS
122  *    CV_RTFUNC_FAIL,
123  *    RTFOUND
124  * CVRootFind return values:
125  *    CV_SUCCESS
126  *    CV_RTFUNC_FAIL,
127  *    RTFOUND
128  */
129
130 #define RTFOUND          +1
131 #define INITROOT         +2
132 #define CLOSERT          +3
133
134 /* SUNDIALS EXTENSION */
135 #define ZERODETACHING    +4
136 #define MASKED           55
137
138 /*
139  * Algorithmic constants
140  * ---------------------
141  *
142  * CVodeGetDky and CVStep
143  *
144  *    FUZZ_FACTOR
145  *
146  * CVHin
147  *
148  *    HLB_FACTOR
149  *    HUB_FACTOR
150  *    H_BIAS
151  *    MAX_ITERS
152  *
153  * CVodeCreate 
154  *
155  *   CORTES
156  *
157  * CVStep
158  *
159  *    THRESH
160  *    ETAMX1
161  *    ETAMX2
162  *    ETAMX3
163  *    ETAMXF
164  *    ETAMIN
165  *    ETACF
166  *    ADDON
167  *    BIAS1
168  *    BIAS2
169  *    BIAS3
170  *    ONEPSM
171  *
172  *    SMALL_NST   nst > SMALL_NST => use ETAMX3 
173  *    MXNCF       max no. of convergence failures during one step try
174  *    MXNEF       max no. of error test failures during one step try
175  *    MXNEF1      max no. of error test failures before forcing a reduction of order
176  *    SMALL_NEF   if an error failure occurs and SMALL_NEF <= nef <= MXNEF1, then
177  *                reset eta =  MIN(eta, ETAMXF)
178  *    LONG_WAIT   number of steps to wait before considering an order change when
179  *                q==1 and MXNEF1 error test failures have occurred
180  *
181  * CVNls
182  *    
183  *    NLS_MAXCOR  maximum no. of corrector iterations for the nonlinear solver
184  *    CRDOWN      constant used in the estimation of the convergence rate (crate)
185  *                of the iterates for the nonlinear equation
186  *    DGMAX       iter == CV_NEWTON, |gamma/gammap-1| > DGMAX => call lsetup
187  *    RDIV        declare divergence if ratio del/delp > RDIV
188  *    MSBP        max no. of steps between lsetup calls
189  *    
190  */
191
192
193 #define FUZZ_FACTOR RCONST(100.0)
194
195 #define HLB_FACTOR RCONST(100.0)
196 #define HUB_FACTOR RCONST(0.1)
197 #define H_BIAS     HALF
198 #define MAX_ITERS  4
199
200 #define CORTES RCONST(0.1)
201
202 #define THRESH RCONST(1.5)
203 #define ETAMX1 RCONST(10000.0) 
204 #define ETAMX2 RCONST(10.0)
205 #define ETAMX3 RCONST(10.0)
206 #define ETAMXF RCONST(0.2)
207 #define ETAMIN RCONST(0.1)
208 #define ETACF  RCONST(0.25)
209 #define ADDON  RCONST(0.000001)
210 #define BIAS1  RCONST(6.0)
211 #define BIAS2  RCONST(6.0)
212 #define BIAS3  RCONST(10.0)
213 #define ONEPSM RCONST(1.000001)
214
215 #define SMALL_NST    10
216 #define MXNCF        10
217 #define MXNEF         7
218 #define MXNEF1        3
219 #define SMALL_NEF     2
220 #define LONG_WAIT    10
221
222 #define NLS_MAXCOR 3
223 #define CRDOWN RCONST(0.3)
224 #define DGMAX  RCONST(0.3)
225
226 #define RDIV      TWO
227 #define MSBP       20
228
229 /*=================================================================*/
230 /*             Private Helper Functions Prototypes                 */
231 /*=================================================================*/
232
233 static booleantype CVCheckNvector(N_Vector tmpl);
234
235 static int CVInitialSetup(CVodeMem cv_mem);
236
237 static booleantype CVAllocVectors(CVodeMem cv_mem, N_Vector tmpl, int tol);
238 static void CVFreeVectors(CVodeMem cv_mem);
239
240 static int CVEwtSetSS(CVodeMem cv_mem, N_Vector ycur, N_Vector weight);
241 static int CVEwtSetSV(CVodeMem cv_mem, N_Vector ycur, N_Vector weight);
242
243 static int CVHin(CVodeMem cv_mem, realtype tout);
244 static int CVHinFixed(CVodeMem cv_mem, realtype tout, realtype *tret);
245 static realtype CVUpperBoundH0(CVodeMem cv_mem, realtype tdist);
246 static int CVYddNorm(CVodeMem cv_mem, realtype hg, realtype *yddnrm);
247
248 static int CVStep(CVodeMem cv_mem);
249 static int CVStepDoPri(CVodeMem cv_mem);
250
251 static int CVsldet(CVodeMem cv_mem);
252
253 static void CVAdjustParams(CVodeMem cv_mem);
254 static void CVAdjustOrder(CVodeMem cv_mem, int deltaq);
255 static void CVAdjustAdams(CVodeMem cv_mem, int deltaq);
256 static void CVAdjustBDF(CVodeMem cv_mem, int deltaq);
257 static void CVIncreaseBDF(CVodeMem cv_mem);
258 static void CVDecreaseBDF(CVodeMem cv_mem);
259
260 static void CVRescale(CVodeMem cv_mem);
261
262 static void CVPredict(CVodeMem cv_mem);
263
264 static void CVSet(CVodeMem cv_mem);
265 static void CVSetAdams(CVodeMem cv_mem);
266 static realtype CVAdamsStart(CVodeMem cv_mem, realtype m[]);
267 static void CVAdamsFinish(CVodeMem cv_mem, realtype m[], realtype M[], realtype hsum);
268 static realtype CVAltSum(int iend, realtype a[], int k);
269 static void CVSetBDF(CVodeMem cv_mem);
270 static void CVSetTqBDF(CVodeMem cv_mem, realtype hsum, realtype alpha0,
271                        realtype alpha0_hat, realtype xi_inv, realtype xistar_inv);
272
273 static int CVNls(CVodeMem cv_mem, int nflag);
274 static int CVNlsFunctional(CVodeMem cv_mem);
275 static int CVNlsNewton(CVodeMem cv_mem, int nflag);
276 static int CVNewtonIteration(CVodeMem cv_mem);
277
278 static int CVHandleNFlag(CVodeMem cv_mem, int *nflagPtr, realtype saved_t,
279                          int *ncfPtr);
280
281 static void CVRestore(CVodeMem cv_mem, realtype saved_t);
282
283 static int CVDoErrorTest(CVodeMem cv_mem, int *nflagPtr,
284                          realtype saved_t, int *nefPtr, realtype *dsmPtr);
285
286 static void CVCompleteStep(CVodeMem cv_mem);
287
288 static void CVPrepareNextStep(CVodeMem cv_mem, realtype dsm);
289 static void CVSetEta(CVodeMem cv_mem);
290 static realtype CVComputeEtaqm1(CVodeMem cv_mem);
291 static realtype CVComputeEtaqp1(CVodeMem cv_mem);
292 static void CVChooseEta(CVodeMem cv_mem);
293 static void CVBDFStab(CVodeMem cv_mem);
294
295 static int  CVHandleFailure(CVodeMem cv_mem,int flag);
296
297 static int CVRcheck1(CVodeMem cv_mem);
298 static int CVRcheck2(CVodeMem cv_mem);
299 static int CVRcheck3(CVodeMem cv_mem);
300 static int CVRootfind(CVodeMem cv_mem);
301
302 /* SUNDIALS STANDARD */
303 static int CVRcheck1Std(CVodeMem cv_mem);
304 static int CVRcheck2Std(CVodeMem cv_mem);
305 static int CVRcheck3Std(CVodeMem cv_mem);
306 static int CVRootfindStd(CVodeMem cv_mem);
307
308 /* SUNDIALS EXTENSION */
309 static int CVRcheck1Ext(CVodeMem cv_mem);
310 static int CVRcheck2Ext(CVodeMem cv_mem);
311 static int CVRcheck3Ext(CVodeMem cv_mem);
312 static int CVRootfindExt(CVodeMem cv_mem);
313
314
315 /* 
316  * =================================================================
317  * EXPORTED FUNCTIONS IMPLEMENTATION
318  * =================================================================
319  */
320
321 /* 
322  * CVodeCreate
323  *
324  * CVodeCreate creates an internal memory block for a problem to 
325  * be solved by CVODE.
326  * If successful, CVodeCreate returns a pointer to the problem memory. 
327  * This pointer should be passed to CVodeMalloc.  
328  * If an initialization error occurs, CVodeCreate prints an error 
329  * message to standard err and returns NULL. 
330  */
331
332 void *CVodeCreate(int lmm, int iter)
333 {
334   int maxord;
335   CVodeMem cv_mem;
336
337   /* Test inputs */
338
339   if ((lmm != CV_ADAMS) && (lmm != CV_BDF) && (lmm != CV_DOPRI)) { /* Integration mode : ADAMS, BDF or DOPRI */
340     CVProcessError(NULL, 0, "CVODE", "CVodeCreate", MSGCV_BAD_LMM);
341     return(NULL);
342   }
343   
344   if ((iter != CV_FUNCTIONAL) && (iter != CV_NEWTON)) {
345     CVProcessError(NULL, 0, "CVODE", "CVodeCreate", MSGCV_BAD_ITER);
346     return(NULL);
347   }
348
349   cv_mem = NULL;
350   cv_mem = (CVodeMem) malloc(sizeof(struct CVodeMemRec));
351   if (cv_mem == NULL) {
352     CVProcessError(NULL, 0, "CVODE", "CVodeCreate", MSGCV_CVMEM_FAIL);
353     return(NULL);
354   }
355
356   maxord = (lmm == CV_ADAMS) ? ADAMS_Q_MAX : BDF_Q_MAX;
357
358   /* If Runge-Kutta is selected, then maxord = 9 to use the 7 extra vectors allocated (zn[2, 3, 4, 5, 6, 7, 8]) */
359   maxord = (lmm == CV_DOPRI) ? 8 : maxord;
360
361   /* copy input parameters into cv_mem */
362   cv_mem->cv_lmm  = lmm;
363   cv_mem->cv_iter = iter;
364
365   /* Set uround */
366   cv_mem->cv_uround = UNIT_ROUNDOFF;
367
368   /* Set default values for integrator optional inputs */
369   cv_mem->cv_f        = NULL;
370   cv_mem->cv_f_data   = NULL;
371   cv_mem->cv_efun     = NULL;
372   cv_mem->cv_e_data   = NULL;
373   cv_mem->cv_ehfun    = CVErrHandler;
374   cv_mem->cv_eh_data  = (void *) cv_mem;
375   cv_mem->cv_errfp    = stderr;
376   cv_mem->cv_qmax     = maxord;
377   cv_mem->cv_mxstep   = MXSTEP_DEFAULT;
378   cv_mem->cv_mxhnil   = MXHNIL_DEFAULT;
379   cv_mem->cv_sldeton  = FALSE;
380   cv_mem->cv_hin      = ZERO;
381   cv_mem->cv_hmin     = HMIN_DEFAULT;
382   cv_mem->cv_hmax_inv = HMAX_INV_DEFAULT;
383   cv_mem->cv_tstopset = FALSE;
384   cv_mem->cv_maxcor   = NLS_MAXCOR;
385   cv_mem->cv_maxnef   = MXNEF;
386   cv_mem->cv_maxncf   = MXNCF;
387   cv_mem->cv_nlscoef  = CORTES;
388
389   /* Initialize root finding variables */
390
391   cv_mem->cv_glo    = NULL;
392   cv_mem->cv_ghi    = NULL;
393   cv_mem->cv_grout  = NULL;
394   cv_mem->cv_iroots = NULL;
395   cv_mem->cv_gfun   = NULL;
396   cv_mem->cv_g_data = NULL;
397   cv_mem->cv_nrtfn  = 0;
398
399   /* Set the saved value qmax_alloc */
400
401   cv_mem->cv_qmax_alloc = maxord;
402   
403   /* Initialize lrw and liw */
404
405   cv_mem->cv_lrw = 58 + 2*L_MAX + NUM_TESTS;
406   cv_mem->cv_liw = 40;
407
408   /* No mallocs have been done yet */
409
410   cv_mem->cv_VabstolMallocDone = FALSE;
411   cv_mem->cv_MallocDone        = FALSE;
412
413   /* Return pointer to CVODE memory block */
414
415   return((void *)cv_mem);
416 }
417
418 /*-----------------------------------------------------------------*/
419
420 #define iter (cv_mem->cv_iter)  
421 #define lmm  (cv_mem->cv_lmm) 
422 #define lrw  (cv_mem->cv_lrw)
423 #define liw  (cv_mem->cv_liw)
424
425 /*-----------------------------------------------------------------*/
426
427 /*
428  * CVodeMalloc
429  * 
430  * CVodeMalloc allocates and initializes memory for a problem. All 
431  * problem inputs are checked for errors. If any error occurs during 
432  * initialization, it is reported to the file whose file pointer is 
433  * errfp and an error flag is returned. Otherwise, it returns CV_SUCCESS
434  */
435
436 int CVodeMalloc(void *cvode_mem, CVRhsFn f, realtype t0, N_Vector y0, 
437                 int itol, realtype reltol, void *abstol)
438 {
439   CVodeMem cv_mem;
440   booleantype nvectorOK, allocOK, neg_abstol;
441   long int lrw1, liw1;
442   int i,k;
443
444   /* Check cvode_mem */
445
446   if (cvode_mem==NULL) {
447     CVProcessError(NULL, CV_MEM_NULL, "CVODE", "CVodeMalloc", MSGCV_NO_MEM);
448     return(CV_MEM_NULL);
449   }
450   cv_mem = (CVodeMem) cvode_mem;
451
452   /* Check for legal input parameters */
453
454   if (y0==NULL) {
455     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeMalloc", MSGCV_NULL_Y0);
456
457         /* SUNDIALS EXTENSION */
458         if (is_sundials_with_extension())
459         {
460                 return(CV_NULL_Y0);
461         }
462         else
463         {
464                 return(CV_ILL_INPUT);
465         }
466   }
467
468   if ((itol != CV_SS) && (itol != CV_SV) && (itol != CV_WF)) {
469     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeMalloc", MSGCV_BAD_ITOL);
470
471         /* SUNDIALS EXTENSION */
472         if (is_sundials_with_extension())
473         {
474                 return(CV_BAD_ITOL);
475         }
476         else
477         {
478                 return(CV_ILL_INPUT);
479         }
480
481   }
482
483   if (f == NULL) {
484     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeMalloc", MSGCV_NULL_F);
485
486         /* SUNDIALS EXTENSION */
487         if (is_sundials_with_extension())
488         {
489                 return(CV_NULL_F);
490         }
491         else
492         {
493                 return(CV_ILL_INPUT);
494         }
495
496   }
497
498   /* Test if all required vector operations are implemented */
499
500   nvectorOK = CVCheckNvector(y0);
501   if(!nvectorOK) {
502     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeMalloc", MSGCV_BAD_NVECTOR);
503
504         /* SUNDIALS EXTENSION */
505         if (is_sundials_with_extension())
506         {
507                 return(CV_BAD_NVECTOR);
508         }
509         else
510         {
511                 return(CV_ILL_INPUT);
512         }
513   }
514
515   /* Test tolerances */
516
517   if (itol != CV_WF) {
518
519     if (abstol == NULL) {
520       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeMalloc", MSGCV_NULL_ABSTOL);
521
522           /* SUNDIALS EXTENSION */
523           if (is_sundials_with_extension())
524           {
525                   return(CV_NULL_ABSTOL);
526           }
527           else
528           {
529                 return(CV_ILL_INPUT);
530           }
531
532     }
533
534     if (reltol < ZERO) {
535       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeMalloc", MSGCV_BAD_RELTOL);
536
537           /* SUNDIALS EXTENSION */
538           if (is_sundials_with_extension())
539           {
540                 return(CV_BAD_RELTOL);
541           }
542           else
543           {
544                 return(CV_ILL_INPUT);
545           }
546
547     }
548
549     if (itol == CV_SS)
550       neg_abstol = (*((realtype *)abstol) < ZERO);
551     else
552       neg_abstol = (N_VMin((N_Vector)abstol) < ZERO);
553
554     if (neg_abstol) {
555       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeMalloc", MSGCV_BAD_ABSTOL);
556
557           /* SUNDIALS EXTENSION */
558           if (is_sundials_with_extension())
559           {
560                   return(CV_BAD_ABSTOL);
561           }
562           else
563           {
564                   return(CV_ILL_INPUT);
565           }
566     }
567
568   }
569
570   /* Set space requirements for one N_Vector */
571
572   if (y0->ops->nvspace != NULL) {
573     N_VSpace(y0, &lrw1, &liw1);
574   } else {
575     lrw1 = 0;
576     liw1 = 0;
577   }
578   cv_mem->cv_lrw1 = lrw1;
579   cv_mem->cv_liw1 = liw1;
580
581   /* Allocate the vectors (using y0 as a template) */
582
583   allocOK = CVAllocVectors(cv_mem, y0, itol);
584   if (!allocOK) {
585     CVProcessError(cv_mem, CV_MEM_FAIL, "CVODE", "CVodeMalloc", MSGCV_MEM_FAIL);
586     return(CV_MEM_FAIL);
587   }
588
589   /* Copy tolerances into memory */
590
591   cv_mem->cv_itol   = itol;
592   cv_mem->cv_reltol = reltol;      
593
594   if (itol == CV_SS)
595     cv_mem->cv_Sabstol = *((realtype *)abstol);
596   else if (itol == CV_SV)
597     N_VScale(ONE, (N_Vector)abstol, cv_mem->cv_Vabstol);
598
599   /* All error checking is complete at this point */
600
601   /* Copy the input parameters into CVODE state */
602
603   cv_mem->cv_f  = f;
604   cv_mem->cv_tn = t0;
605
606   /* Set step parameters */
607
608   cv_mem->cv_q      = 1;
609   cv_mem->cv_L      = 2;
610   cv_mem->cv_qwait  = cv_mem->cv_L;
611   cv_mem->cv_etamax = ETAMX1;
612
613   cv_mem->cv_qu    = 0;
614   cv_mem->cv_hu    = ZERO;
615   cv_mem->cv_tolsf = ONE;
616
617   /* Set the linear solver addresses to NULL.
618      (We check != NULL later, in CVode, if using CV_NEWTON.) */
619
620   cv_mem->cv_linit  = NULL;
621   cv_mem->cv_lsetup = NULL;
622   cv_mem->cv_lsolve = NULL;
623   cv_mem->cv_lfree  = NULL;
624   cv_mem->cv_lmem   = NULL;
625
626   /* Initialize zn[0] in the history array */
627
628   N_VScale(ONE, y0, cv_mem->cv_zn[0]);
629
630   /* Initialize all the counters */
631
632   cv_mem->cv_nst     = 0;
633   cv_mem->cv_nfe     = 0;
634   cv_mem->cv_ncfn    = 0;
635   cv_mem->cv_netf    = 0;
636   cv_mem->cv_nni     = 0;
637   cv_mem->cv_nsetups = 0;
638   cv_mem->cv_nhnil   = 0;
639   cv_mem->cv_nstlp   = 0;
640   cv_mem->cv_nscon   = 0;
641   cv_mem->cv_nge     = 0;
642
643   /* Initialize other integrator optional outputs */
644
645   cv_mem->cv_h0u      = ZERO;
646   cv_mem->cv_next_h   = ZERO;
647   cv_mem->cv_next_q   = 0;
648
649   /* Initialize Stablilty Limit Detection data */
650   /* NOTE: We do this even if stab lim det was not
651      turned on yet. This way, the user can turn it
652      on at any time */
653
654   cv_mem->cv_nor = 0;
655   for (i = 1; i <= 5; i++)
656     for (k = 1; k <= 3; k++) 
657       cv_mem->cv_ssdat[i-1][k-1] = ZERO;
658
659   /* Problem has been successfully initialized */
660
661   cv_mem->cv_MallocDone = TRUE;
662
663   return(CV_SUCCESS);
664 }
665
666 /*-----------------------------------------------------------------*/
667
668 #define lrw1 (cv_mem->cv_lrw1)
669 #define liw1 (cv_mem->cv_liw1)
670
671 /*-----------------------------------------------------------------*/
672
673 /*
674  * CVodeReInit
675  *
676  * CVodeReInit re-initializes CVODE's memory for a problem, assuming
677  * it has already been allocated in a prior CVodeMalloc call.
678  * All problem specification inputs are checked for errors.
679  * If any error occurs during initialization, it is reported to the
680  * file whose file pointer is errfp.
681  * The return value is CV_SUCCESS = 0 if no errors occurred, or
682  * a negative value otherwise.
683  */
684
685 int CVodeReInit(void *cvode_mem, CVRhsFn f, realtype t0, N_Vector y0, 
686                 int itol, realtype reltol, void *abstol)
687 {
688   CVodeMem cv_mem;
689   booleantype neg_abstol;
690   int i,k;
691  
692   /* Check cvode_mem */
693
694   if (cvode_mem==NULL) {
695     CVProcessError(NULL, CV_MEM_NULL, "CVODE", "CVodeReInit", MSGCV_NO_MEM);
696     return(CV_MEM_NULL);
697   }
698   cv_mem = (CVodeMem) cvode_mem;
699
700   /* Check if cvode_mem was allocated */
701
702   if (cv_mem->cv_MallocDone == FALSE) {
703     CVProcessError(cv_mem, CV_NO_MALLOC, "CVODE", "CVodeReInit", MSGCV_NO_MALLOC);
704     return(CV_NO_MALLOC);
705   }
706
707 /* Check for legal input parameters */
708
709   if (y0 == NULL) {
710     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeReInit", MSGCV_NULL_Y0);
711
712         /* SUNDIALS EXTENSION */
713         if (is_sundials_with_extension())
714         {
715                 return(CV_NULL_Y0);
716         }
717         else
718         {
719                 return(CV_ILL_INPUT);
720         }
721   }
722   
723   if ((itol != CV_SS) && (itol != CV_SV) && (itol != CV_WF)) {
724     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeReInit", MSGCV_BAD_ITOL);
725
726         /* SUNDIALS EXTENSION */
727         if (is_sundials_with_extension())
728         {
729                 return(CV_BAD_ITOL);
730         }
731         else
732         {
733                 return(CV_ILL_INPUT);
734         }
735   }
736
737   if (f == NULL) {
738     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeReInit", MSGCV_NULL_F);
739
740         /* SUNDIALS EXTENSION */
741         if (is_sundials_with_extension())
742         {
743                 return(CV_NULL_F);
744         }
745         else
746         {
747                 return(CV_ILL_INPUT);
748         }
749   }
750
751   /* Test tolerances */
752
753   if (itol != CV_WF) {
754
755     if (abstol == NULL) {
756       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeReInit", MSGCV_NULL_ABSTOL);
757
758           /* SUNDIALS EXTENSION */
759           if (is_sundials_with_extension())
760           {
761                   return(CV_NULL_ABSTOL);
762           }
763           else
764           {
765                   return(CV_ILL_INPUT);
766           }
767     }
768
769     if (reltol < ZERO) {
770       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeReInit", MSGCV_BAD_RELTOL);
771
772           /* SUNDIALS EXTENSION */
773           if (is_sundials_with_extension())
774           {
775                   return(CV_BAD_RELTOL);
776           }
777           else
778           {
779                   return(CV_ILL_INPUT);
780           }
781     }
782     
783     if (itol == CV_SS) {
784       neg_abstol = (*((realtype *)abstol) < ZERO);
785     } else {
786       neg_abstol = (N_VMin((N_Vector)abstol) < ZERO);
787     }
788     
789     if (neg_abstol) {
790       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeReInit", MSGCV_BAD_ABSTOL);
791
792           /* SUNDIALS EXTENSION */
793           if (is_sundials_with_extension())
794           {
795                 return(CV_BAD_ABSTOL);
796           }
797           else
798           {
799                 return(CV_ILL_INPUT);
800           }
801     }
802
803   }
804
805   /* Copy tolerances into memory */
806
807   if ( (itol != CV_SV) && (cv_mem->cv_VabstolMallocDone) ) {
808     N_VDestroy(cv_mem->cv_Vabstol);
809     lrw -= lrw1;
810     liw -= liw1;
811     cv_mem->cv_VabstolMallocDone = FALSE;
812   }
813
814   if ( (itol == CV_SV) && !(cv_mem->cv_VabstolMallocDone) ) {
815     cv_mem->cv_Vabstol = NULL;
816     cv_mem->cv_Vabstol = N_VClone(y0);
817     lrw += lrw1;
818     liw += liw1;
819     cv_mem->cv_VabstolMallocDone = TRUE;
820   }
821
822   cv_mem->cv_itol   = itol;
823   cv_mem->cv_reltol = reltol;    
824   
825   if (itol == CV_SS)
826     cv_mem->cv_Sabstol = *((realtype *)abstol);
827   else if (itol == CV_SV)
828     N_VScale(ONE, (N_Vector)abstol, cv_mem->cv_Vabstol);
829   
830   /* Copy the input parameters into CVODE state */
831
832   cv_mem->cv_f = f;
833   cv_mem->cv_tn = t0;
834   
835   /* Set step parameters */
836
837   cv_mem->cv_q      = 1;
838   cv_mem->cv_L      = 2;
839   cv_mem->cv_qwait  = cv_mem->cv_L;
840   cv_mem->cv_etamax = ETAMX1;
841
842   cv_mem->cv_qu    = 0;
843   cv_mem->cv_hu    = ZERO;
844   cv_mem->cv_tolsf = ONE;
845
846   /* Initialize zn[0] in the history array */
847
848   N_VScale(ONE, y0, cv_mem->cv_zn[0]);
849  
850   /* Initialize all the counters */
851
852   cv_mem->cv_nst     = 0;
853   cv_mem->cv_nfe     = 0;
854   cv_mem->cv_ncfn    = 0;
855   cv_mem->cv_netf    = 0;
856   cv_mem->cv_nni     = 0;
857   cv_mem->cv_nsetups = 0;
858   cv_mem->cv_nhnil   = 0;
859   cv_mem->cv_nstlp   = 0;
860   cv_mem->cv_nscon   = 0;
861   cv_mem->cv_nge     = 0;
862
863   /* Initialize other integrator optional outputs */
864
865   cv_mem->cv_h0u      = ZERO;
866   cv_mem->cv_next_h   = ZERO;
867   cv_mem->cv_next_q   = 0;
868
869   /* Initialize Stablilty Limit Detection data */
870
871   cv_mem->cv_nor = 0;
872   for (i = 1; i <= 5; i++)
873     for (k = 1; k <= 3; k++) 
874       cv_mem->cv_ssdat[i-1][k-1] = ZERO;
875   
876   /* Problem has been successfully re-initialized */
877
878   return(CV_SUCCESS);
879 }
880
881 /*-----------------------------------------------------------------*/
882
883 #define gfun   (cv_mem->cv_gfun)
884 #define g_data (cv_mem->cv_g_data) 
885 #define glo    (cv_mem->cv_glo)
886 #define ghi    (cv_mem->cv_ghi)
887 #define grout  (cv_mem->cv_grout)
888 #define iroots (cv_mem->cv_iroots)
889
890 /*-----------------------------------------------------------------*/
891
892 /*
893  * CVodeRootInit
894  *
895  * CVodeRootInit initializes a rootfinding problem to be solved
896  * during the integration of the ODE system.  It loads the root
897  * function pointer and the number of root functions, and allocates
898  * workspace memory.  The return value is CV_SUCCESS = 0 if no errors
899  * occurred, or a negative value otherwise.
900  */
901
902 int CVodeRootInit(void *cvode_mem, int nrtfn, CVRootFn g, void *gdata)
903 {
904   CVodeMem cv_mem;
905   int nrt;
906
907   /* Check cvode_mem pointer */
908   if (cvode_mem == NULL) {
909     CVProcessError(NULL, CV_MEM_NULL, "CVODE", "CVodeRootInit", MSGCV_NO_MEM);
910     return(CV_MEM_NULL);
911   }
912   cv_mem = (CVodeMem) cvode_mem;
913
914   nrt = (nrtfn < 0) ? 0 : nrtfn;
915
916   /* If rerunning CVodeRootInit() with a different number of root
917      functions (changing number of gfun components), then free
918      currently held memory resources */
919   if ((nrt != cv_mem->cv_nrtfn) && (cv_mem->cv_nrtfn > 0)) {
920     free(glo); glo = NULL;
921     free(ghi); ghi = NULL;
922     free(grout); grout = NULL;
923     free(iroots); iroots = NULL;
924
925     lrw -= 3* (cv_mem->cv_nrtfn);
926     liw -= cv_mem->cv_nrtfn;
927   }
928
929   /* If CVodeRootInit() was called with nrtfn == 0, then set cv_nrtfn to
930      zero and cv_gfun to NULL before returning */
931   if (nrt == 0) {
932     cv_mem->cv_nrtfn = nrt;
933     gfun = NULL;
934     g_data = NULL;
935     return(CV_SUCCESS);
936   }
937
938   /* Store user's data pointer */
939   g_data = gdata;
940
941   /* If rerunning CVodeRootInit() with the same number of root functions
942      (not changing number of gfun components), then check if the root
943      function argument has changed */
944   /* If g != NULL then return as currently reserved memory resources
945      will suffice */
946   if (nrt == cv_mem->cv_nrtfn) {
947     if (g != gfun) {
948       if (g == NULL) {
949         free(glo); glo = NULL;
950         free(ghi); ghi = NULL;
951         free(grout); grout = NULL;
952         free(iroots); iroots = NULL;
953
954         lrw -= 3*nrt;
955         liw -= nrt;
956
957         CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeRootInit", MSGCV_NULL_G);
958
959                 /* SUNDIALS EXTENSION */
960                 if (is_sundials_with_extension())
961                 {
962                         return(CV_NULL_G);
963                 }
964                 else
965                 {
966                         return(CV_ILL_INPUT);
967                 }
968       }
969       else {
970         gfun = g;
971         return(CV_SUCCESS);
972       }
973     }
974     else return(CV_SUCCESS);
975   }
976
977   /* Set variable values in CVode memory block */
978   cv_mem->cv_nrtfn = nrt;
979   if (g == NULL) {
980     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVodeRootInit", MSGCV_NULL_G);
981
982         /* SUNDIALS EXTENSION */
983         if (is_sundials_with_extension())
984         {
985                 return(CV_NULL_G);
986         }
987         else
988         {
989                 return(CV_ILL_INPUT);
990         }
991
992   }
993   else gfun = g;
994
995   /* Allocate necessary memory and return */
996   glo = NULL;
997   glo = (realtype *) malloc(nrt*sizeof(realtype));
998   if (glo == NULL) {
999     CVProcessError(cv_mem, CV_MEM_FAIL, "CVODE", "CVodeRootInit", MSGCV_MEM_FAIL);
1000     return(CV_MEM_FAIL);
1001   }
1002
1003   ghi = NULL;
1004   ghi = (realtype *) malloc(nrt*sizeof(realtype));
1005   if (ghi == NULL) {
1006     free(glo); glo = NULL;
1007     CVProcessError(cv_mem, CV_MEM_FAIL, "CVODE", "CVodeRootInit", MSGCV_MEM_FAIL);
1008     return(CV_MEM_FAIL);
1009   }
1010
1011   grout = NULL;
1012   grout = (realtype *) malloc(nrt*sizeof(realtype));
1013   if (grout == NULL) {
1014     free(glo); glo = NULL;
1015     free(ghi); ghi = NULL;
1016     CVProcessError(cv_mem, CV_MEM_FAIL, "CVODE", "CVodeRootInit", MSGCV_MEM_FAIL);
1017     return(CV_MEM_FAIL);
1018   }
1019
1020   iroots = NULL;
1021   iroots = (int *) malloc(nrt*sizeof(int));
1022   if (iroots == NULL) {
1023     free(glo); glo = NULL; 
1024     free(ghi); ghi = NULL;
1025     free(grout); grout = NULL;
1026     CVProcessError(cv_mem, CV_MEM_FAIL, "CVODE", "CVodeRootInit", MSGCV_MEM_FAIL);
1027     return(CV_MEM_FAIL);
1028   }
1029
1030   lrw += 3*nrt;
1031   liw += nrt;
1032
1033   return(CV_SUCCESS);
1034 }
1035
1036 /* 
1037  * =================================================================
1038  * Readibility Constants
1039  * =================================================================
1040  */
1041
1042 #define f              (cv_mem->cv_f)      
1043 #define f_data         (cv_mem->cv_f_data) 
1044 #define efun           (cv_mem->cv_efun)
1045 #define e_data         (cv_mem->cv_e_data) 
1046 #define qmax           (cv_mem->cv_qmax)
1047 #define mxstep         (cv_mem->cv_mxstep)
1048 #define mxhnil         (cv_mem->cv_mxhnil)
1049 #define sldeton        (cv_mem->cv_sldeton)
1050 #define hin            (cv_mem->cv_hin)
1051 #define hmin           (cv_mem->cv_hmin)
1052 #define hmax_inv       (cv_mem->cv_hmax_inv)
1053 #define istop          (cv_mem->cv_istop)
1054 #define tstop          (cv_mem->cv_tstop)
1055 #define tstopset       (cv_mem->cv_tstopset)
1056 #define maxnef         (cv_mem->cv_maxnef)
1057 #define maxncf         (cv_mem->cv_maxncf)
1058 #define maxcor         (cv_mem->cv_maxcor)
1059 #define nlscoef        (cv_mem->cv_nlscoef)
1060 #define itol           (cv_mem->cv_itol)         
1061 #define reltol         (cv_mem->cv_reltol)       
1062 #define Sabstol        (cv_mem->cv_Sabstol)
1063 #define Vabstol        (cv_mem->cv_Vabstol)
1064
1065 #define uround         (cv_mem->cv_uround)  
1066 #define zn             (cv_mem->cv_zn) 
1067 #define ewt            (cv_mem->cv_ewt)  
1068 #define y              (cv_mem->cv_y)
1069 #define acor           (cv_mem->cv_acor)
1070 #define tempv          (cv_mem->cv_tempv)
1071 #define ftemp          (cv_mem->cv_ftemp) 
1072 #define q              (cv_mem->cv_q)
1073 #define qprime         (cv_mem->cv_qprime)
1074 #define next_q         (cv_mem->cv_next_q)
1075 #define qwait          (cv_mem->cv_qwait)
1076 #define L              (cv_mem->cv_L)
1077 #define h              (cv_mem->cv_h)
1078 #define hprime         (cv_mem->cv_hprime)
1079 #define next_h         (cv_mem->cv_next_h)
1080 #define eta            (cv_mem->cv_eta) 
1081 #define etaqm1         (cv_mem->cv_etaqm1) 
1082 #define etaq           (cv_mem->cv_etaq) 
1083 #define etaqp1         (cv_mem->cv_etaqp1) 
1084 #define nscon          (cv_mem->cv_nscon)
1085 #define hscale         (cv_mem->cv_hscale)
1086 #define tn             (cv_mem->cv_tn)
1087 #define tau            (cv_mem->cv_tau)
1088 #define tq             (cv_mem->cv_tq)
1089 #define l              (cv_mem->cv_l)
1090 #define rl1            (cv_mem->cv_rl1)
1091 #define gamma          (cv_mem->cv_gamma) 
1092 #define gammap         (cv_mem->cv_gammap) 
1093 #define gamrat         (cv_mem->cv_gamrat)
1094 #define crate          (cv_mem->cv_crate)
1095 #define acnrm          (cv_mem->cv_acnrm)
1096 #define mnewt          (cv_mem->cv_mnewt)
1097 #define etamax         (cv_mem->cv_etamax)
1098 #define nst            (cv_mem->cv_nst)
1099 #define nfe            (cv_mem->cv_nfe)
1100 #define ncfn           (cv_mem->cv_ncfn)
1101 #define netf           (cv_mem->cv_netf)
1102 #define nni            (cv_mem->cv_nni)
1103 #define nsetups        (cv_mem->cv_nsetups)
1104 #define nhnil          (cv_mem->cv_nhnil)
1105 #define linit          (cv_mem->cv_linit)
1106 #define lsetup         (cv_mem->cv_lsetup)
1107 #define lsolve         (cv_mem->cv_lsolve) 
1108 #define lfree          (cv_mem->cv_lfree) 
1109 #define lmem           (cv_mem->cv_lmem) 
1110 #define qu             (cv_mem->cv_qu)          
1111 #define nstlp          (cv_mem->cv_nstlp)  
1112 #define h0u            (cv_mem->cv_h0u)
1113 #define hu             (cv_mem->cv_hu)         
1114 #define saved_tq5      (cv_mem->cv_saved_tq5)  
1115 #define indx_acor      (cv_mem->cv_indx_acor)
1116 #define jcur           (cv_mem->cv_jcur)         
1117 #define tolsf          (cv_mem->cv_tolsf)      
1118 #define setupNonNull   (cv_mem->cv_setupNonNull) 
1119 #define nor            (cv_mem->cv_nor)
1120 #define ssdat          (cv_mem->cv_ssdat)
1121
1122 #define nrtfn          (cv_mem->cv_nrtfn)
1123 #define tlo            (cv_mem->cv_tlo)
1124 #define thi            (cv_mem->cv_thi)
1125 #define tretlast       (cv_mem->cv_tretlast)
1126 #define toutc          (cv_mem->cv_toutc)
1127 #define trout          (cv_mem->cv_trout)
1128 #define ttol           (cv_mem->cv_ttol)
1129 #define taskc          (cv_mem->cv_taskc)
1130 #define irfnd          (cv_mem->cv_irfnd)
1131 #define nge            (cv_mem->cv_nge)
1132
1133 /*-----------------------------------------------------------------*/
1134
1135 /*
1136  * CVode
1137  *
1138  * This routine is the main driver of the CVODE package. 
1139  *
1140  * It integrates over a time interval defined by the user, by calling
1141  * CVStep to do internal time steps.
1142  *
1143  * The first time that CVode is called for a successfully initialized
1144  * problem, it computes a tentative initial step size h.
1145  *
1146  * CVode supports four modes, specified by itask: CV_NORMAL, CV_ONE_STEP,
1147  * CV_NORMAL_TSTOP, and CV_ONE_STEP_TSTOP.
1148  * In the CV_NORMAL mode, the solver steps until it reaches or passes tout
1149  * and then interpolates to obtain y(tout).
1150  * In the CV_ONE_STEP mode, it takes one internal step and returns.
1151  * CV_NORMAL_TSTOP and CV_ONE_STEP_TSTOP are similar to CV_NORMAL and CV_ONE_STEP,
1152  * respectively, but the integration never proceeds past tstop (which
1153  * must have been defined through a call to CVodeSetStopTime).
1154  */
1155
1156 int CVode(void *cvode_mem, realtype tout, N_Vector yout, 
1157           realtype *tret, int itask)
1158 {
1159   CVodeMem cv_mem;
1160   long int nstloc;
1161   int retval, hflag, kflag, istate, ier, task, irfndp;
1162   int ewtsetOK;
1163   realtype troundoff, rh, nrm;
1164
1165   /*
1166    * -------------------------------------
1167    * 1. Check and process inputs
1168    * -------------------------------------
1169    */
1170
1171   /* Check if cvode_mem exists */
1172   if (cvode_mem == NULL) {
1173     CVProcessError(NULL, CV_MEM_NULL, "CVODE", "CVode", MSGCV_NO_MEM);
1174     return(CV_MEM_NULL);
1175   }
1176   cv_mem = (CVodeMem) cvode_mem;
1177
1178   /* Check if cvode_mem was allocated */
1179   if (cv_mem->cv_MallocDone == FALSE) {
1180     CVProcessError(cv_mem, CV_NO_MALLOC, "CVODE", "CVode", MSGCV_NO_MALLOC);
1181     return(CV_NO_MALLOC);
1182   }
1183   
1184   /* Check for yout != NULL */
1185   if ((y = yout) == NULL) {
1186     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_YOUT_NULL);
1187
1188         /* SUNDIALS EXTENSION */
1189         if (is_sundials_with_extension())
1190         {
1191                 return(CV_YOUT_NULL);
1192         }
1193         else
1194         {
1195                 return(CV_ILL_INPUT);
1196         }
1197
1198   }
1199
1200   /* Check for tret != NULL */
1201   if (tret == NULL) {
1202     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_TRET_NULL);
1203
1204         /* SUNDIALS EXTENSION */
1205         if (is_sundials_with_extension())
1206         {
1207                 return(CV_TRET_NULL);
1208         }
1209         else
1210         {
1211                 return(CV_ILL_INPUT);
1212         }
1213   }
1214
1215   /* Check for valid itask */
1216   if ((itask != CV_NORMAL)       && 
1217       (itask != CV_ONE_STEP)     &&
1218       (itask != CV_NORMAL_TSTOP) &&
1219       (itask != CV_ONE_STEP_TSTOP) ) {
1220     CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_BAD_ITASK);
1221
1222         /* SUNDIALS EXTENSION */
1223         if (is_sundials_with_extension())
1224         {
1225                 return(CV_BAD_ITASK);
1226         }
1227         else
1228         {
1229                 return(CV_ILL_INPUT);
1230         }
1231
1232   }
1233
1234   /* Split itask into task and istop */
1235   if ((itask == CV_NORMAL_TSTOP) || (itask == CV_ONE_STEP_TSTOP)) {
1236     if ( tstopset == FALSE ) {
1237       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_NO_TSTOP);
1238
1239           /* SUNDIALS EXTENSION */
1240           if (is_sundials_with_extension())
1241           {
1242                   return(CV_NO_TSTOP);
1243           }
1244           else
1245           {
1246                   return(CV_ILL_INPUT);
1247           }
1248
1249     }
1250     istop = TRUE;
1251   } else {
1252     istop = FALSE;
1253   }
1254   if ((itask == CV_NORMAL) || (itask == CV_NORMAL_TSTOP)) {
1255     task = CV_NORMAL; toutc = tout;
1256   } else {
1257     task = CV_ONE_STEP;
1258   }
1259   taskc = task;
1260
1261   /*
1262    * ----------------------------------------
1263    * 2. Initializations performed only at
1264    *    the first step (nst=0):
1265    *    - initial setup
1266    *    - initialize Nordsieck history array
1267    *    - compute initial step size
1268    *    - check for approach to tstop
1269    *    - check for approach to a root
1270    * ----------------------------------------
1271    */
1272
1273   if (nst == 0) {
1274
1275     ier = CVInitialSetup(cv_mem);
1276     if (ier!= CV_SUCCESS) return(ier);
1277     
1278     /* Call f at (t0,y0), set zn[1] = y'(t0), 
1279        set initial h (from H0 or CVHin), and scale zn[1] by h.
1280        Also check for zeros of root function g at and near t0.    */
1281     
1282     retval = f(tn, zn[0], zn[1], f_data); 
1283     nfe++;
1284     if (retval < 0) {
1285       CVProcessError(cv_mem, CV_RHSFUNC_FAIL, "CVODE", "CVode", MSGCV_RHSFUNC_FAILED, tn);
1286       return(CV_RHSFUNC_FAIL);
1287     }
1288     if (retval > 0) {
1289       CVProcessError(cv_mem, CV_FIRST_RHSFUNC_ERR, "CVODE", "CVode", MSGCV_RHSFUNC_FIRST);
1290       return(CV_FIRST_RHSFUNC_ERR);
1291     }
1292
1293     /* Set initial h (from H0 or CVHin). 
1294        If Dormand-Price is selected, then we choose to set h to hmax (= 1/hmax_inv). */
1295
1296     if ((lmm == CV_ADAMS) || (lmm == CV_BDF)) {
1297     h = hin;
1298     if ( (h != ZERO) && ((tout-tn)*h < ZERO) ) {
1299       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_BAD_H0);
1300
1301           /* SUNDIALS EXTENSION */
1302           if (is_sundials_with_extension())
1303           {
1304                   return(CV_BAD_H0);
1305           }
1306           else
1307           {
1308                   return(CV_ILL_INPUT);
1309           }
1310
1311     }
1312     if (h == ZERO) {
1313       hflag = CVHin(cv_mem, tout);
1314       if (hflag != CV_SUCCESS) {
1315         istate = CVHandleFailure(cv_mem, hflag);
1316         return(istate);
1317       }
1318     }
1319     rh = ABS(h)*hmax_inv;
1320     if (rh > ONE) h /= rh;
1321     if (ABS(h) < hmin) h *= hmin/ABS(h);
1322         }
1323
1324     else { /* Compute the fixed step size h, and set the max number of steps */
1325       mxstep = CVHinFixed(cv_mem, tout, tret);
1326     }
1327
1328     /* Check for approach to tstop */
1329
1330     if (istop) {
1331       if ( (tstop - tn)*h < ZERO ) {
1332         CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_BAD_TSTOP, tn);
1333
1334                 /* SUNDIALS EXTENSION */
1335                 if (is_sundials_with_extension())
1336                 {
1337                         return(CV_BAD_TSTOP);
1338                 }
1339                 else
1340                 {
1341                         return(CV_ILL_INPUT);
1342                 }
1343
1344       }
1345       if ( (tn + h - tstop)*h > ZERO ) 
1346         h = (tstop - tn)*(ONE-FOUR*uround);
1347     }
1348
1349     /* Scale zn[1] by h.*/
1350
1351     hscale = h; 
1352     h0u    = h;
1353     hprime = h;
1354
1355     N_VScale(h, zn[1], zn[1]);
1356
1357     /* Check for zeros of root function g at and near t0. */
1358
1359     if (nrtfn > 0) {
1360
1361       retval = CVRcheck1(cv_mem);
1362
1363       if (retval == INITROOT) {
1364         CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVRcheck1", MSGCV_BAD_INIT_ROOT);
1365
1366                 /* SUNDIALS EXTENSION */
1367                 if (is_sundials_with_extension())
1368                 {
1369                         return(CV_BAD_INIT_ROOT);
1370                 }
1371                 else
1372                 {
1373                         return(CV_ILL_INPUT);
1374                 }
1375
1376       } else if (retval == CV_RTFUNC_FAIL) {
1377         CVProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODE", "CVRcheck1", MSGCV_RTFUNC_FAILED, tn);
1378         return(CV_RTFUNC_FAIL);
1379       }
1380
1381     }
1382
1383   } /* end of first call block */
1384
1385   /*
1386    * ------------------------------------------------------
1387    * 3. At following steps, perform stop tests:
1388    *    - check for root in last step
1389    *    - check if we passed tstop
1390    *    - check if we passed tout (NORMAL mode)
1391    *    - check if current tn was returned (ONE_STEP mode)
1392    *    - check if we are close to tstop
1393    *      (adjust step size if needed)
1394    * -------------------------------------------------------
1395    */
1396
1397   if (nst > 0) {
1398
1399     /* Estimate an infinitesimal time interval to be used as
1400        a roundoff for time quantities (based on current time 
1401        and step size) */
1402     troundoff = FUZZ_FACTOR*uround*(ABS(tn) + ABS(h));
1403
1404     /* First, check for a root in the last step taken, other than the
1405        last root found, if any.  If task = CV_ONE_STEP and y(tn) was not
1406        returned because of an intervening root, return y(tn) now.     */
1407     if (nrtfn > 0) {
1408
1409       irfndp = irfnd;
1410       
1411       retval = CVRcheck2(cv_mem);
1412
1413       if (retval == CLOSERT) {
1414         CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVRcheck2", MSGCV_CLOSE_ROOTS, tlo);
1415
1416                 /* SUNDIALS EXTENSION */
1417                 if (is_sundials_with_extension())
1418                 {
1419                         return(CV_CLOSE_ROOTS);
1420                 }
1421                 else
1422                 {
1423                         return(CV_ILL_INPUT);
1424                 }
1425
1426       } else if (retval == CV_RTFUNC_FAIL) {
1427         CVProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODE", "CVRcheck2", MSGCV_RTFUNC_FAILED, tlo);
1428         return(CV_RTFUNC_FAIL);
1429       } else if (retval == RTFOUND) {
1430         tretlast = *tret = tlo;
1431         return(CV_ROOT_RETURN);
1432       }
1433
1434       /* If tn is distinct from tretlast (within roundoff),
1435          check remaining interval for roots */
1436       if ( ABS(tn - tretlast) > troundoff ) {
1437
1438         retval = CVRcheck3(cv_mem);
1439
1440         if (retval == CV_SUCCESS) {     /* no root found */
1441           irfnd = 0;
1442           if ((irfndp == 1) && (task == CV_ONE_STEP)) {
1443             tretlast = *tret = tn;
1444             N_VScale(ONE, zn[0], yout);
1445             return(CV_SUCCESS);
1446           }
1447         } else if (retval == RTFOUND) {  /* a new root was found */
1448           irfnd = 1;
1449           tretlast = *tret = tlo;
1450           return(CV_ROOT_RETURN);
1451         } else if (retval == CV_RTFUNC_FAIL) {  /* g failed */
1452           CVProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODE", "CVRcheck3", MSGCV_RTFUNC_FAILED, tlo);
1453           return(CV_RTFUNC_FAIL);
1454         }
1455
1456                 /* SUNDIALS EXTENSION */
1457                 if (is_sundials_with_extension()) 
1458                 {
1459                 if (retval == ZERODETACHING) {  /* Zero detaching */
1460            irfnd = 1;
1461            tretlast = *tret = tlo;
1462           return(CV_ZERO_DETACH_RETURN);
1463         }
1464                 }
1465       }
1466
1467     } /* end of root stop check */
1468     /* Test for tn past tstop */
1469     if ( istop && ((tstop - tn)*h < ZERO) ) {
1470       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_BAD_TSTOP, tn);
1471
1472           /* SUNDIALS EXTENSION */
1473           if (is_sundials_with_extension())
1474           {
1475                   return(CV_BAD_TSTOP);
1476           }
1477           else
1478           {
1479                   return(CV_ILL_INPUT);
1480           }
1481
1482     }
1483
1484     /* In CV_NORMAL mode, test if tout was reached */
1485     if ( (task == CV_NORMAL) && ((tn-tout)*h >= ZERO) ) {
1486       tretlast = *tret = tout;
1487       ier =  CVodeGetDky(cv_mem, tout, 0, yout);
1488       if (ier != CV_SUCCESS) {
1489         CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_BAD_TOUT, tout);
1490
1491                 /* SUNDIALS EXTENSION */
1492                 if (is_sundials_with_extension())
1493                 {
1494                         return(CV_BAD_TOUT);
1495                 }
1496                 else
1497                 {
1498                         return(CV_ILL_INPUT);
1499                 }
1500       }
1501       return(CV_SUCCESS);
1502     }
1503
1504     /* In CV_ONE_STEP mode, test if tn was returned */
1505     if ( task == CV_ONE_STEP && ABS(tn - tretlast) > troundoff ) {
1506       tretlast = *tret = tn;
1507       N_VScale(ONE, zn[0], yout);
1508       return(CV_SUCCESS);
1509     }
1510
1511     /* Test for tn at tstop or near tstop */
1512     if ( istop ) {
1513
1514       if ( ABS(tn - tstop) <= troundoff) {
1515         ier =  CVodeGetDky(cv_mem, tstop, 0, yout);
1516         if (ier != CV_SUCCESS) {
1517           CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_BAD_TSTOP, tn);
1518
1519                   /* SUNDIALS EXTENSION */
1520                   if (is_sundials_with_extension())
1521                   {
1522                         return(CV_BAD_TSTOP);
1523                   }
1524                   else
1525                   {
1526                         return(CV_ILL_INPUT);
1527                   }
1528
1529         }
1530         tretlast = *tret = tstop;
1531         return(CV_TSTOP_RETURN);
1532       }
1533       
1534       /* If next step would overtake tstop, adjust stepsize */
1535       if ( (tn + hprime - tstop)*h > ZERO ) {
1536         hprime = (tstop - tn)*(ONE-FOUR*uround);
1537         eta = hprime/h;
1538       }
1539
1540     } /* end of istop tests block */
1541     if (lmm==CV_DOPRI) mxstep = CVHinFixed(cv_mem, tout, tret);
1542   } /* end stopping tests block */  
1543
1544   /*
1545    * --------------------------------------------------
1546    * 4. Looping point for internal steps
1547    *
1548    *    4.1. check for errors (too many steps, too much
1549    *         accuracy requested, step size too small)
1550    *    4.2. take a new step (call CVStep)
1551    *    4.3. stop on error 
1552    *    4.4. perform stop tests:
1553    *         - check for root in last step
1554    *         - check if tout was passed
1555    *         - check if close to tstop
1556    *         - check if in ONE_STEP mode (must return)
1557    * --------------------------------------------------
1558    */
1559
1560   nstloc = 0;
1561   for(;;) {
1562    
1563     next_h = h;
1564     next_q = q;
1565     
1566     /* Reset and check ewt */
1567     if (nst > 0) {
1568
1569       ewtsetOK = efun(zn[0], ewt, e_data);
1570
1571       if (ewtsetOK != 0) {
1572
1573         if (itol == CV_WF)
1574           CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_EWT_NOW_FAIL, tn);
1575         else
1576           CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVode", MSGCV_EWT_NOW_BAD, tn);
1577         
1578         istate = CV_ILL_INPUT;
1579         tretlast = *tret = tn;
1580         N_VScale(ONE, zn[0], yout);
1581         break;
1582
1583       }
1584     }
1585     
1586     /* Check for too many steps */
1587     if (nstloc >= mxstep) {
1588       CVProcessError(cv_mem, CV_TOO_MUCH_WORK, "CVODE", "CVode", MSGCV_MAX_STEPS, tn);
1589       istate = CV_TOO_MUCH_WORK;
1590       tretlast = *tret = tn;
1591       N_VScale(ONE, zn[0], yout);
1592       break;
1593     }
1594
1595     /* Check for too much accuracy requested */
1596     nrm = N_VWrmsNorm(zn[0], ewt);
1597     tolsf = uround * nrm;
1598     if (tolsf > ONE) {
1599       CVProcessError(cv_mem, CV_TOO_MUCH_ACC, "CVODE", "CVode", MSGCV_TOO_MUCH_ACC, tn);
1600       istate = CV_TOO_MUCH_ACC;
1601       tretlast = *tret = tn;
1602       N_VScale(ONE, zn[0], yout);
1603       tolsf *= TWO;
1604       break;
1605     } else {
1606       tolsf = ONE;
1607     }
1608
1609     /* Check for h below roundoff level in tn */
1610     if (tn + h == tn) {
1611       nhnil++;
1612       if (nhnil <= mxhnil) 
1613         CVProcessError(cv_mem, CV_WARNING, "CVODE", "CVode", MSGCV_HNIL, tn, h);
1614       if (nhnil == mxhnil) 
1615         CVProcessError(cv_mem, CV_WARNING, "CVODE", "CVode", MSGCV_HNIL_DONE);
1616     }
1617
1618     /* Call CVStep to take a step */
1619     if (lmm == CV_DOPRI)
1620     kflag = CVStepDoPri(cv_mem);
1621         else
1622         kflag = CVStep(cv_mem);
1623
1624     /* Process failed step cases, and exit loop */
1625     if (kflag != CV_SUCCESS) {
1626       istate = CVHandleFailure(cv_mem, kflag);
1627       tretlast = *tret = tn;
1628       N_VScale(ONE, zn[0], yout);
1629       break;
1630     }
1631     
1632     nstloc++;
1633
1634     /* Check for root in last step taken. */
1635     if (nrtfn > 0) {
1636
1637       retval = CVRcheck3(cv_mem);
1638
1639       if (retval == RTFOUND) {  /* A new root was found */
1640         irfnd = 1;
1641         istate = CV_ROOT_RETURN;
1642         tretlast = *tret = tlo;
1643         break;
1644       } else if (retval == CV_RTFUNC_FAIL) { /* g failed */
1645         CVProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODE", "CVRcheck3", MSGCV_RTFUNC_FAILED, tlo);
1646         istate = CV_RTFUNC_FAIL;
1647         break;
1648       }
1649
1650           /* SUNDIALS EXTENSION */
1651           if (is_sundials_with_extension()) 
1652           {
1653           if (retval == ZERODETACHING) {  /* Zero detaching */
1654         irfnd = 1;
1655         istate = CV_ZERO_DETACH_RETURN;
1656         tretlast = *tret = tlo;
1657         break;
1658       }
1659           }
1660
1661     }
1662
1663     /* In NORMAL mode, check if tout reached */
1664     if ( (task == CV_NORMAL) &&  (tn-tout)*h >= ZERO ) {
1665       istate = CV_SUCCESS;
1666       tretlast = *tret = tout;
1667       (void) CVodeGetDky(cv_mem, tout, 0, yout);
1668       next_q = qprime;
1669       next_h = hprime;
1670       break;
1671     }
1672
1673     /* Check if tn is at tstop or near tstop */
1674     if ( istop ) {
1675
1676       troundoff = FUZZ_FACTOR*uround*(ABS(tn) + ABS(h));
1677       if ( ABS(tn - tstop) <= troundoff) {
1678         (void) CVodeGetDky(cv_mem, tstop, 0, yout);
1679         tretlast = *tret = tstop;
1680         istate = CV_TSTOP_RETURN;
1681         break;
1682       }
1683
1684       if ( (tn + hprime - tstop)*h > ZERO ) {
1685         hprime = (tstop - tn)*(ONE-FOUR*uround);
1686         eta = hprime/h;
1687       }
1688
1689     }
1690
1691     /* In ONE_STEP mode, copy y and exit loop */
1692     if (task == CV_ONE_STEP) {
1693       istate = CV_SUCCESS;
1694       tretlast = *tret = tn;
1695       N_VScale(ONE, zn[0], yout);
1696       next_q = qprime;
1697       next_h = hprime;
1698       break;
1699     }
1700
1701   } /* end looping for internal steps */
1702
1703   return(istate);
1704 }
1705
1706 /*-----------------------------------------------------------------*/
1707
1708 /*
1709  * CVodeGetDky
1710  *
1711  * This routine computes the k-th derivative of the interpolating
1712  * polynomial at the time t and stores the result in the vector dky.
1713  * The formula is:
1714  *         q 
1715  *  dky = SUM c(j,k) * (t - tn)^(j-k) * h^(-j) * zn[j] , 
1716  *        j=k 
1717  * where c(j,k) = j*(j-1)*...*(j-k+1), q is the current order, and
1718  * zn[j] is the j-th column of the Nordsieck history array.
1719  *
1720  * This function is called by CVode with k = 0 and t = tout, but
1721  * may also be called directly by the user.
1722  */
1723
1724 int CVodeGetDky(void *cvode_mem, realtype t, int k, N_Vector dky)
1725 {
1726   realtype s, c, r;
1727   realtype tfuzz, tp, tn1;
1728   int i, j;
1729   CVodeMem cv_mem;
1730   
1731   /* Check all inputs for legality */
1732  
1733   if (cvode_mem == NULL) {
1734     CVProcessError(NULL, CV_MEM_NULL, "CVODE", "CVodeGetDky", MSGCV_NO_MEM);
1735     return(CV_MEM_NULL);
1736   }
1737   cv_mem = (CVodeMem) cvode_mem;
1738
1739   if (dky == NULL) {
1740     CVProcessError(cv_mem, CV_BAD_DKY, "CVODE", "CVodeGetDky", MSGCV_NULL_DKY);
1741     return(CV_BAD_DKY);
1742   }
1743
1744   if ((k < 0) || (k > q)) {
1745     CVProcessError(cv_mem, CV_BAD_K, "CVODE", "CVodeGetDky", MSGCV_BAD_K);
1746     return(CV_BAD_K);
1747   }
1748   
1749   /* Allow for some slack */
1750   tfuzz = FUZZ_FACTOR * uround * (ABS(tn) + ABS(hu));
1751   if (hu < ZERO) tfuzz = -tfuzz;
1752   tp = tn - hu - tfuzz;
1753   tn1 = tn + tfuzz;
1754   if ((t-tp)*(t-tn1) > ZERO) {
1755     CVProcessError(cv_mem, CV_BAD_T, "CVODE", "CVodeGetDky", MSGCV_BAD_T, t, tn-hu, tn);
1756     return(CV_BAD_T);
1757   }
1758
1759   /* Sum the differentiated interpolating polynomial */
1760
1761   s = (t - tn) / h;
1762   for (j=q; j >= k; j--) {
1763     c = ONE;
1764     for (i=j; i >= j-k+1; i--) c *= i;
1765     if (j == q) {
1766       N_VScale(c, zn[q], dky);
1767     } else {
1768       N_VLinearSum(c, zn[j], s, dky, dky);
1769     }
1770   }
1771   if (k == 0) return(CV_SUCCESS);
1772   r = RPowerI(h,-k);
1773   N_VScale(r, dky, dky);
1774   return(CV_SUCCESS);
1775 }
1776
1777 /*
1778  * CVodeFree
1779  *
1780  * This routine frees the problem memory allocated by CVodeMalloc.
1781  * Such memory includes all the vectors allocated by CVAllocVectors,
1782  * and the memory lmem for the linear solver (deallocated by a call
1783  * to lfree).
1784  */
1785
1786 void CVodeFree(void **cvode_mem)
1787 {
1788   CVodeMem cv_mem;
1789
1790   if (*cvode_mem == NULL) return;
1791
1792   cv_mem = (CVodeMem) (*cvode_mem);
1793   
1794   CVFreeVectors(cv_mem);
1795
1796   if (iter == CV_NEWTON && lfree != NULL) lfree(cv_mem);
1797
1798   if (nrtfn > 0) {
1799     free(glo); glo = NULL;
1800     free(ghi); ghi = NULL;
1801     free(grout); grout = NULL;
1802     free(iroots); iroots = NULL;
1803   }
1804
1805   free(*cvode_mem);
1806   *cvode_mem = NULL;
1807 }
1808
1809 /* 
1810  * =================================================================
1811  *  Private Functions Implementation
1812  * =================================================================
1813  */
1814
1815 /*
1816  * CVCheckNvector
1817  * This routine checks if all required vector operations are present.
1818  * If any of them is missing it returns FALSE.
1819  */
1820
1821 static booleantype CVCheckNvector(N_Vector tmpl)
1822 {
1823   if((tmpl->ops->nvclone     == NULL) ||
1824      (tmpl->ops->nvdestroy   == NULL) ||
1825      (tmpl->ops->nvlinearsum == NULL) ||
1826      (tmpl->ops->nvconst     == NULL) ||
1827      (tmpl->ops->nvprod      == NULL) ||
1828      (tmpl->ops->nvdiv       == NULL) ||
1829      (tmpl->ops->nvscale     == NULL) ||
1830      (tmpl->ops->nvabs       == NULL) ||
1831      (tmpl->ops->nvinv       == NULL) ||
1832      (tmpl->ops->nvaddconst  == NULL) ||
1833      (tmpl->ops->nvmaxnorm   == NULL) ||
1834      (tmpl->ops->nvwrmsnorm  == NULL) ||
1835      (tmpl->ops->nvmin       == NULL))
1836     return(FALSE);
1837   else
1838     return(TRUE);
1839 }
1840
1841 /*
1842  * CVAllocVectors
1843  *
1844  * This routine allocates the CVODE vectors ewt, acor, tempv, ftemp, and
1845  * zn[0], ..., zn[maxord]. If itol=CV_SV, it also allocates space for Vabstol.
1846  * If all memory allocations are successful, CVAllocVectors returns TRUE. 
1847  * Otherwise all allocated memory is freed and CVAllocVectors returns FALSE.
1848  * This routine also sets the optional outputs lrw and liw, which are
1849  * (respectively) the lengths of the real and integer work spaces
1850  * allocated here.
1851  */
1852
1853 static booleantype CVAllocVectors(CVodeMem cv_mem, N_Vector tmpl, int tol)
1854 {
1855   int i, j;
1856
1857   /* Allocate ewt, acor, tempv, ftemp */
1858   
1859   ewt = NULL;
1860   ewt = N_VClone(tmpl);
1861   if (ewt == NULL) return(FALSE);
1862
1863   acor = NULL;
1864   acor = N_VClone(tmpl);
1865   if (acor == NULL) {
1866     N_VDestroy(ewt);
1867     return(FALSE);
1868   }
1869
1870   tempv = NULL;
1871   tempv = N_VClone(tmpl);
1872   if (tempv == NULL) {
1873     N_VDestroy(ewt);
1874     N_VDestroy(acor);
1875     return(FALSE);
1876   }
1877
1878   ftemp = NULL;
1879   ftemp = N_VClone(tmpl);
1880   if (ftemp == NULL) {
1881     N_VDestroy(tempv);
1882     N_VDestroy(ewt);
1883     N_VDestroy(acor);
1884     return(FALSE);
1885   }
1886
1887   /* Allocate zn[0] ... zn[qmax] */
1888
1889   for (j=0; j <= qmax; j++) {
1890     zn[j] = NULL;
1891     zn[j] = N_VClone(tmpl);
1892     if (zn[j] == NULL) {
1893       N_VDestroy(ewt);
1894       N_VDestroy(acor);
1895       N_VDestroy(tempv);
1896       N_VDestroy(ftemp);
1897       for (i=0; i < j; i++) N_VDestroy(zn[i]);
1898       return(FALSE);
1899     }
1900   }
1901
1902   /* Update solver workspace lengths  */
1903   lrw += (qmax + 5)*lrw1;
1904   liw += (qmax + 5)*liw1;
1905
1906   if (tol == CV_SV) {
1907     Vabstol = NULL;
1908     Vabstol = N_VClone(tmpl);
1909     if (Vabstol == NULL) {
1910       N_VDestroy(ewt);
1911       N_VDestroy(acor);
1912       N_VDestroy(tempv);
1913       N_VDestroy(ftemp);
1914       for (i=0; i <= qmax; i++) N_VDestroy(zn[i]);
1915       return(FALSE);
1916     }
1917     lrw += lrw1;
1918     liw += liw1;
1919     cv_mem->cv_VabstolMallocDone = TRUE;
1920   }
1921
1922   /* Store the value of qmax used here */
1923   cv_mem->cv_qmax_alloc = qmax;
1924
1925   return(TRUE);
1926 }
1927
1928 /*  
1929  * CVFreeVectors
1930  *
1931  * This routine frees the CVODE vectors allocated in CVAllocVectors.
1932  */
1933
1934 static void CVFreeVectors(CVodeMem cv_mem)
1935 {
1936   int j, maxord;
1937   
1938   maxord = cv_mem->cv_qmax_alloc;
1939
1940   N_VDestroy(ewt);
1941   N_VDestroy(acor);
1942   N_VDestroy(tempv);
1943   N_VDestroy(ftemp);
1944   for(j=0; j <= maxord; j++) N_VDestroy(zn[j]);
1945
1946   lrw -= (maxord + 5)*lrw1;
1947   liw -= (maxord + 5)*liw1;
1948
1949   if (cv_mem->cv_VabstolMallocDone) {
1950     N_VDestroy(Vabstol);
1951     lrw -= lrw1;
1952     liw -= liw1;
1953   }
1954 }
1955
1956 /*  
1957  * CVInitialSetup
1958  *
1959  * This routine performs input consistency checks at the first step.
1960  * If needed, it also checks the linear solver module and calls the
1961  * linear solver initialization routine.
1962  */
1963
1964 static int CVInitialSetup(CVodeMem cv_mem)
1965 {
1966   int ier;
1967   int ewtsetOK;
1968
1969   /* Did the user provide efun? */
1970
1971   if (itol != CV_WF) {
1972     efun = CVEwtSet;
1973     e_data = (void *)cv_mem;
1974   } else {
1975     if (efun == NULL) {
1976       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVInitialSetup", MSGCV_NO_EFUN);
1977
1978           /* SUNDIALS EXTENSION */
1979           if (is_sundials_with_extension())
1980           {
1981                   return(CV_NO_EFUN);
1982           }
1983           else
1984           {
1985                   return(CV_ILL_INPUT);
1986           }
1987
1988     }
1989   }
1990
1991   ewtsetOK = efun(zn[0], ewt, e_data);
1992   if (ewtsetOK != 0) {
1993
1994     if (itol == CV_WF){ 
1995       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVInitialSetup", MSGCV_EWT_FAIL);
1996
1997           /* SUNDIALS EXTENSION */
1998           if (is_sundials_with_extension())
1999           {
2000                 return(CV_EWT_FAIL);
2001           }
2002           else
2003           {
2004                 return(CV_ILL_INPUT);
2005           }
2006     }else{ 
2007       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVInitialSetup", MSGCV_BAD_EWT);
2008
2009           /* SUNDIALS EXTENSION */
2010           if (is_sundials_with_extension())
2011           {
2012                   return(CV_BAD_EWT);
2013           }
2014           else
2015           {
2016                   return(CV_ILL_INPUT);
2017           }
2018
2019     }
2020   }
2021   
2022   /* Check if lsolve function exists (if needed)
2023      and call linit function (if it exists) */
2024
2025   if (iter == CV_NEWTON) {
2026     if (lsolve == NULL) {
2027       CVProcessError(cv_mem, CV_ILL_INPUT, "CVODE", "CVInitialSetup", MSGCV_LSOLVE_NULL);
2028
2029           /* SUNDIALS EXTENSION */
2030           if (is_sundials_with_extension())
2031           {
2032                 return(CV_LSOLVE_NULL);
2033           }
2034           else
2035           {
2036                 return(CV_ILL_INPUT);
2037           }
2038
2039     }
2040     if (linit != NULL) {
2041       ier = linit(cv_mem);
2042       if (ier != 0) {
2043         CVProcessError(cv_mem, CV_LINIT_FAIL, "CVODE", "CVInitialSetup", MSGCV_LINIT_FAIL);
2044         return(CV_LINIT_FAIL);
2045       }
2046     }
2047   }
2048
2049   return(CV_SUCCESS);
2050 }
2051
2052 /* 
2053  * -----------------------------------------------------------------
2054  * PRIVATE FUNCTIONS FOR CVODE
2055  * -----------------------------------------------------------------
2056  */
2057
2058 /*
2059  * CVHin
2060  *
2061  * This routine computes a tentative initial step size h0. 
2062  * If tout is too close to tn (= t0), then CVHin returns CV_TOO_CLOSE
2063  * and h remains uninitialized. 
2064  * If the RHS function fails unrecoverably, CVHin returns CV_RHSFUNC_FAIL.
2065  * If the RHS function fails recoverably too many times and recovery is
2066  * not possible, CVHin returns CV_REPTD_RHSFUNC_ERR.
2067  * Otherwise, CVHin sets h to the chosen value h0 and returns CV_SUCCESS.
2068  *
2069  * The algorithm used seeks to find h0 as a solution of
2070  *       (WRMS norm of (h0^2 ydd / 2)) = 1, 
2071  * where ydd = estimated second derivative of y.
2072  *
2073  * We start with an initial estimate equal to the geometric mean of the
2074  * lower and upper bounds on the step size.
2075  *
2076  * Loop up to MAX_ITERS times to find h0.
2077  * Stop if new and previous values differ by a factor < 2.
2078  * Stop if hnew/hg > 2 after one iteration, as this probably means
2079  * that the ydd value is bad because of cancellation error.        
2080  *  
2081  * For each new proposed hg, we allow MAX_ITERS attempts to
2082  * resolve a possible recoverable failure from f() by reducing
2083  * the proposed stepsize by a factor of 0.2. If a legal stepsize
2084  * still cannot be found, fall back on a previous value if possible,
2085  * or else return CV_REPTD_RHSFUNC_ERR.
2086  *
2087  * Finally, we apply a bias (0.5) and verify that h0 is within bounds.
2088  */
2089
2090 static int CVHin(CVodeMem cv_mem, realtype tout)
2091 {
2092   int retval = 0, sign = 0, count1 = 0, count2 = 0;
2093   realtype tdiff = 0., tdist = 0., tround = 0., hlb = 0., hub = 0.;
2094   realtype hg = 0., hgs = 0., hs = 0., hnew = 0., hrat = 0., h0 = 0., yddnrm = 0.;
2095   booleantype hgOK, hnewOK;
2096
2097   /* If tout is too close to tn, give up */
2098   
2099   if ((tdiff = tout-tn) == ZERO) return(CV_TOO_CLOSE);
2100   
2101   sign = (tdiff > ZERO) ? 1 : -1;
2102   tdist = ABS(tdiff);
2103   tround = uround * MAX(ABS(tn), ABS(tout));
2104
2105   if (tdist < TWO*tround) return(CV_TOO_CLOSE);
2106   
2107   /* 
2108      Set lower and upper bounds on h0, and take geometric mean 
2109      as first trial value.
2110      Exit with this value if the bounds cross each other.
2111   */
2112
2113   hlb = HLB_FACTOR * tround;
2114   hub = CVUpperBoundH0(cv_mem, tdist);
2115
2116   hg  = RSqrt(hlb*hub);
2117
2118   if (hub < hlb) {
2119     if (sign == -1) h = -hg;
2120     else            h =  hg;
2121     return(CV_SUCCESS);
2122   }
2123   
2124   /* Outer loop */
2125
2126   hnewOK = FALSE;
2127   hs = hg;         /* safeguard against 'uninitialized variable' warning */
2128
2129   for(count1 = 1; count1 <= MAX_ITERS; count1++) {
2130
2131     /* Attempts to estimate ydd */
2132
2133     hgOK = FALSE;
2134
2135     for (count2 = 1; count2 <= MAX_ITERS; count2++) {
2136       hgs = hg*sign;
2137       retval = CVYddNorm(cv_mem, hgs, &yddnrm);
2138       /* If f() failed unrecoverably, give up */
2139       if (retval < 0) return(CV_RHSFUNC_FAIL);
2140       /* If successful, we can use ydd */
2141       if (retval == CV_SUCCESS) {hgOK = TRUE; break;}
2142       /* f() failed recoverably; cut step size and test it again */
2143       hg *= POINT2;
2144     }
2145
2146     /* If f() failed recoverably MAX_ITERS times */
2147
2148     if (!hgOK) {
2149       /* Exit if this is the first or second pass. No recovery possible */
2150       if (count1 <= 2) return(CV_REPTD_RHSFUNC_ERR);
2151       /* We have a fall-back option. The value hs is a previous hnew which
2152          passed through f(). Use it and break */
2153       hnew = hs;
2154       break;
2155     }
2156
2157     /* The proposed step size is feasible. Save it. */
2158     hs = hg;
2159
2160     /* If the stopping criteria was met, or if this is the last pass, stop */
2161     if ( (hnewOK) || (count1 == MAX_ITERS))  {hnew = hg; break;}
2162
2163     /* Propose new step size */
2164     hnew = (yddnrm*hub*hub > TWO) ? RSqrt(TWO/yddnrm) : RSqrt(hg*hub);
2165     hrat = hnew/hg;
2166     
2167     /* Accept hnew if it does not differ from hg by more than a factor of 2 */
2168     if ((hrat > HALF) && (hrat < TWO)) {
2169       hnewOK = TRUE;
2170     }
2171
2172     /* After one pass, if ydd seems to be bad, use fall-back value. */
2173     if ((count1 > 1) && (hrat > TWO)) {
2174       hnew = hg;
2175       hnewOK = TRUE;
2176     }
2177
2178     /* Send this value back through f() */
2179     hg = hnew;
2180
2181   }
2182
2183   /* Apply bounds, bias factor, and attach sign */
2184
2185   h0 = H_BIAS*hnew;
2186   if (h0 < hlb) h0 = hlb;
2187   if (h0 > hub) h0 = hub;
2188   if (sign == -1) h0 = -h0;
2189   h = h0;
2190
2191   return(CV_SUCCESS);
2192 }
2193
2194 /*
2195  * CVHinFixed
2196  *
2197  * This routine computes the fixed step size h.
2198  * The objective is to approach hmax (= 1/hmax_inv) with h by trying to split the time interval (t-*told) into hmax-long parts.
2199  * - if t-*told is smaller than hmax, then set h = t-*told (one iteration)
2200  * - if it is divisible by hmax, then set h = hmax.
2201  * - if it is not, then "add an integration point" by setting h < hmax, just enough to fit the interval
2202  *
2203  * Dormand-Price being a fixed step size method, we know the maximum number of steps to take.
2204  * This procedure returns that number (minus 2 because nstloc starts at 0).
2205  */
2206
2207 static int CVHinFixed(CVodeMem cv_mem, realtype tout, realtype *tret)
2208 {
2209   long int n_points;
2210   realtype interval_size, hmax, test_div, floor_test;
2211
2212   hmax = 1./hmax_inv;
2213   interval_size = tout-*tret;
2214
2215   if (interval_size <= hmax) {  /* "Small" interval, h is the size of it */
2216     n_points = 2;
2217     h = interval_size;
2218   }
2219   else {
2220     test_div = interval_size/hmax;
2221     floor_test = FLOOR(test_div);
2222     if (test_div-floor_test <= TINY) {  /* t-*told divisible by hmax, cutting the interval into hmax-long parts */
2223       n_points = floor_test+1;
2224       h = interval_size/(n_points-1);
2225     }
2226     else {  /* Adding a point and h < hmax to fit the interval */
2227       n_points = floor_test+2;
2228       h = interval_size/(n_points-1);
2229     }
2230   }
2231   return(n_points-1);
2232 }
2233
2234 /*
2235  * CVUpperBoundH0
2236  *
2237  * This routine sets an upper bound on abs(h0) based on
2238  * tdist = tn - t0 and the values of y[i]/y'[i].
2239  */
2240
2241 static realtype CVUpperBoundH0(CVodeMem cv_mem, realtype tdist)
2242 {
2243   realtype hub_inv, hub;
2244   N_Vector temp1, temp2;
2245
2246   /* 
2247    * Bound based on |y0|/|y0'| -- allow at most an increase of
2248    * HUB_FACTOR in y0 (based on a forward Euler step). The weight 
2249    * factor is used as a safeguard against zero components in y0. 
2250    */
2251
2252   temp1 = tempv;
2253   temp2 = acor;
2254
2255   N_VAbs(zn[0], temp2);
2256   efun(zn[0], temp1, e_data);
2257   N_VInv(temp1, temp1);
2258   N_VLinearSum(HUB_FACTOR, temp2, ONE, temp1, temp1);
2259
2260   N_VAbs(zn[1], temp2);
2261
2262   N_VDiv(temp2, temp1, temp1);
2263   hub_inv = N_VMaxNorm(temp1);
2264
2265   /*
2266    * bound based on tdist -- allow at most a step of magnitude
2267    * HUB_FACTOR * tdist
2268    */
2269
2270   hub = HUB_FACTOR*tdist;
2271
2272   /* Use the smaler of the two */
2273
2274   if (hub*hub_inv > ONE) hub = ONE/hub_inv;
2275
2276   return(hub);
2277 }
2278
2279 /*
2280  * CVYddNorm
2281  *
2282  * This routine computes an estimate of the second derivative of y
2283  * using a difference quotient, and returns its WRMS norm.
2284  */
2285
2286 static int CVYddNorm(CVodeMem cv_mem, realtype hg, realtype *yddnrm)
2287 {
2288   int retval;
2289
2290   N_VLinearSum(hg, zn[1], ONE, zn[0], y);
2291   retval = f(tn+hg, y, tempv, f_data);
2292   nfe++;
2293   if (retval < 0) return(CV_RHSFUNC_FAIL);
2294   if (retval > 0) return(RHSFUNC_RECVR);
2295
2296   N_VLinearSum(ONE, tempv, -ONE, zn[1], tempv);
2297   N_VScale(ONE/hg, tempv, tempv);
2298
2299   *yddnrm = N_VWrmsNorm(tempv, ewt);
2300
2301   return(CV_SUCCESS);
2302 }
2303
2304 /* 
2305  * CVStep
2306  *
2307  * This routine performs one internal cvode step, from tn to tn + h.
2308  * It calls other routines to do all the work.
2309  *
2310  * The main operations done here are as follows:
2311  * - preliminary adjustments if a new step size was chosen;
2312  * - prediction of the Nordsieck history array zn at tn + h;
2313  * - setting of multistep method coefficients and test quantities;
2314  * - solution of the nonlinear system;
2315  * - testing the local error;
2316  * - updating zn and other state data if successful;
2317  * - resetting stepsize and order for the next step.
2318  * - if SLDET is on, check for stability, reduce order if necessary.
2319  * On a failure in the nonlinear system solution or error test, the
2320  * step may be reattempted, depending on the nature of the failure.
2321  */
2322
2323 static int CVStep(CVodeMem cv_mem)
2324 {
2325   realtype saved_t, dsm;
2326   int ncf, nef;
2327   int nflag, kflag, eflag;
2328   
2329   saved_t = tn;
2330   ncf = nef = 0;
2331   nflag = FIRST_CALL;
2332
2333   if ((nst > 0) && (hprime != h)) CVAdjustParams(cv_mem);
2334   
2335   /* Looping point for attempts to take a step */
2336   for(;;) {  
2337
2338     CVPredict(cv_mem);  
2339     CVSet(cv_mem);
2340
2341     nflag = CVNls(cv_mem, nflag);
2342     kflag = CVHandleNFlag(cv_mem, &nflag, saved_t, &ncf);
2343
2344     /* Go back in loop if we need to predict again (nflag=PREV_CONV_FAIL)*/
2345     if (kflag == PREDICT_AGAIN) continue;
2346
2347     /* Return if nonlinear solve failed and recovery not possible. */
2348     if (kflag != DO_ERROR_TEST) return(kflag);
2349
2350     /* Perform error test (nflag=CV_SUCCESS) */
2351     eflag = CVDoErrorTest(cv_mem, &nflag, saved_t, &nef, &dsm);
2352
2353     /* Go back in loop if we need to predict again (nflag=PREV_ERR_FAIL) */
2354     if (eflag == TRY_AGAIN)  continue;
2355
2356     /* Return if error test failed and recovery not possible. */
2357     if (eflag != CV_SUCCESS) return(eflag);
2358
2359     /* Error test passed (eflag=CV_SUCCESS), break from loop */
2360     break;
2361
2362   }
2363
2364   /* Nonlinear system solve and error test were both successful.
2365      Update data, and consider change of step and/or order.       */
2366
2367   CVCompleteStep(cv_mem); 
2368
2369   CVPrepareNextStep(cv_mem, dsm); 
2370
2371   /* If Stablilty Limit Detection is turned on, call stability limit
2372      detection routine for possible order reduction. */
2373
2374   if (sldeton) CVBDFStab(cv_mem);
2375
2376   etamax = (nst <= SMALL_NST) ? ETAMX2 : ETAMX3;
2377
2378   /*  Finally, we rescale the acor array to be the 
2379       estimated local error vector. */
2380
2381   N_VScale(ONE/tq[2], acor, acor);
2382   return(CV_SUCCESS);
2383       
2384 }
2385
2386 /*
2387  * CVStepDoPri
2388  *
2389  * This routine performs one internal cvode step using the Dormand-Price method, from tn to tn + h.
2390  * Proceed to computing the K[i] coefficients, build the final solution, increment tn and return CV_SUCCESS.
2391  *
2392  * In order to temporarily store the results, we use tempv and ftemp.
2393  * While ftempv is just going be an accumulator, tempv will represent the "correction solution".
2394  */
2395
2396 static int CVStepDoPri(CVodeMem cv_mem)
2397 {
2398   int retval;
2399
2400   /* Constant coefficients */
2401   realtype c2, c3, c4, c5, b1, b3, b4, b5, b6, d1, d3, d4, d5, d6, d7;
2402   realtype a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a73, a74, a75, a76;
2403   a21 =  0.2;                   c2 =  0.2;
2404   a31 =  0.075;                 c3 =  0.3;
2405   a32 =  0.225;                 c4 =  0.8;
2406   a41 =  0.9777777777777778;    c5 =  0.8888888888888888;
2407   a42 = -3.7333333333333333;    b1 =  0.0911458333333333;
2408   a43 =  3.5555555555555556;    b3 =  0.4492362982929021;
2409   a51 =  2.9525986892242036;    b4 =  0.6510416666666667;
2410   a52 = -11.595793324188386;    b5 = -0.3223761792452830;
2411   a53 =  9.8228928516994361;    b6 =  0.1309523809523810;
2412   a54 = -0.2908093278463649;    d1 =  0.0899131944444444;
2413   a61 =  2.8462752525252525;    d3 =  0.4534890685834082;
2414   a62 = -10.757575757575758;    d4 =  0.6140625;
2415   a63 =  8.9064227177434725;    d5 = -0.2715123820754717;
2416   a64 =  0.2784090909090909;    d6 =  0.0890476190476190;
2417   a65 = -0.2735313036020583;    d7 =  0.025;
2418   a71 =  0.0911458333333333;
2419   a73 =  0.4492362982929021;
2420   a74 =  0.6510416666666667;
2421   a75 = -0.3223761792452830;
2422   a76 =  0.1309523809523810;
2423
2424   /* K1, K2 */
2425   retval = f(tn, zn[0], zn[2], f_data);                    /* zn[2] = K1 = f(Tn, Yn), */
2426   N_VLinearSum_Serial (a21*h, zn[2], 1., zn[0], y);        /* y = K1*a21*h + Yn */
2427   retval = f(tn + c2*h, y, zn[3], f_data);                 /* zn[3] = K2 = f(Tn+c2*h, Yn + a21*h*K1), */
2428   N_VScale(b1, zn[2], y);                                  /* y = b1*K1 + 0*K2, */
2429   N_VScale(d1, zn[2], tempv);                              /* y = d1*K1 + 0*K2, */
2430
2431   /* K3 */
2432   N_VLinearSum_Serial (a31*h, zn[2], 1., zn[0], ftemp);    /* ftemp = Yn + K1*a31*h, */
2433   N_VLinearSum_Serial (a32*h, zn[3], 1., ftemp, ftemp);    /* ftemp += K2*a32*h, */
2434   retval = f(tn + c3*h, ftemp, zn[4], f_data);             /* zn[4] = K3 = f(Tn+c3*h, Yn + a31*h*K1 + a32*h*K2), */
2435   N_VLinearSum_Serial (1., y, b3, zn[4], y);               /* y = b1*K1 + 0*K2 + b3*K3, */
2436   N_VLinearSum_Serial (1., tempv, d3, zn[4], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3, */
2437
2438   /* K4 */
2439   N_VLinearSum_Serial (a41*h, zn[2], 1., zn[0], ftemp);    /* ftemp = Yn + K1*a41*h, */
2440   N_VLinearSum_Serial (a42*h, zn[3], 1., ftemp, ftemp);    /* ftemp += K2*a42*h, */
2441   N_VLinearSum_Serial (a43*h, zn[4], 1., ftemp, ftemp);    /* ftemp += K3*a43*h, */
2442   retval = f(tn + c4*h, ftemp, zn[5], f_data);             /* zn[5] = K4 = f(Tn+c4*h, Yn + a41*h*K1 + a42*h*K2 + a43*h*K3), */
2443   N_VLinearSum_Serial (1., y, b4, zn[5], y);               /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4, */
2444   N_VLinearSum_Serial (1., tempv, d4, zn[5], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4, */
2445
2446   /* K5 */
2447   N_VLinearSum_Serial (a51*h, zn[2], 1., zn[0], ftemp);    /* ftemp = Yn + K1*a51*h, */
2448   N_VLinearSum_Serial (a52*h, zn[3], 1., ftemp, ftemp);    /* ftemp += K2*a52*h, */
2449   N_VLinearSum_Serial (a53*h, zn[4], 1., ftemp, ftemp);    /* ftemp += K3*a53*h, */
2450   N_VLinearSum_Serial (a54*h, zn[5], 1., ftemp, ftemp);    /* ftemp += K4*a54*h, */
2451   retval = f(tn + c5*h, ftemp, zn[6], f_data);             /* zn[6] = K5 = f(Tn+c5*h, Yn + a51*h*K1 + a52*h*K2 + a53*h*K3 + a54*h*K4), */
2452   N_VLinearSum_Serial (1., y, b5, zn[6], y);               /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4 + b5*K5, */
2453   N_VLinearSum_Serial (1., tempv, d5, zn[6], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5, */
2454
2455   /* K6 */
2456   N_VLinearSum_Serial (a61*h, zn[2], 1., zn[0], ftemp);    /* ftemp = Yn + K1*a61*h, */
2457   N_VLinearSum_Serial (a62*h, zn[3], 1., ftemp, ftemp);    /* ftemp += K2*a62*h, */
2458   N_VLinearSum_Serial (a63*h, zn[4], 1., ftemp, ftemp);    /* ftemp += K3*a63*h, */
2459   N_VLinearSum_Serial (a64*h, zn[5], 1., ftemp, ftemp);    /* ftemp += K3*a64*h, */
2460   N_VLinearSum_Serial (a65*h, zn[6], 1., ftemp, ftemp);    /* ftemp += K3*a65*h, */
2461   retval = f(tn + h, ftemp, zn[7], f_data);                /* zn[7] = K6 = f(Tn+h, Yn + a61*h*K1 + a62*h*K2 + a63*h*K3 + a64*h*K4 + a65*h*K5), */
2462   N_VLinearSum_Serial (1., y, b6, zn[7], y);               /* y = b1*K1 + 0*K2 + b3*K3 + b4*K4 + b5*K5 + b6*K6, */
2463   N_VLinearSum_Serial (1., tempv, d6, zn[7], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5 + d6*K6, */
2464
2465   /* K7 */
2466   N_VLinearSum_Serial (a71*h, zn[2], 1., zn[0], ftemp);    /* ftemp = Yn + K1*a71*h, */
2467   N_VLinearSum_Serial (a73*h, zn[4], 1., ftemp, ftemp);    /* ftemp += K3*a73*h, */
2468   N_VLinearSum_Serial (a74*h, zn[5], 1., ftemp, ftemp);    /* ftemp += K3*a74*h, */
2469   N_VLinearSum_Serial (a75*h, zn[6], 1., ftemp, ftemp);    /* ftemp += K3*a75*h, */
2470   N_VLinearSum_Serial (a76*h, zn[7], 1., ftemp, ftemp);    /* ftemp += K3*a76*h, */
2471   retval = f(tn + h, ftemp, zn[8], f_data);                /* zn[8] = K7 = f(Tn+h, Yn + a71*h*K1 + 0*h*K2 + a73*h*K3 + a74*h*K4 + a75*h*K5 + a76*h*K6), */
2472   N_VLinearSum_Serial (1., tempv, d7, zn[8], tempv);       /* tempv = d1*K1 + 0*K2 + d3*K3 + d4*K4 + d5*K5 + d6*K6 + d7*K7, */
2473
2474   /* Yn+1 */
2475   N_VLinearSum_Serial(1., zn[0], h, y, zn[0]);             /* zn[0] = Yn+1 = Yn + h*y */
2476
2477   /* Check for errors in the evaluations of f thanks to retval */
2478   if (retval < 0) {
2479     CVProcessError(cv_mem, CV_RHSFUNC_FAIL, "Dormand-Price", "CVStepDoPri", MSGCV_RHSFUNC_FAILED, tn);
2480     return(CV_RHSFUNC_FAIL);
2481   }
2482   if (retval > 0) {
2483     CVProcessError(cv_mem, CV_FIRST_RHSFUNC_ERR, "Dormand-Price", "CVStepDoPri", MSGCV_RHSFUNC_FIRST);
2484     return(CV_FIRST_RHSFUNC_ERR);
2485   }
2486
2487   /* Increment tn => take a step. Increment solver calls as well */
2488   tn += h;
2489   nst++;
2490
2491   /* Update the Nordsieck history array */
2492   retval = f(tn, zn[0], zn[1], f_data);    /* zn[1] = y'(tn) */
2493
2494   return(CV_SUCCESS);
2495 }
2496
2497 /*
2498  * CVAdjustParams
2499  *
2500  * This routine is called when a change in step size was decided upon,
2501  * and it handles the required adjustments to the history array zn.
2502  * If there is to be a change in order, we call CVAdjustOrder and reset
2503  * q, L = q+1, and qwait.  Then in any case, we call CVRescale, which
2504  * resets h and rescales the Nordsieck array.
2505  */
2506
2507 static void CVAdjustParams(CVodeMem cv_mem)
2508 {
2509   if (qprime != q) {
2510     CVAdjustOrder(cv_mem, qprime-q);
2511     q = qprime;
2512     L = q+1;
2513     qwait = L;
2514   }
2515   CVRescale(cv_mem);
2516 }
2517
2518 /*
2519  * CVAdjustOrder
2520  *
2521  * This routine is a high level routine which handles an order
2522  * change by an amount deltaq (= +1 or -1). If a decrease in order
2523  * is requested and q==2, then the routine returns immediately.
2524  * Otherwise CVAdjustAdams or CVAdjustBDF is called to handle the
2525  * order change (depending on the value of lmm).
2526  */
2527
2528 static void CVAdjustOrder(CVodeMem cv_mem, int deltaq)
2529 {
2530   if ((q==2) && (deltaq != 1)) return;
2531   
2532   switch(lmm){
2533   case CV_ADAMS: 
2534     CVAdjustAdams(cv_mem, deltaq);
2535     break;
2536   case CV_BDF:   
2537     CVAdjustBDF(cv_mem, deltaq);
2538     break;
2539   }
2540 }
2541
2542 /*
2543  * CVAdjustAdams
2544  *
2545  * This routine adjusts the history array on a change of order q by
2546  * deltaq, in the case that lmm == CV_ADAMS.
2547  */
2548
2549 static void CVAdjustAdams(CVodeMem cv_mem, int deltaq)
2550 {
2551   int i, j;
2552   realtype xi, hsum;
2553
2554   /* On an order increase, set new column of zn to zero and return */
2555   
2556   if (deltaq==1) {
2557     N_VConst(ZERO, zn[L]);
2558     return;
2559   }
2560
2561   /*
2562    * On an order decrease, each zn[j] is adjusted by a multiple of zn[q].
2563    * The coeffs. in the adjustment are the coeffs. of the polynomial:
2564    *        x
2565    * q * INT { u * ( u + xi_1 ) * ... * ( u + xi_{q-2} ) } du 
2566    *        0
2567    * where xi_j = [t_n - t_(n-j)]/h => xi_0 = 0
2568    */
2569
2570   for (i=0; i <= qmax; i++) l[i] = ZERO;
2571   l[1] = ONE;
2572   hsum = ZERO;
2573   for (j=1; j <= q-2; j++) {
2574     hsum += tau[j];
2575     xi = hsum / hscale;
2576     for (i=j+1; i >= 1; i--) l[i] = l[i]*xi + l[i-1];
2577   }
2578   
2579   for (j=1; j <= q-2; j++) l[j+1] = q * (l[j] / (j+1));
2580   
2581   for (j=2; j < q; j++)
2582     N_VLinearSum(-l[j], zn[q], ONE, zn[j], zn[j]);
2583 }
2584
2585 /*
2586  * CVAdjustBDF
2587  *
2588  * This is a high level routine which handles adjustments to the
2589  * history array on a change of order by deltaq in the case that 
2590  * lmm == CV_BDF.  CVAdjustBDF calls CVIncreaseBDF if deltaq = +1 and 
2591  * CVDecreaseBDF if deltaq = -1 to do the actual work.
2592  */
2593
2594 static void CVAdjustBDF(CVodeMem cv_mem, int deltaq)
2595 {
2596   switch(deltaq) {
2597   case 1 : 
2598     CVIncreaseBDF(cv_mem);
2599     return;
2600   case -1: 
2601     CVDecreaseBDF(cv_mem);
2602     return;
2603   }
2604 }
2605
2606 /*
2607  * CVIncreaseBDF
2608  *
2609  * This routine adjusts the history array on an increase in the 
2610  * order q in the case that lmm == CV_BDF.  
2611  * A new column zn[q+1] is set equal to a multiple of the saved 
2612  * vector (= acor) in zn[indx_acor].  Then each zn[j] is adjusted by
2613  * a multiple of zn[q+1].  The coefficients in the adjustment are the 
2614  * coefficients of the polynomial x*x*(x+xi_1)*...*(x+xi_j),
2615  * where xi_j = [t_n - t_(n-j)]/h.
2616  */
2617
2618 static void CVIncreaseBDF(CVodeMem cv_mem)
2619 {
2620   realtype alpha0, alpha1, prod, xi, xiold, hsum, A1;
2621   int i, j;
2622   
2623   for (i=0; i <= qmax; i++) l[i] = ZERO;
2624   l[2] = alpha1 = prod = xiold = ONE;
2625   alpha0 = -ONE;
2626   hsum = hscale;
2627   if (q > 1) {
2628     for (j=1; j < q; j++) {
2629       hsum += tau[j+1];
2630       xi = hsum / hscale;
2631       prod *= xi;
2632       alpha0 -= ONE / (j+1);
2633       alpha1 += ONE / xi;
2634       for (i=j+2; i >= 2; i--) l[i] = l[i]*xiold + l[i-1];
2635       xiold = xi;
2636     }
2637   }
2638   A1 = (-alpha0 - alpha1) / prod;
2639   N_VScale(A1, zn[indx_acor], zn[L]);
2640   for (j=2; j <= q; j++) {
2641     N_VLinearSum(l[j], zn[L], ONE, zn[j], zn[j]);
2642   }  
2643 }
2644
2645 /*
2646  * CVDecreaseBDF
2647  *
2648  * This routine adjusts the history array on a decrease in the 
2649  * order q in the case that lmm == CV_BDF.  
2650  * Each zn[j] is adjusted by a multiple of zn[q].  The coefficients
2651  * in the adjustment are the coefficients of the polynomial
2652  *   x*x*(x+xi_1)*...*(x+xi_j), where xi_j = [t_n - t_(n-j)]/h.
2653  */
2654
2655 static void CVDecreaseBDF(CVodeMem cv_mem)
2656 {
2657   realtype hsum, xi;
2658   int i, j;
2659   
2660   for (i=0; i <= qmax; i++) l[i] = ZERO;
2661   l[2] = ONE;
2662   hsum = ZERO;
2663   for(j=1; j <= q-2; j++) {
2664     hsum += tau[j];
2665     xi = hsum /hscale;
2666     for (i=j+2; i >= 2; i--) l[i] = l[i]*xi + l[i-1];
2667   }
2668   
2669   for(j=2; j < q; j++)
2670     N_VLinearSum(-l[j], zn[q], ONE, zn[j], zn[j]);
2671 }
2672
2673 /*
2674  * CVRescale
2675  *
2676  * This routine rescales the Nordsieck array by multiplying the
2677  * jth column zn[j] by eta^j, j = 1, ..., q.  Then the value of
2678  * h is rescaled by eta, and hscale is reset to h.
2679  */
2680
2681 static void CVRescale(CVodeMem cv_mem)
2682 {
2683   int j;
2684   realtype factor;
2685   
2686   factor = eta;
2687   for (j=1; j <= q; j++) {
2688     N_VScale(factor, zn[j], zn[j]);
2689     factor *= eta;
2690   }
2691   h = hscale * eta;
2692   next_h = h;
2693   hscale = h;
2694   nscon = 0;
2695 }
2696
2697 /*
2698  * CVPredict
2699  *
2700  * This routine advances tn by the tentative step size h, and computes
2701  * the predicted array z_n(0), which is overwritten on zn.  The
2702  * prediction of zn is done by repeated additions.
2703  * In TSTOP mode, it is possible for tn + h to be past tstop by roundoff,
2704  * and in that case, we reset tn (after incrementing by h) to tstop.
2705  */
2706
2707 static void CVPredict(CVodeMem cv_mem)
2708 {
2709   int j, k;
2710   
2711   tn += h;
2712   if (istop) {
2713     if ((tn - tstop)*h > ZERO) tn = tstop;
2714   }
2715   for (k = 1; k <= q; k++)
2716     for (j = q; j >= k; j--) 
2717       N_VLinearSum(ONE, zn[j-1], ONE, zn[j], zn[j-1]); 
2718 }
2719
2720 /*
2721  * CVSet
2722  *
2723  * This routine is a high level routine which calls CVSetAdams or
2724  * CVSetBDF to set the polynomial l, the test quantity array tq, 
2725  * and the related variables  rl1, gamma, and gamrat.
2726  */
2727
2728 static void CVSet(CVodeMem cv_mem)
2729 {
2730   switch(lmm) {
2731   case CV_ADAMS: 
2732     CVSetAdams(cv_mem);
2733     break;
2734   case CV_BDF:
2735     CVSetBDF(cv_mem);
2736     break;
2737   }
2738   rl1 = ONE / l[1];
2739   gamma = h * rl1;
2740   if (nst == 0) gammap = gamma;
2741   gamrat = (nst > 0) ? gamma / gammap : ONE;  /* protect x / x != 1.0 */
2742 }
2743
2744 /*
2745  * CVSetAdams
2746  *
2747  * This routine handles the computation of l and tq for the
2748  * case lmm == CV_ADAMS.
2749  *
2750  * The components of the array l are the coefficients of a
2751  * polynomial Lambda(x) = l_0 + l_1 x + ... + l_q x^q, given by
2752  *                          q-1
2753  * (d/dx) Lambda(x) = c * PRODUCT (1 + x / xi_i) , where
2754  *                          i=1
2755  *  Lambda(-1) = 0, Lambda(0) = 1, and c is a normalization factor.
2756  * Here xi_i = [t_n - t_(n-i)] / h.
2757  *
2758  * The array tq is set to test quantities used in the convergence
2759  * test, the error test, and the selection of h at a new order.
2760  */
2761
2762 static void CVSetAdams(CVodeMem cv_mem)
2763 {
2764   realtype m[L_MAX], M[3], hsum;
2765   
2766   if (q == 1) {
2767     l[0] = l[1] = tq[1] = tq[5] = ONE;
2768     tq[2] = TWO;
2769     tq[3] = TWELVE;
2770     tq[4] = nlscoef * tq[2];       /* = 0.1 * tq[2] */
2771     return;
2772   }
2773   
2774   hsum = CVAdamsStart(cv_mem, m);
2775   
2776   M[0] = CVAltSum(q-1, m, 1);
2777   M[1] = CVAltSum(q-1, m, 2);
2778   
2779   CVAdamsFinish(cv_mem, m, M, hsum);
2780 }
2781
2782 /*
2783  * CVAdamsStart
2784  *
2785  * This routine generates in m[] the coefficients of the product
2786  * polynomial needed for the Adams l and tq coefficients for q > 1.
2787  */
2788
2789 static realtype CVAdamsStart(CVodeMem cv_mem, realtype m[])
2790 {
2791   realtype hsum, xi_inv, sum;
2792   int i, j;
2793   
2794   hsum = h;
2795   m[0] = ONE;
2796   for (i=1; i <= q; i++) m[i] = ZERO;
2797   for (j=1; j < q; j++) {
2798     if ((j==q-1) && (qwait == 1)) {
2799       sum = CVAltSum(q-2, m, 2);
2800       tq[1] = m[q-2] / (q * sum);
2801     }
2802     xi_inv = h / hsum;
2803     for (i=j; i >= 1; i--) m[i] += m[i-1] * xi_inv;
2804     hsum += tau[j];
2805     /* The m[i] are coefficients of product(1 to j) (1 + x/xi_i) */
2806   }
2807   return(hsum);
2808 }
2809
2810 /*
2811  * CVAdamsFinish
2812  *
2813  * This routine completes the calculation of the Adams l and tq.
2814  */
2815
2816 static void CVAdamsFinish(CVodeMem cv_mem, realtype m[], realtype M[], realtype hsum)
2817 {
2818   int i;
2819   realtype M0_inv, xi, xi_inv;
2820   
2821   M0_inv = ONE / M[0];
2822   
2823   l[0] = ONE;
2824   for (i=1; i <= q; i++) l[i] = M0_inv * (m[i-1] / i);
2825   xi = hsum / h;
2826   xi_inv = ONE / xi;
2827   
2828   tq[2] = xi * M[0] / M[1];
2829   tq[5] = xi / l[q];
2830
2831   if (qwait == 1) {
2832     for (i=q; i >= 1; i--) m[i] += m[i-1] * xi_inv;
2833     M[2] = CVAltSum(q, m, 2);
2834     tq[3] = L * M[0] / M[2];
2835   }
2836
2837   tq[4] = nlscoef * tq[2];
2838 }
2839
2840 /*  
2841  * CVAltSum
2842  *
2843  * CVAltSum returns the value of the alternating sum
2844  *   sum (i= 0 ... iend) [ (-1)^i * (a[i] / (i + k)) ].
2845  * If iend < 0 then CVAltSum returns 0.
2846  * This operation is needed to compute the integral, from -1 to 0,
2847  * of a polynomial x^(k-1) M(x) given the coefficients of M(x).
2848  */
2849
2850 static realtype CVAltSum(int iend, realtype a[], int k)
2851 {
2852   int i, sign;
2853   realtype sum;
2854   
2855   if (iend < 0) return(ZERO);
2856   
2857   sum = ZERO;
2858   sign = 1;
2859   for (i=0; i <= iend; i++) {
2860     sum += sign * (a[i] / (i+k));
2861     sign = -sign;
2862   }
2863   return(sum);
2864 }
2865
2866 /*
2867  * CVSetBDF
2868  *
2869  * This routine computes the coefficients l and tq in the case
2870  * lmm == CV_BDF.  CVSetBDF calls CVSetTqBDF to set the test
2871  * quantity array tq. 
2872  * 
2873  * The components of the array l are the coefficients of a
2874  * polynomial Lambda(x) = l_0 + l_1 x + ... + l_q x^q, given by
2875  *                                 q-1
2876  * Lambda(x) = (1 + x / xi*_q) * PRODUCT (1 + x / xi_i) , where
2877  *                                 i=1
2878  *  xi_i = [t_n - t_(n-i)] / h.
2879  *
2880  * The array tq is set to test quantities used in the convergence
2881  * test, the error test, and the selection of h at a new order.
2882  */
2883
2884 static void CVSetBDF(CVodeMem cv_mem)
2885 {
2886   realtype alpha0, alpha0_hat, xi_inv, xistar_inv, hsum;
2887   int i,j;
2888   
2889   l[0] = l[1] = xi_inv = xistar_inv = ONE;
2890   for (i=2; i <= q; i++) l[i] = ZERO;
2891   alpha0 = alpha0_hat = -ONE;
2892   hsum = h;
2893   if (q > 1) {
2894     for (j=2; j < q; j++) {
2895       hsum += tau[j-1];
2896       xi_inv = h / hsum;
2897       alpha0 -= ONE / j;
2898       for(i=j; i >= 1; i--) l[i] += l[i-1]*xi_inv;
2899       /* The l[i] are coefficients of product(1 to j) (1 + x/xi_i) */
2900     }
2901     
2902     /* j = q */
2903     alpha0 -= ONE / q;
2904     xistar_inv = -l[1] - alpha0;
2905     hsum += tau[q-1];
2906     xi_inv = h / hsum;
2907     alpha0_hat = -l[1] - xi_inv;
2908     for (i=q; i >= 1; i--) l[i] += l[i-1]*xistar_inv;
2909   }
2910
2911   CVSetTqBDF(cv_mem, hsum, alpha0, alpha0_hat, xi_inv, xistar_inv);
2912 }
2913
2914 /*
2915  * CVSetTqBDF
2916  *
2917  * This routine sets the test quantity array tq in the case
2918  * lmm == CV_BDF.
2919  */
2920
2921 static void CVSetTqBDF(CVodeMem cv_mem, realtype hsum, realtype alpha0,
2922                        realtype alpha0_hat, realtype xi_inv, realtype xistar_inv)
2923 {
2924   realtype A1, A2, A3, A4, A5, A6;
2925   realtype C, CPrime, CPrimePrime;
2926   
2927   A1 = ONE - alpha0_hat + alpha0;
2928   A2 = ONE + q * A1;
2929   tq[2] = ABS(alpha0 * (A2 / A1));
2930   tq[5] = ABS((A2) / (l[q] * xi_inv/xistar_inv));
2931   if (qwait == 1) {
2932     C = xistar_inv / l[q];
2933     A3 = alpha0 + ONE / q;
2934     A4 = alpha0_hat + xi_inv;
2935     CPrime = A3 / (ONE - A4 + A3);
2936     tq[1] = ABS(CPrime / C);
2937     hsum += tau[q];
2938     xi_inv = h / hsum;
2939     A5 = alpha0 - (ONE / (q+1));
2940     A6 = alpha0_hat - xi_inv;
2941     CPrimePrime = A2 / (ONE - A6 + A5);
2942     tq[3] = ABS(CPrimePrime * xi_inv * (q+2) * A5);
2943   }
2944   tq[4] = nlscoef * tq[2];
2945 }
2946
2947 /*
2948  * CVNls
2949  *
2950  * This routine attempts to solve the nonlinear system associated
2951  * with a single implicit step of the linear multistep method.
2952  * Depending on iter, it calls CVNlsFunctional or CVNlsNewton
2953  * to do the work.
2954  */
2955
2956 static int CVNls(CVodeMem cv_mem, int nflag)
2957 {
2958   int flag = CV_SUCCESS;
2959
2960   switch(iter) {
2961   case CV_FUNCTIONAL: 
2962     flag = CVNlsFunctional(cv_mem);
2963     break;
2964   case CV_NEWTON:
2965     flag = CVNlsNewton(cv_mem, nflag);
2966     break;
2967   }
2968
2969   return(flag);
2970 }
2971
2972 /*
2973  * CVNlsFunctional
2974  *
2975  * This routine attempts to solve the nonlinear system using 
2976  * functional iteration (no matrices involved).
2977  *
2978  * Possible return values are:
2979  *
2980  *   CV_SUCCESS      --->  continue with error test
2981  *
2982  *   CV_RHSFUNC_FAIL --->  halt the integration
2983  *
2984  *   CONV_FAIL       -+
2985  *   RHSFUNC_RECVR   -+->  predict again or stop if too many
2986  *
2987  */
2988
2989 static int CVNlsFunctional(CVodeMem cv_mem)
2990 {
2991   int retval, m;
2992   realtype del, delp, dcon;
2993
2994   /* Initialize counter and evaluate f at predicted y */
2995   
2996   crate = ONE;
2997   m = 0;
2998
2999   retval = f(tn, zn[0], tempv, f_data);
3000   nfe++;
3001   if (retval < 0) return(CV_RHSFUNC_FAIL);
3002   if (retval > 0) return(RHSFUNC_RECVR);
3003
3004   N_VConst(ZERO, acor);
3005
3006   /* Initialize delp to avoid compiler warning message */
3007   del = delp = ZERO;
3008
3009   /* Loop until convergence; accumulate corrections in acor */
3010
3011   for(;;) {
3012
3013     nni++;
3014
3015     /* Correct y directly from the last f value */
3016     N_VLinearSum(h, tempv, -ONE, zn[1], tempv);
3017     N_VScale(rl1, tempv, tempv);
3018     N_VLinearSum(ONE, zn[0], ONE, tempv, y);
3019     /* Get WRMS norm of current correction to use in convergence test */
3020     N_VLinearSum(ONE, tempv, -ONE, acor, acor);
3021     del = N_VWrmsNorm(acor, ewt);
3022     N_VScale(ONE, tempv, acor);
3023     
3024     /* Test for convergence.  If m > 0, an estimate of the convergence
3025        rate constant is stored in crate, and used in the test.        */
3026     if (m > 0) crate = MAX(CRDOWN * crate, del / delp);
3027     dcon = del * MIN(ONE, crate) / tq[4];
3028     if (dcon <= ONE) {
3029       acnrm = (m == 0) ? del : N_VWrmsNorm(acor, ewt);
3030       return(CV_SUCCESS);  /* Convergence achieved */
3031     }
3032
3033     /* Stop at maxcor iterations or if iter. seems to be diverging */
3034     m++;
3035     if ((m==maxcor) || ((m >= 2) && (del > RDIV * delp))) return(CONV_FAIL);
3036
3037     /* Save norm of correction, evaluate f, and loop again */
3038     delp = del;
3039
3040     retval = f(tn, y, tempv, f_data);
3041     nfe++;
3042     if (retval < 0) return(CV_RHSFUNC_FAIL);
3043     if (retval > 0) return(RHSFUNC_RECVR);
3044
3045   }
3046 }
3047
3048 /*
3049  * CVNlsNewton
3050  *
3051  * This routine handles the Newton iteration. It calls lsetup if 
3052  * indicated, calls CVNewtonIteration to perform the iteration, and 
3053  * retries a failed attempt at Newton iteration if that is indicated.
3054  *
3055  * Possible return values:
3056  *
3057  *   CV_SUCCESS       ---> continue with error test
3058  *
3059  *   CV_RHSFUNC_FAIL  -+  
3060  *   CV_LSETUP_FAIL    |-> halt the integration 
3061  *   CV_LSOLVE_FAIL   -+
3062  *
3063  *   CONV_FAIL        -+
3064  *   RHSFUNC_RECVR    -+-> predict again or stop if too many
3065  *
3066  */
3067
3068 static int CVNlsNewton(CVodeMem cv_mem, int nflag)
3069 {
3070   N_Vector vtemp1, vtemp2, vtemp3;
3071   int convfail, retval, ier;
3072   booleantype callSetup;
3073   
3074   vtemp1 = acor;  /* rename acor as vtemp1 for readability  */
3075   vtemp2 = y;     /* rename y as vtemp2 for readability     */
3076   vtemp3 = tempv; /* rename tempv as vtemp3 for readability */
3077   
3078   /* Set flag convfail, input to lsetup for its evaluation decision */
3079   convfail = ((nflag == FIRST_CALL) || (nflag == PREV_ERR_FAIL)) ?
3080     CV_NO_FAILURES : CV_FAIL_OTHER;
3081
3082   /* Decide whether or not to call setup routine (if one exists) */
3083   if (setupNonNull) {      
3084     callSetup = (nflag == PREV_CONV_FAIL) || (nflag == PREV_ERR_FAIL) ||
3085       (nst == 0) || (nst >= nstlp + MSBP) || (ABS(gamrat-ONE) > DGMAX);
3086   } else {  
3087     crate = ONE;
3088     callSetup = FALSE;
3089   }
3090   
3091   /* Looping point for the solution of the nonlinear system.
3092      Evaluate f at the predicted y, call lsetup if indicated, and
3093      call CVNewtonIteration for the Newton iteration itself.      */
3094   
3095   for(;;) {
3096
3097     retval = f(tn, zn[0], ftemp, f_data);
3098     nfe++; 
3099     if (retval < 0) return(CV_RHSFUNC_FAIL);
3100     if (retval > 0) return(RHSFUNC_RECVR);
3101
3102     if (callSetup) {
3103       ier = lsetup(cv_mem, convfail, zn[0], ftemp, &jcur, 
3104                    vtemp1, vtemp2, vtemp3);
3105       nsetups++;
3106       callSetup = FALSE;
3107       gamrat = crate = ONE; 
3108       gammap = gamma;
3109       nstlp = nst;
3110       /* Return if lsetup failed */
3111       if (ier < 0) return(CV_LSETUP_FAIL);
3112       if (ier > 0) return(CONV_FAIL);
3113     }
3114
3115     /* Set acor to zero and load prediction into y vector */
3116     N_VConst(ZERO, acor);
3117     N_VScale(ONE, zn[0], y);
3118
3119     /* Do the Newton iteration */
3120     ier = CVNewtonIteration(cv_mem);
3121
3122     /* If there is a convergence failure and the Jacobian-related 
3123        data appears not to be current, loop again with a call to lsetup
3124        in which convfail=CV_FAIL_BAD_J.  Otherwise return.                 */
3125     if (ier != TRY_AGAIN) return(ier);
3126     
3127     callSetup = TRUE;
3128     convfail = CV_FAIL_BAD_J;
3129   }
3130 }
3131
3132 /*
3133  * CVNewtonIteration
3134  *
3135  * This routine performs the Newton iteration. If the iteration succeeds,
3136  * it returns the value CV_SUCCESS. If not, it may signal the CVNlsNewton 
3137  * routine to call lsetup again and reattempt the iteration, by
3138  * returning the value TRY_AGAIN. (In this case, CVNlsNewton must set 
3139  * convfail to CV_FAIL_BAD_J before calling setup again). 
3140  * Otherwise, this routine returns one of the appropriate values 
3141  * CV_LSOLVE_FAIL, CV_RHSFUNC_FAIL, CONV_FAIL, or RHSFUNC_RECVR back 
3142  * to CVNlsNewton.
3143  */
3144
3145 static int CVNewtonIteration(CVodeMem cv_mem)
3146 {
3147   int m, retval;
3148   realtype del, delp, dcon;
3149   N_Vector b;
3150
3151   mnewt = m = 0;
3152
3153   /* Initialize delp to avoid compiler warning message */
3154   del = delp = ZERO;
3155
3156   /* Looping point for Newton iteration */
3157   for(;;) {
3158
3159     /* Evaluate the residual of the nonlinear system*/
3160     N_VLinearSum(rl1, zn[1], ONE, acor, tempv);
3161     N_VLinearSum(gamma, ftemp, -ONE, tempv, tempv);
3162
3163     /* Call the lsolve function */
3164     b = tempv;
3165     retval = lsolve(cv_mem, b, ewt, y, ftemp); 
3166     nni++;
3167     
3168     if (retval < 0) return(CV_LSOLVE_FAIL);
3169     
3170     /* If lsolve had a recoverable failure and Jacobian data is
3171        not current, signal to try the solution again            */
3172     if (retval > 0) { 
3173       if ((!jcur) && (setupNonNull)) return(TRY_AGAIN);
3174       else                           return(CONV_FAIL);
3175     }
3176
3177     /* Get WRMS norm of correction; add correction to acor and y */
3178     del = N_VWrmsNorm(b, ewt);
3179     N_VLinearSum(ONE, acor, ONE, b, acor);
3180     N_VLinearSum(ONE, zn[0], ONE, acor, y);
3181     
3182     /* Test for convergence.  If m > 0, an estimate of the convergence
3183        rate constant is stored in crate, and used in the test.        */
3184     if (m > 0) {
3185       crate = MAX(CRDOWN * crate, del/delp);
3186     }
3187     dcon = del * MIN(ONE, crate) / tq[4];
3188     
3189     if (dcon <= ONE) {
3190       acnrm = (m==0) ? del : N_VWrmsNorm(acor, ewt);
3191       jcur = FALSE;
3192       return(CV_SUCCESS); /* Nonlinear system was solved successfully */
3193     }
3194     
3195     mnewt = ++m;
3196     
3197     /* Stop at maxcor iterations or if iter. seems to be diverging.
3198        If still not converged and Jacobian data is not current, 
3199        signal to try the solution again                            */
3200     if ((m == maxcor) || ((m >= 2) && (del > RDIV*delp))) {
3201       if ((!jcur) && (setupNonNull)) return(TRY_AGAIN);
3202       else                           return(CONV_FAIL);
3203     }
3204     
3205     /* Save norm of correction, evaluate f, and loop again */
3206     delp = del;
3207     retval = f(tn, y, ftemp, f_data);
3208     nfe++;
3209     if (retval < 0) return(CV_RHSFUNC_FAIL);
3210     if (retval > 0) {
3211       if ((!jcur) && (setupNonNull)) return(TRY_AGAIN);
3212       else                           return(RHSFUNC_RECVR);
3213     }
3214
3215   } /* end loop */
3216 }
3217
3218 /*
3219  * CVHandleFlag
3220  *
3221  * This routine takes action on the return value nflag = *nflagPtr
3222  * returned by CVNls, as follows:
3223  *
3224  * If CVNls succeeded in solving the nonlinear system, then
3225  * CVHandleNFlag returns the constant DO_ERROR_TEST, which tells CVStep
3226  * to perform the error test.
3227  *
3228  * If the nonlinear system was not solved successfully, then ncfn and
3229  * ncf = *ncfPtr are incremented and Nordsieck array zn is restored.
3230  *
3231  * If the solution of the nonlinear system failed due to an
3232  * unrecoverable failure by setup, we return the value CV_LSETUP_FAIL.
3233  * 
3234  * If it failed due to an unrecoverable failure in solve, then we return
3235  * the value CV_LSOLVE_FAIL.
3236  *
3237  * If it failed due to an unrecoverable failure in rhs, then we return
3238  * the value CV_RHSFUNC_FAIL.
3239  *
3240  * Otherwise, a recoverable failure occurred when solving the 
3241  * nonlinear system (CVNls returned nflag == CONV_FAIL or RHSFUNC_RECVR). 
3242  * In this case, if ncf is now equal to maxncf or |h| = hmin, 
3243  * we return the value CV_CONV_FAILURE (if nflag=CONV_FAIL) or
3244  * CV_REPTD_RHSFUNC_ERR (if nflag=RHSFUNC_RECVR).
3245  * If not, we set *nflagPtr = PREV_CONV_FAIL and return the value
3246  * PREDICT_AGAIN, telling CVStep to reattempt the step.
3247  *
3248  */
3249
3250 static int CVHandleNFlag(CVodeMem cv_mem, int *nflagPtr, realtype saved_t,
3251                          int *ncfPtr)
3252 {
3253   int nflag;
3254   
3255   nflag = *nflagPtr;
3256   
3257   if (nflag == CV_SUCCESS) return(DO_ERROR_TEST);
3258
3259   /* The nonlinear soln. failed; increment ncfn and restore zn */
3260   ncfn++;
3261   CVRestore(cv_mem, saved_t);
3262   
3263   /* Return if lsetup, lsolve, or rhs failed unrecoverably */
3264   if (nflag == CV_LSETUP_FAIL)  return(CV_LSETUP_FAIL);
3265   if (nflag == CV_LSOLVE_FAIL)  return(CV_LSOLVE_FAIL);
3266   if (nflag == CV_RHSFUNC_FAIL) return(CV_RHSFUNC_FAIL);
3267   
3268   /* At this point, nflag = CONV_FAIL or RHSFUNC_RECVR; increment ncf */
3269   
3270   (*ncfPtr)++;
3271   etamax = ONE;
3272
3273   /* If we had maxncf failures or |h| = hmin, 
3274      return CV_CONV_FAILURE or CV_REPTD_RHSFUNC_ERR. */
3275
3276   if ((ABS(h) <= hmin*ONEPSM) || (*ncfPtr == maxncf)) {
3277     if (nflag == CONV_FAIL)     return(CV_CONV_FAILURE);
3278     if (nflag == RHSFUNC_RECVR) return(CV_REPTD_RHSFUNC_ERR);    
3279   }
3280
3281   /* Reduce step size; return to reattempt the step */
3282
3283   eta = MAX(ETACF, hmin / ABS(h));
3284   *nflagPtr = PREV_CONV_FAIL;
3285   CVRescale(cv_mem);
3286
3287   return(PREDICT_AGAIN);
3288 }
3289
3290 /*
3291  * CVRestore
3292  *
3293  * This routine restores the value of tn to saved_t and undoes the
3294  * prediction.  After execution of CVRestore, the Nordsieck array zn has
3295  * the same values as before the call to CVPredict.
3296  */
3297
3298 static void CVRestore(CVodeMem cv_mem, realtype saved_t)
3299 {
3300   int j, k;
3301   
3302   tn = saved_t;
3303   for (k = 1; k <= q; k++)
3304     for (j = q; j >= k; j--)
3305       N_VLinearSum(ONE, zn[j-1], -ONE, zn[j], zn[j-1]);
3306 }
3307
3308 /*
3309  * CVDoErrorTest
3310  *
3311  * This routine performs the local error test. 
3312  * The weighted local error norm dsm is loaded into *dsmPtr, and 
3313  * the test dsm ?<= 1 is made.
3314  *
3315  * If the test passes, CVDoErrorTest returns CV_SUCCESS. 
3316  *
3317  * If the test fails, we undo the step just taken (call CVRestore) and 
3318  *
3319  *   - if maxnef error test failures have occurred or if ABS(h) = hmin,
3320  *     we return CV_ERR_FAILURE.
3321  *
3322  *   - if more than MXNEF1 error test failures have occurred, an order
3323  *     reduction is forced. If already at order 1, restart by reloading 
3324  *     zn from scratch. If f() fails we return either CV_RHSFUNC_FAIL
3325  *     or CV_UNREC_RHSFUNC_ERR (no recovery is possible at this stage).
3326  *
3327  *   - otherwise, set *nflagPtr to PREV_ERR_FAIL, and return TRY_AGAIN. 
3328  *
3329  */
3330
3331 static booleantype CVDoErrorTest(CVodeMem cv_mem, int *nflagPtr,
3332                                 realtype saved_t, int *nefPtr, realtype *dsmPtr)
3333 {
3334   realtype dsm;
3335   int retval;
3336   
3337   dsm = acnrm / tq[2];
3338
3339   /* If est. local error norm dsm passes test, return CV_SUCCESS */  
3340   *dsmPtr = dsm; 
3341   if (dsm <= ONE) return(CV_SUCCESS);
3342   
3343   /* Test failed; increment counters, set nflag, and restore zn array */
3344   (*nefPtr)++;
3345   netf++;
3346   *nflagPtr = PREV_ERR_FAIL;
3347   CVRestore(cv_mem, saved_t);
3348
3349   /* At maxnef failures or |h| = hmin, return CV_ERR_FAILURE */
3350   if ((ABS(h) <= hmin*ONEPSM) || (*nefPtr == maxnef)) return(CV_ERR_FAILURE);
3351
3352   /* Set etamax = 1 to prevent step size increase at end of this step */
3353   etamax = ONE;
3354
3355   /* Set h ratio eta from dsm, rescale, and return for retry of step */
3356   if (*nefPtr <= MXNEF1) {
3357     eta = ONE / (RPowerR(BIAS2*dsm,ONE/L) + ADDON);
3358     eta = MAX(ETAMIN, MAX(eta, hmin / ABS(h)));
3359     if (*nefPtr >= SMALL_NEF) eta = MIN(eta, ETAMXF);
3360     CVRescale(cv_mem);
3361     return(TRY_AGAIN);
3362   }
3363   
3364   /* After MXNEF1 failures, force an order reduction and retry step */
3365   if (q > 1) {
3366     eta = MAX(ETAMIN, hmin / ABS(h));
3367     CVAdjustOrder(cv_mem,-1);
3368     L = q;
3369     q--;
3370     qwait = L;
3371     CVRescale(cv_mem);
3372     return(TRY_AGAIN);
3373   }
3374
3375   /* If already at order 1, restart: reload zn from scratch */
3376
3377   eta = MAX(ETAMIN, hmin / ABS(h));
3378   h *= eta;
3379   next_h = h;
3380   hscale = h;
3381   qwait = LONG_WAIT;
3382   nscon = 0;
3383
3384   retval = f(tn, zn[0], tempv, f_data);
3385   nfe++;
3386   if (retval < 0)  return(CV_RHSFUNC_FAIL);
3387   if (retval > 0)  return(CV_UNREC_RHSFUNC_ERR);
3388
3389   N_VScale(h, tempv, zn[1]);
3390
3391   return(TRY_AGAIN);
3392 }
3393
3394 /* 
3395  * =================================================================
3396  *  Private Functions Implementation after succesful step
3397  * =================================================================
3398  */
3399
3400 /*
3401  * CVCompleteStep
3402  *
3403  * This routine performs various update operations when the solution
3404  * to the nonlinear system has passed the local error test. 
3405  * We increment the step counter nst, record the values hu and qu,
3406  * update the tau array, and apply the corrections to the zn array.
3407  * The tau[i] are the last q values of h, with tau[1] the most recent.
3408  * The counter qwait is decremented, and if qwait == 1 (and q < qmax)
3409  * we save acor and tq[5] for a possible order increase.
3410  */
3411
3412 static void CVCompleteStep(CVodeMem cv_mem)
3413 {
3414   int i, j;
3415   
3416   nst++;
3417   nscon++;
3418   hu = h;
3419   qu = q;
3420
3421   for (i=q; i >= 2; i--)  tau[i] = tau[i-1];
3422   if ((q==1) && (nst > 1)) tau[2] = tau[1];
3423   tau[1] = h;
3424
3425   for (j=0; j <= q; j++) 
3426     N_VLinearSum(l[j], acor, ONE, zn[j], zn[j]);
3427   qwait--;
3428   if ((qwait == 1) && (q != qmax)) {
3429     N_VScale(ONE, acor, zn[qmax]);
3430     saved_tq5 = tq[5];
3431     indx_acor = qmax;
3432   }
3433 }
3434
3435 /*
3436  * CVprepareNextStep
3437  *
3438  * This routine handles the setting of stepsize and order for the
3439  * next step -- hprime and qprime.  Along with hprime, it sets the
3440  * ratio eta = hprime/h.  It also updates other state variables 
3441  * related to a change of step size or order. 
3442  */
3443
3444  static void CVPrepareNextStep(CVodeMem cv_mem, realtype dsm)
3445 {
3446   /* If etamax = 1, defer step size or order changes */
3447   if (etamax == ONE) {
3448     qwait = MAX(qwait, 2);
3449     qprime = q;
3450     hprime = h;
3451     eta = ONE;
3452     return;
3453   }
3454
3455   /* etaq is the ratio of new to old h at the current order */  
3456   etaq = ONE /(RPowerR(BIAS2*dsm,ONE/L) + ADDON);
3457   
3458   /* If no order change, adjust eta and acor in CVSetEta and return */
3459   if (qwait != 0) {
3460     eta = etaq;
3461     qprime = q;
3462     CVSetEta(cv_mem);
3463     return;
3464   }
3465   
3466   /* If qwait = 0, consider an order change.   etaqm1 and etaqp1 are 
3467      the ratios of new to old h at orders q-1 and q+1, respectively.
3468      CVChooseEta selects the largest; CVSetEta adjusts eta and acor */
3469   qwait = 2;
3470   etaqm1 = CVComputeEtaqm1(cv_mem);
3471   etaqp1 = CVComputeEtaqp1(cv_mem);  
3472   CVChooseEta(cv_mem); 
3473   CVSetEta(cv_mem);
3474 }
3475
3476 /*
3477  * CVsetEta
3478  *
3479  * This routine adjusts the value of eta according to the various
3480  * heuristic limits and the optional input hmax.  It also resets
3481  * etamax to be the estimated local error vector.
3482  */
3483
3484 static void CVSetEta(CVodeMem cv_mem)
3485 {
3486
3487   /* If eta below the threshhold THRESH, reject a change of step size */
3488   if (eta < THRESH) {
3489     eta = ONE;
3490     hprime = h;
3491   } else {
3492     /* Limit eta by etamax and hmax, then set hprime */
3493     eta = MIN(eta, etamax);
3494     eta /= MAX(ONE, ABS(h)*hmax_inv*eta);
3495     hprime = h * eta;
3496     if (qprime < q) nscon = 0;
3497   }
3498   
3499   /* Reset etamax for the next step size change, and scale acor */
3500 }
3501
3502 /*
3503  * CVComputeEtaqm1
3504  *
3505  * This routine computes and returns the value of etaqm1 for a
3506  * possible decrease in order by 1.
3507  */
3508
3509 static realtype CVComputeEtaqm1(CVodeMem cv_mem)
3510 {
3511   realtype ddn;
3512   
3513   etaqm1 = ZERO;
3514   if (q > 1) {
3515     ddn = N_VWrmsNorm(zn[q], ewt) / tq[1];
3516     etaqm1 = ONE/(RPowerR(BIAS1*ddn, ONE/q) + ADDON);
3517   }
3518   return(etaqm1);
3519 }
3520
3521 /*
3522  * CVComputeEtaqp1
3523  *
3524  * This routine computes and returns the value of etaqp1 for a
3525  * possible increase in order by 1.
3526  */
3527
3528 static realtype CVComputeEtaqp1(CVodeMem cv_mem)
3529 {
3530   realtype dup, cquot;
3531   
3532   etaqp1 = ZERO;
3533   if (q != qmax) {
3534     cquot = (tq[5] / saved_tq5) * RPowerI(h/tau[2], L);
3535     N_VLinearSum(-cquot, zn[qmax], ONE, acor, tempv);
3536     dup = N_VWrmsNorm(tempv, ewt) /tq[3];
3537     etaqp1 = ONE / (RPowerR(BIAS3*dup, ONE/(L+1)) + ADDON);
3538   }
3539   return(etaqp1);
3540 }
3541
3542 /*
3543  * CVChooseEta
3544  * Given etaqm1, etaq, etaqp1 (the values of eta for qprime =
3545  * q - 1, q, or q + 1, respectively), this routine chooses the 
3546  * maximum eta value, sets eta to that value, and sets qprime to the
3547  * corresponding value of q.  If there is a tie, the preference
3548  * order is to (1) keep the same order, then (2) decrease the order,
3549  * and finally (3) increase the order.  If the maximum eta value
3550  * is below the threshhold THRESH, the order is kept unchanged and
3551  * eta is set to 1.
3552  */
3553
3554 static void CVChooseEta(CVodeMem cv_mem)
3555 {
3556   realtype etam;
3557   
3558   etam = MAX(etaqm1, MAX(etaq, etaqp1));
3559   
3560   if (etam < THRESH) {
3561     eta = ONE;
3562     qprime = q;
3563     return;
3564   }
3565
3566   if (etam == etaq) {
3567
3568     eta = etaq;
3569     qprime = q;
3570
3571   } else if (etam == etaqm1) {
3572
3573     eta = etaqm1;
3574     qprime = q - 1;
3575
3576   } else {
3577
3578     eta = etaqp1;
3579     qprime = q + 1;
3580
3581     if (lmm == CV_BDF) {
3582
3583       /* 
3584        * Store Delta_n in zn[qmax] to be used in order increase 
3585        *
3586        * This happens at the last step of order q before an increase
3587        * to order q+1, so it represents Delta_n in the ELTE at q+1
3588        */
3589
3590       N_VScale(ONE, acor, zn[qmax]);
3591
3592     }
3593
3594   }
3595
3596 }
3597
3598 /*
3599  * CVHandleFailure
3600  *
3601  * This routine prints error messages for all cases of failure by
3602  * CVHin and CVStep. It returns to CVode the value that CVode is 
3603  * to return to the user.
3604  */
3605
3606 static int CVHandleFailure(CVodeMem cv_mem, int flag)
3607 {
3608
3609   /* Set vector of  absolute weighted local errors */
3610   /*
3611   N_VProd(acor, ewt, tempv);
3612   N_VAbs(tempv, tempv);
3613   */
3614
3615   /* Depending on kflag, print error message and return error flag */
3616   switch (flag) {
3617   case CV_ERR_FAILURE: 
3618     CVProcessError(cv_mem, CV_ERR_FAILURE, "CVODE", "CVode", MSGCV_ERR_FAILS, tn, h);
3619     break;
3620   case CV_CONV_FAILURE:
3621     CVProcessError(cv_mem, CV_CONV_FAILURE, "CVODE", "CVode", MSGCV_CONV_FAILS, tn, h);
3622     break;
3623   case CV_LSETUP_FAIL:
3624     CVProcessError(cv_mem, CV_LSETUP_FAIL, "CVODE", "CVode", MSGCV_SETUP_FAILED, tn);
3625     break;
3626   case CV_LSOLVE_FAIL:
3627     CVProcessError(cv_mem, CV_LSOLVE_FAIL, "CVODE", "CVode", MSGCV_SOLVE_FAILED, tn);
3628     break;
3629   case CV_RHSFUNC_FAIL:
3630     CVProcessError(cv_mem, CV_RHSFUNC_FAIL, "CVODE", "CVode", MSGCV_RHSFUNC_FAILED, tn);
3631     break;
3632   case CV_UNREC_RHSFUNC_ERR:
3633     CVProcessError(cv_mem, CV_UNREC_RHSFUNC_ERR, "CVODE", "CVode", MSGCV_RHSFUNC_UNREC, tn);
3634     break;
3635   case CV_REPTD_RHSFUNC_ERR:
3636     CVProcessError(cv_mem, CV_REPTD_RHSFUNC_ERR, "CVODE", "CVode", MSGCV_RHSFUNC_REPTD, tn);
3637     break;
3638   case CV_RTFUNC_FAIL:    
3639     CVProcessError(cv_mem, CV_RTFUNC_FAIL, "CVODE", "CVode", MSGCV_RTFUNC_FAILED, tn);
3640     break;
3641   case CV_TOO_CLOSE:
3642     CVProcessError(cv_mem, CV_TOO_CLOSE, "CVODE", "CVode", MSGCV_TOO_CLOSE);
3643   default:
3644     return(CV_SUCCESS);   
3645   }
3646
3647   return(flag);
3648
3649 }
3650
3651 /* 
3652  * =================================================================
3653  * BDF Stability Limit Detection                       
3654  * =================================================================
3655  */
3656
3657 /*
3658  * CVBDFStab
3659  *
3660  * This routine handles the BDF Stability Limit Detection Algorithm
3661  * STALD.  It is called if lmm = CV_BDF and the SLDET option is on.
3662  * If the order is 3 or more, the required norm data is saved.
3663  * If a decision to reduce order has not already been made, and
3664  * enough data has been saved, CVsldet is called.  If it signals
3665  * a stability limit violation, the order is reduced, and the step
3666  * size is reset accordingly.
3667  */
3668
3669 void CVBDFStab(CVodeMem cv_mem)
3670 {
3671   int i,k, ldflag, factorial;
3672   realtype sq, sqm1, sqm2;
3673       
3674   /* If order is 3 or greater, then save scaled derivative data,
3675      push old data down in i, then add current values to top.    */
3676
3677   if (q >= 3) {
3678     for (k = 1; k <= 3; k++)
3679       { for (i = 5; i >= 2; i--) ssdat[i][k] = ssdat[i-1][k]; }
3680     factorial = 1;
3681     for (i = 1; i <= q-1; i++) factorial *= i;
3682     sq = factorial*q*(q+1)*acnrm/tq[5];
3683     sqm1 = factorial*q*N_VWrmsNorm(zn[q], ewt);
3684     sqm2 = factorial*N_VWrmsNorm(zn[q-1], ewt);
3685     ssdat[1][1] = sqm2*sqm2;
3686     ssdat[1][2] = sqm1*sqm1;
3687     ssdat[1][3] = sq*sq;
3688   }  
3689
3690   if (qprime >= q) {
3691
3692     /* If order is 3 or greater, and enough ssdat has been saved,
3693        nscon >= q+5, then call stability limit detection routine.  */
3694
3695     if ( (q >= 3) && (nscon >= q+5) ) {
3696       ldflag = CVsldet(cv_mem);
3697       if (ldflag > 3) {
3698         /* A stability limit violation is indicated by
3699            a return flag of 4, 5, or 6.
3700            Reduce new order.                     */
3701         qprime = q-1;
3702         eta = etaqm1; 
3703         eta = MIN(eta,etamax);
3704         eta = eta/MAX(ONE,ABS(h)*hmax_inv*eta);
3705         hprime = h*eta;
3706         nor = nor + 1;
3707       }
3708     }
3709   }
3710   else {
3711     /* Otherwise, let order increase happen, and 
3712        reset stability limit counter, nscon.     */
3713     nscon = 0;
3714   }
3715 }
3716
3717 /*
3718  * CVsldet
3719  *
3720  * This routine detects stability limitation using stored scaled 
3721  * derivatives data. CVsldet returns the magnitude of the
3722  * dominate characteristic root, rr. The presents of a stability
3723  * limit is indicated by rr > "something a little less then 1.0",  
3724  * and a positive kflag. This routine should only be called if
3725  * order is greater than or equal to 3, and data has been collected
3726  * for 5 time steps. 
3727  * 
3728  * Returned values:
3729  *    kflag = 1 -> Found stable characteristic root, normal matrix case
3730  *    kflag = 2 -> Found stable characteristic root, quartic solution
3731  *    kflag = 3 -> Found stable characteristic root, quartic solution,
3732  *                 with Newton correction
3733  *    kflag = 4 -> Found stability violation, normal matrix case
3734  *    kflag = 5 -> Found stability violation, quartic solution
3735  *    kflag = 6 -> Found stability violation, quartic solution,
3736  *                 with Newton correction
3737  *
3738  *    kflag < 0 -> No stability limitation, 
3739  *                 or could not compute limitation.
3740  *
3741  *    kflag = -1 -> Min/max ratio of ssdat too small.
3742  *    kflag = -2 -> For normal matrix case, vmax > vrrt2*vrrt2
3743  *    kflag = -3 -> For normal matrix case, The three ratios
3744  *                  are inconsistent.
3745  *    kflag = -4 -> Small coefficient prevents elimination of quartics.  
3746  *    kflag = -5 -> R value from quartics not consistent.
3747  *    kflag = -6 -> No corrected root passes test on qk values
3748  *    kflag = -7 -> Trouble solving for sigsq.
3749  *    kflag = -8 -> Trouble solving for B, or R via B.
3750  *    kflag = -9 -> R via sigsq[k] disagrees with R from data.
3751  */
3752
3753 static int CVsldet(CVodeMem cv_mem)
3754 {
3755   int i, k, j, it, kmin, kflag = 0;
3756   realtype rat[5][4], rav[4], qkr[4], sigsq[4], smax[4], ssmax[4];
3757   realtype drr[4], rrc[4],sqmx[4], qjk[4][4], vrat[5], qc[6][4], qco[6][4];
3758   realtype rr, rrcut, vrrtol, vrrt2, sqtol, rrtol;
3759   realtype smink, smaxk, sumrat, sumrsq, vmin, vmax, drrmax, adrr;
3760   realtype tem, sqmax, saqk, qp, s, sqmaxk, saqj, sqmin;
3761   realtype rsa, rsb, rsc, rsd, rd1a, rd1b, rd1c;
3762   realtype rd2a, rd2b, rd3a, cest1, corr1; 
3763   realtype ratp, ratm, qfac1, qfac2, bb, rrb;
3764
3765   /* The following are cutoffs and tolerances used by this routine */
3766
3767   rrcut  = RCONST(0.98);
3768   vrrtol = RCONST(1.0e-4);
3769   vrrt2  = RCONST(5.0e-4);
3770   sqtol  = RCONST(1.0e-3);
3771   rrtol  = RCONST(1.0e-2);
3772   
3773   rr = ZERO;
3774   
3775   /*  Index k corresponds to the degree of the interpolating polynomial. */
3776   /*      k = 1 -> q-1          */
3777   /*      k = 2 -> q            */
3778   /*      k = 3 -> q+1          */
3779   
3780   /*  Index i is a backward-in-time index, i = 1 -> current time, */
3781   /*      i = 2 -> previous step, etc    */
3782   
3783   /* get maxima, minima, and variances, and form quartic coefficients  */
3784   
3785   for (k=1; k<=3; k++) {
3786     smink = ssdat[1][k];
3787     smaxk = ZERO;
3788     
3789     for (i=1; i<=5; i++) {
3790       smink = MIN(smink,ssdat[i][k]);
3791       smaxk = MAX(smaxk,ssdat[i][k]);
3792     }
3793     
3794     if (smink < TINY*smaxk) {
3795       kflag = -1;  
3796       return(kflag);
3797     }
3798     smax[k] = smaxk;
3799     ssmax[k] = smaxk*smaxk;
3800     
3801     sumrat = ZERO;
3802     sumrsq = ZERO;
3803     for (i=1; i<=4; i++) {
3804       rat[i][k] = ssdat[i][k]/ssdat[i+1][k];
3805       sumrat = sumrat + rat[i][k];
3806       sumrsq = sumrsq + rat[i][k]*rat[i][k];
3807     } 
3808     rav[k] = FOURTH*sumrat;
3809     vrat[k] = ABS(FOURTH*sumrsq - rav[k]*rav[k]);
3810     
3811     qc[5][k] = ssdat[1][k]*ssdat[3][k] - ssdat[2][k]*ssdat[2][k];
3812     qc[4][k] = ssdat[2][k]*ssdat[3][k] - ssdat[1][k]*ssdat[4][k];
3813     qc[3][k] = ZERO;
3814     qc[2][k] = ssdat[2][k]*ssdat[5][k] - ssdat[3][k]*ssdat[4][k];
3815     qc[1][k] = ssdat[4][k]*ssdat[4][k] - ssdat[3][k]*ssdat[5][k];
3816     
3817     for (i=1; i<=5; i++) {
3818       qco[i][k] = qc[i][k];
3819     }
3820   }                            /* End of k loop */
3821   
3822   /* Isolate normal or nearly-normal matrix case. Three quartic will
3823      have common or nearly-common roots in this case. 
3824      Return a kflag = 1 if this procedure works. If three root 
3825      differ more than vrrt2, return error kflag = -3.    */
3826   
3827   vmin = MIN(vrat[1],MIN(vrat[2],vrat[3]));
3828   vmax = MAX(vrat[1],MAX(vrat[2],vrat[3]));
3829   
3830   if(vmin < vrrtol*vrrtol) {
3831     if (vmax > vrrt2*vrrt2) {
3832       kflag = -2;  
3833       return(kflag);
3834     } else {
3835       rr = (rav[1] + rav[2] + rav[3])/THREE;
3836       
3837       drrmax = ZERO;
3838       for(k = 1;k<=3;k++) {
3839         adrr = ABS(rav[k] - rr);
3840         drrmax = MAX(drrmax, adrr);
3841       }
3842       if (drrmax > vrrt2) {
3843         kflag = -3;    
3844       }
3845       
3846       kflag = 1;
3847
3848       /*  can compute charactistic root, drop to next section   */
3849       
3850     }
3851   } else {
3852
3853     /* use the quartics to get rr. */
3854     
3855     if (ABS(qco[1][1]) < TINY*ssmax[1]) {
3856       kflag = -4;    
3857       return(kflag);
3858     }
3859     
3860     tem = qco[1][2]/qco[1][1];
3861     for(i=2; i<=5; i++) {
3862       qco[i][2] = qco[i][2] - tem*qco[i][1];
3863     }
3864
3865     qco[1][2] = ZERO;
3866     tem = qco[1][3]/qco[1][1];
3867     for(i=2; i<=5; i++) {
3868       qco[i][3] = qco[i][3] - tem*qco[i][1];
3869     }
3870     qco[1][3] = ZERO;
3871     
3872     if (ABS(qco[2][2]) < TINY*ssmax[2]) {
3873       kflag = -4;    
3874       return(kflag);
3875     }
3876     
3877     tem = qco[2][3]/qco[2][2];
3878     for(i=3; i<=5; i++) {
3879       qco[i][3] = qco[i][3] - tem*qco[i][2];
3880     }
3881     
3882     if (ABS(qco[4][3]) < TINY*ssmax[3]) {
3883       kflag = -4;    
3884       return(kflag);
3885     }
3886     
3887     rr = -qco[5][3]/qco[4][3];
3888     
3889     if (rr < TINY || rr > HUN) {
3890       kflag = -5;   
3891       return(kflag);
3892     }
3893     
3894     for(k=1; k<=3; k++) {
3895       qkr[k] = qc[5][k] + rr*(qc[4][k] + rr*rr*(qc[2][k] + rr*qc[1][k]));
3896     }  
3897     
3898     sqmax = ZERO;
3899     for(k=1; k<=3; k++) {
3900       saqk = ABS(qkr[k])/ssmax[k];
3901       if (saqk > sqmax) sqmax = saqk;
3902     } 
3903     
3904     if (sqmax < sqtol) {
3905       kflag = 2;
3906       
3907       /*  can compute charactistic root, drop to "given rr,etc"   */
3908       
3909     } else {
3910
3911       /* do Newton corrections to improve rr.  */
3912       
3913       for(it=1; it<=3; it++) {
3914         for(k=1; k<=3; k++) {
3915           qp = qc[4][k] + rr*rr*(THREE*qc[2][k] + rr*FOUR*qc[1][k]);
3916           drr[k] = ZERO;
3917           if (ABS(qp) > TINY*ssmax[k]) drr[k] = -qkr[k]/qp;
3918           rrc[k] = rr + drr[k];
3919         } 
3920         
3921         for(k=1; k<=3; k++) {
3922           s = rrc[k];
3923           sqmaxk = ZERO;
3924           for(j=1; j<=3; j++) {
3925             qjk[j][k] = qc[5][j] + s*(qc[4][j] + 
3926                                       s*s*(qc[2][j] + s*qc[1][j]));
3927             saqj = ABS(qjk[j][k])/ssmax[j];
3928             if (saqj > sqmaxk) sqmaxk = saqj;
3929           } 
3930           sqmx[k] = sqmaxk;
3931         } 
3932
3933         sqmin = sqmx[1]; kmin = 1;
3934         for(k=2; k<=3; k++) {
3935           if (sqmx[k] < sqmin) {
3936             kmin = k;
3937             sqmin = sqmx[k];
3938           }
3939         } 
3940         rr = rrc[kmin];
3941         
3942         if (sqmin < sqtol) {
3943           kflag = 3;
3944           /*  can compute charactistic root   */
3945           /*  break out of Newton correction loop and drop to "given rr,etc" */ 
3946           break;
3947         } else {
3948           for(j=1; j<=3; j++) {
3949             qkr[j] = qjk[j][kmin];
3950           }
3951         }     
3952       }          /*  end of Newton correction loop  */ 
3953       
3954       if (sqmin > sqtol) {
3955         kflag = -6;
3956         return(kflag);
3957       }
3958     }     /*  end of if (sqmax < sqtol) else   */
3959   }      /*  end of if(vmin < vrrtol*vrrtol) else, quartics to get rr. */
3960   
3961   /* given rr, find sigsq[k] and verify rr.  */
3962   /* All positive kflag drop to this section  */
3963   
3964   for(k=1; k<=3; k++) {
3965     rsa = ssdat[1][k];
3966     rsb = ssdat[2][k]*rr;
3967     rsc = ssdat[3][k]*rr*rr;
3968     rsd = ssdat[4][k]*rr*rr*rr;
3969     rd1a = rsa - rsb;
3970     rd1b = rsb - rsc;
3971     rd1c = rsc - rsd;
3972     rd2a = rd1a - rd1b;
3973     rd2b = rd1b - rd1c;
3974     rd3a = rd2a - rd2b;
3975     
3976     if (ABS(rd1b) < TINY*smax[k]) {
3977       kflag = -7;
3978       return(kflag);
3979     }
3980     
3981     cest1 = -rd3a/rd1b;
3982     if (cest1 < TINY || cest1 > FOUR) {
3983       kflag = -7;
3984       return(kflag);
3985     }
3986     corr1 = (rd2b/cest1)/(rr*rr);
3987     sigsq[k] = ssdat[3][k] + corr1;
3988   }
3989   
3990   if (sigsq[2] < TINY) {
3991     kflag = -8;
3992     return(kflag);
3993   }
3994   
3995   ratp = sigsq[3]/sigsq[2];
3996   ratm = sigsq[1]/sigsq[2];
3997   qfac1 = FOURTH*(q*q - ONE);
3998   qfac2 = TWO/(q - ONE);
3999   bb = ratp*ratm - ONE - qfac1*ratp;
4000   tem = ONE - qfac2*bb;
4001   
4002   if (ABS(tem) < TINY) {
4003     kflag = -8;
4004     return(kflag);
4005   }
4006   
4007   rrb = ONE/tem;
4008   
4009   if (ABS(rrb - rr) > rrtol) {
4010     kflag = -9;
4011     return(kflag);
4012   }
4013   
4014   /* Check to see if rr is above cutoff rrcut  */
4015   if (rr > rrcut) {
4016     if (kflag == 1) kflag = 4;
4017     if (kflag == 2) kflag = 5;
4018     if (kflag == 3) kflag = 6;
4019   }
4020   
4021   /* All positive kflag returned at this point  */
4022   
4023   return(kflag);
4024   
4025 }
4026
4027 /* 
4028  * =================================================================
4029  * Root finding   
4030  * =================================================================
4031  */
4032
4033 /*-----------------------------------------------------------------*/
4034
4035 /* 
4036  * CVRcheck1
4037  *
4038  * This routine completes the initialization of rootfinding memory
4039  * information, and checks whether g has a zero both at and very near
4040  * the initial point of the IVP.
4041  *
4042  * This routine returns an int equal to:
4043  *  INITROOT = -1 if a close pair of zeros was found, and
4044  *  CV_SUCCESS     =  0 otherwise.
4045  */
4046
4047 static int CVRcheck1Std(CVodeMem cv_mem)
4048 {
4049   int i, retval;
4050   realtype smallh, hratio;
4051   booleantype zroot;
4052
4053   for (i = 0; i < nrtfn; i++) iroots[i] = 0;
4054   tlo = tn;
4055   ttol = (ABS(tn) + ABS(h))*uround*HUN;
4056
4057   /* Evaluate g at initial t and check for zero values. */
4058   retval = gfun(tlo, zn[0], glo, g_data);
4059   nge = 1;
4060   if (retval != 0) return(CV_RTFUNC_FAIL);
4061
4062   zroot = FALSE;
4063   for (i = 0; i < nrtfn; i++) {
4064     if (ABS(glo[i]) == ZERO) zroot = TRUE;
4065   }
4066   if (!zroot) return(CV_SUCCESS);
4067
4068   /* Some g_i is zero at t0; look at g at t0+(small increment). */
4069   hratio = MAX(ttol/ABS(h), TENTH);
4070   smallh = hratio*h;
4071   tlo += smallh;
4072   N_VLinearSum(ONE, zn[0], hratio, zn[1], y);
4073   retval = gfun(tlo, y, glo, g_data);
4074   nge++;
4075   if (retval != 0) return(CV_RTFUNC_FAIL);
4076
4077   zroot = FALSE;
4078   for (i = 0; i < nrtfn; i++) {
4079     if (ABS(glo[i]) == ZERO) {
4080       zroot = TRUE;
4081       iroots[i] = 1;
4082     }
4083   }
4084   if (zroot) return(INITROOT);
4085   return(CV_SUCCESS);
4086
4087 }
4088
4089 /*
4090  * CVRcheck2
4091  *
4092  * This routine checks for exact zeros of g at the last root found,
4093  * if the last return was a root.  It then checks for a close
4094  * pair of zeros (an error condition), and for a new root at a
4095  * nearby point.  The left endpoint (tlo) of the search interval
4096  * is adjusted if necessary to assure that all g_i are nonzero
4097  * there, before returning to do a root search in the interval.
4098  *
4099  * On entry, tlo = tretlast is the last value of tret returned by
4100  * CVode.  This may be the previous tn, the previous tout value, or
4101  * the last root location.
4102  *
4103  * This routine returns an int equal to:
4104  *      CLOSERT = -2 if a close pair of zeros was found,
4105  *      RTFOUND =  1 if a new zero of g was found near tlo, or
4106  *      CV_SUCCESS    =  0 otherwise.
4107  */
4108
4109 static int CVRcheck2Std(CVodeMem cv_mem)
4110 {
4111   int i, retval;
4112   realtype smallh, hratio;
4113   booleantype zroot;
4114
4115   if (irfnd == 0) return(CV_SUCCESS);
4116
4117   (void) CVodeGetDky(cv_mem, tlo, 0, y);
4118   retval = gfun(tlo, y, glo, g_data);
4119   nge++;
4120   if (retval != 0) return(CV_RTFUNC_FAIL);
4121
4122   zroot = FALSE;
4123   for (i = 0; i < nrtfn; i++) iroots[i] = 0;
4124   for (i = 0; i < nrtfn; i++) {
4125     if (ABS(glo[i]) == ZERO) {
4126       zroot = TRUE;
4127       iroots[i] = 1;
4128     }
4129   }
4130   if (!zroot) return(CV_SUCCESS);
4131
4132   /* One or more g_i has a zero at tlo.  Check g at tlo+smallh. */
4133   ttol = (ABS(tn) + ABS(h))*uround*HUN;
4134   smallh = (h > ZERO) ? ttol : -ttol;
4135   tlo += smallh;
4136   if ( (tlo - tn)*h >= ZERO) {
4137     hratio = smallh/h;
4138     N_VLinearSum(ONE, y, hratio, zn[1], y);
4139   } else {
4140     (void) CVodeGetDky(cv_mem, tlo, 0, y);
4141   }
4142   retval = gfun(tlo, y, glo, g_data);
4143   nge++;
4144   if (retval != 0) return(CV_RTFUNC_FAIL);
4145
4146   zroot = FALSE;
4147   for (i = 0; i < nrtfn; i++) {
4148     if (ABS(glo[i]) == ZERO) {
4149       if (iroots[i] == 1) return(CLOSERT);
4150       zroot = TRUE;
4151       iroots[i] = 1;
4152     }
4153   }
4154   if (zroot) return(RTFOUND);
4155   return(CV_SUCCESS);
4156
4157 }
4158
4159 /*
4160  * CVRcheck3
4161  *
4162  * This routine interfaces to CVRootfind to look for a root of g
4163  * between tlo and either tn or tout, whichever comes first.
4164  * Only roots beyond tlo in the direction of integration are sought.
4165  *
4166  * This routine returns an int equal to:
4167  *      RTFOUND =  1 if a root of g was found, or
4168  *      CV_SUCCESS    =  0 otherwise.
4169  */
4170
4171 static int CVRcheck3Std(CVodeMem cv_mem)
4172 {
4173   int i, retval, ier;
4174
4175   /* Set thi = tn or tout, whichever comes first; set y = y(thi). */
4176   if (taskc == CV_ONE_STEP) {
4177     thi = tn;
4178     N_VScale(ONE, zn[0], y);
4179   }
4180   if (taskc == CV_NORMAL) {
4181     if ( (toutc - tn)*h >= ZERO) {
4182       thi = tn; 
4183       N_VScale(ONE, zn[0], y);
4184     } else {
4185       thi = toutc;
4186       (void) CVodeGetDky(cv_mem, thi, 0, y);
4187     }
4188   }
4189
4190   /* Set ghi = g(thi) and call CVRootfind to search (tlo,thi) for roots. */
4191   retval = gfun(thi, y, ghi, g_data);
4192   nge++;
4193   if (retval != 0) return(CV_RTFUNC_FAIL);
4194
4195   ttol = (ABS(tn) + ABS(h))*uround*HUN;
4196   ier = CVRootfind(cv_mem);
4197   tlo = trout;
4198   for (i = 0; i < nrtfn; i++) glo[i] = grout[i];
4199
4200   /* If no root found, return CV_SUCCESS. */  
4201   if (ier == CV_SUCCESS) return(CV_SUCCESS);
4202
4203   /* If a root was found, interpolate to get y(trout) and return.  */
4204   (void) CVodeGetDky(cv_mem, trout, 0, y);
4205   return(RTFOUND);
4206
4207 }
4208
4209 /*
4210  * CVRootFind
4211  *
4212  * This routine solves for a root of g(t) between tlo and thi, if
4213  * one exists.  Only roots of odd multiplicity (i.e. with a change
4214  * of sign in one of the g_i), or exact zeros, are found.
4215  * Here the sign of tlo - thi is arbitrary, but if multiple roots
4216  * are found, the one closest to tlo is returned.
4217  *
4218  * The method used is the Illinois algorithm, a modified secant method.
4219  * Reference: Kathie L. Hiebert and Lawrence F. Shampine, Implicitly
4220  * Defined Output Points for Solutions of ODEs, Sandia National
4221  * Laboratory Report SAND80-0180, February 1980.
4222  *
4223  * This routine uses the following parameters for communication:
4224  *
4225  * nrtfn    = number of functions g_i, or number of components of
4226  *            the vector-valued function g(t).  Input only.
4227  *
4228  * gfun     = user-defined function for g(t).  Its form is
4229  *            (void) gfun(t, y, gt, g_data)
4230  *
4231  * nge      = cumulative counter for gfun calls.
4232  *
4233  * ttol     = a convergence tolerance for trout.  Input only.
4234  *            When a root at trout is found, it is located only to
4235  *            within a tolerance of ttol.  Typically, ttol should
4236  *            be set to a value on the order of
4237  *               100 * UROUND * max (ABS(tlo), ABS(thi))
4238  *            where UROUND is the unit roundoff of the machine.
4239  *
4240  * tlo, thi = endpoints of the interval in which roots are sought.
4241  *            On input, and must be distinct, but tlo - thi may
4242  *            be of either sign.  The direction of integration is
4243  *            assumed to be from tlo to thi.  On return, tlo and thi
4244  *            are the endpoints of the final relevant interval.
4245  *
4246  * glo, ghi = arrays of length nrtfn containing the vectors g(tlo)
4247  *            and g(thi) respectively.  Input and output.  On input,
4248  *            none of the glo[i] should be zero.
4249  *
4250  * trout    = root location, if a root was found, or thi if not.
4251  *            Output only.  If a root was found other than an exact
4252  *            zero of g, trout is the endpoint thi of the final
4253  *            interval bracketing the root, with size at most ttol.
4254  *
4255  * grout    = array of length nrtfn containing g(trout) on return.
4256  *
4257  * iroots   = int array of length nrtfn with root information.
4258  *            Output only.  If a root was found, iroots indicates
4259  *            which components g_i have a root at trout.  For
4260  *            i = 0, ..., nrtfn-1, iroots[i] = 1 if g_i has a root
4261  *            and iroots[i] = 0 otherwise.
4262  *
4263  * This routine returns an int equal to:
4264  *      RTFOUND =  1 if a root of g was found, or
4265  *      CV_SUCCESS    =  0 otherwise.
4266  */
4267
4268 static int CVRootfindStd(CVodeMem cv_mem)
4269 {
4270   realtype alpha, tmid, gfrac, maxfrac, fracint, fracsub;
4271   int i, retval, imax, side, sideprev;
4272   booleantype zroot, sgnchg;
4273
4274   imax = 0;
4275
4276   /* First check for change in sign in ghi or for a zero in ghi. */
4277   maxfrac = ZERO;
4278   zroot = FALSE;
4279   sgnchg = FALSE;
4280   for (i = 0;  i < nrtfn; i++) {
4281     if (ABS(ghi[i]) == ZERO) {
4282       zroot = TRUE;
4283     } else {
4284       if (glo[i]*ghi[i] < ZERO) {
4285         gfrac = ABS(ghi[i]/(ghi[i] - glo[i]));
4286         if (gfrac > maxfrac) {
4287           sgnchg = TRUE;
4288           maxfrac = gfrac;
4289           imax = i;
4290         }
4291       }
4292     }
4293   }
4294
4295   /* If no sign change was found, reset trout and grout.  Then return
4296      CV_SUCCESS if no zero was found, or set iroots and return RTFOUND.  */ 
4297   if (!sgnchg) {
4298     trout = thi;
4299     for (i = 0; i < nrtfn; i++) grout[i] = ghi[i];
4300     if (!zroot) return(CV_SUCCESS);
4301     for (i = 0; i < nrtfn; i++) {
4302       iroots[i] = 0;
4303       if (ABS(ghi[i]) == ZERO) iroots[i] = 1;
4304     }
4305     return(RTFOUND);
4306   }
4307
4308   /* Initialize alpha to avoid compiler warning */
4309   alpha = ONE;
4310
4311   /* A sign change was found.  Loop to locate nearest root. */
4312
4313   side = 0;  sideprev = -1;
4314   for(;;) {                                    /* Looping point */
4315
4316     /* Set weight alpha.
4317        On the first two passes, set alpha = 1.  Thereafter, reset alpha
4318        according to the side (low vs high) of the subinterval in which
4319        the sign change was found in the previous two passes.
4320        If the sides were opposite, set alpha = 1.
4321        If the sides were the same, then double alpha (if high side),
4322        or halve alpha (if low side).
4323        The next guess tmid is the secant method value if alpha = 1, but
4324        is closer to tlo if alpha < 1, and closer to thi if alpha > 1.    */
4325
4326     if (sideprev == side) {
4327       alpha = (side == 2) ? alpha*TWO : alpha*HALF;
4328     } else {
4329       alpha = ONE;
4330     }
4331
4332     /* Set next root approximation tmid and get g(tmid).
4333        If tmid is too close to tlo or thi, adjust it inward,
4334        by a fractional distance that is between 0.1 and 0.5.  */
4335     tmid = thi - (thi - tlo)*ghi[imax]/(ghi[imax] - alpha*glo[imax]);
4336     if (ABS(tmid - tlo) < HALF*ttol) {
4337       fracint = ABS(thi - tlo)/ttol;
4338       fracsub = (fracint > FIVE) ? TENTH : HALF/fracint;
4339       tmid = tlo + fracsub*(thi - tlo);
4340     }
4341     if (ABS(thi - tmid) < HALF*ttol) {
4342       fracint = ABS(thi - tlo)/ttol;
4343       fracsub = (fracint > FIVE) ? TENTH : HALF/fracint;
4344       tmid = thi - fracsub*(thi - tlo);
4345     }
4346
4347     (void) CVodeGetDky(cv_mem, tmid, 0, y);
4348     retval = gfun(tmid, y, grout, g_data);
4349     nge++;
4350     if (retval != 0) return(CV_RTFUNC_FAIL);
4351
4352     /* Check to see in which subinterval g changes sign, and reset imax.
4353        Set side = 1 if sign change is on low side, or 2 if on high side.  */  
4354     maxfrac = ZERO;
4355     zroot = FALSE;
4356     sgnchg = FALSE;
4357     sideprev = side;
4358     for (i = 0;  i < nrtfn; i++) {
4359       if (ABS(grout[i]) == ZERO) {
4360         zroot = TRUE;
4361       } else {
4362         if (glo[i]*grout[i] < ZERO) {
4363           gfrac = ABS(grout[i]/(grout[i] - glo[i]));
4364           if (gfrac > maxfrac) {
4365             sgnchg = TRUE;
4366             maxfrac = gfrac;
4367             imax = i;
4368           }
4369         }
4370       }
4371     }
4372     if (sgnchg) {
4373       /* Sign change found in (tlo,tmid); replace thi with tmid. */
4374       thi = tmid;
4375       for (i = 0; i < nrtfn; i++) ghi[i] = grout[i];
4376       side = 1;
4377       /* Stop at root thi if converged; otherwise loop. */
4378       if (ABS(thi - tlo) <= ttol) break;
4379       continue;  /* Return to looping point. */
4380     }
4381
4382     if (zroot) {
4383       /* No sign change in (tlo,tmid), but g = 0 at tmid; return root tmid. */
4384       thi = tmid;
4385       for (i = 0; i < nrtfn; i++) ghi[i] = grout[i];
4386       break;
4387     }
4388
4389     /* No sign change in (tlo,tmid), and no zero at tmid.
4390        Sign change must be in (tmid,thi).  Replace tlo with tmid. */
4391     tlo = tmid;
4392     for (i = 0; i < nrtfn; i++) glo[i] = grout[i];
4393     side = 2;
4394     /* Stop at root thi if converged; otherwise loop back. */
4395     if (ABS(thi - tlo) <= ttol) break;
4396
4397   } /* End of root-search loop */
4398
4399   /* Reset trout and grout, set iroots, and return RTFOUND. */
4400   trout = thi;
4401   for (i = 0; i < nrtfn; i++) {
4402     grout[i] = ghi[i];
4403     iroots[i] = 0;
4404     if (ABS(ghi[i]) == ZERO) iroots[i] = 1;
4405     if (glo[i]*ghi[i] < ZERO) iroots[i] = 1;
4406   }
4407   return(RTFOUND);
4408 }
4409
4410 /* 
4411  * =================================================================
4412  * Internal EWT function
4413  * =================================================================
4414  */
4415
4416 /*
4417  * CVEwtSet
4418  *
4419  * This routine is responsible for setting the error weight vector ewt,
4420  * according to tol_type, as follows:
4421  *
4422  * (1) ewt[i] = 1 / (reltol * ABS(ycur[i]) + *abstol), i=0,...,neq-1
4423  *     if tol_type = CV_SS
4424  * (2) ewt[i] = 1 / (reltol * ABS(ycur[i]) + abstol[i]), i=0,...,neq-1
4425  *     if tol_type = CV_SV
4426  *
4427  * CVEwtSet returns 0 if ewt is successfully set as above to a
4428  * positive vector and -1 otherwise. In the latter case, ewt is
4429  * considered undefined.
4430  *
4431  * All the real work is done in the routines CVEwtSetSS, CVEwtSetSV.
4432  */
4433
4434 int CVEwtSet(N_Vector ycur, N_Vector weight, void *data)
4435 {
4436   CVodeMem cv_mem;
4437   int flag = 0;
4438
4439   /* data points to cv_mem here */
4440
4441   cv_mem = (CVodeMem) data;
4442
4443   switch(itol) {
4444   case CV_SS: 
4445     flag = CVEwtSetSS(cv_mem, ycur, weight);
4446     break;
4447   case CV_SV: 
4448     flag = CVEwtSetSV(cv_mem, ycur, weight);
4449     break;
4450   }
4451   
4452   return(flag);
4453 }
4454
4455 /*
4456  * CVEwtSetSS
4457  *
4458  * This routine sets ewt as decribed above in the case tol_type = CV_SS.
4459  * It tests for non-positive components before inverting. CVEwtSetSS
4460  * returns 0 if ewt is successfully set to a positive vector
4461  * and -1 otherwise. In the latter case, ewt is considered undefined.
4462  */
4463
4464 static int CVEwtSetSS(CVodeMem cv_mem, N_Vector ycur, N_Vector weight)
4465 {
4466   N_VAbs(ycur, tempv);
4467   N_VScale(reltol, tempv, tempv);
4468   N_VAddConst(tempv, Sabstol, tempv);
4469   if (N_VMin(tempv) <= ZERO) return(-1);
4470   N_VInv(tempv, weight);
4471   return(0);
4472 }
4473
4474 /*
4475  * CVEwtSetSV
4476  *
4477  * This routine sets ewt as decribed above in the case tol_type = CV_SV.
4478  * It tests for non-positive components before inverting. CVEwtSetSV
4479  * returns 0 if ewt is successfully set to a positive vector
4480  * and -1 otherwise. In the latter case, ewt is considered undefined.
4481  */
4482
4483 static int CVEwtSetSV(CVodeMem cv_mem, N_Vector ycur, N_Vector weight)
4484 {
4485   N_VAbs(ycur, tempv);
4486   N_VLinearSum(reltol, tempv, ONE, Vabstol, tempv);
4487   if (N_VMin(tempv) <= ZERO) return(-1);
4488   N_VInv(tempv, weight);
4489   return(0);
4490 }
4491
4492 /* 
4493  * =================================================================
4494  * CVODE Error Handling function   
4495  * =================================================================
4496  */
4497
4498 /* 
4499  * CVProcessError is a high level error handling function
4500  * - if cv_mem==NULL it prints the error message to stderr
4501  * - otherwise, it sets-up and calls the error hadling function 
4502  *   pointed to by cv_ehfun
4503  */
4504
4505 #define ehfun    (cv_mem->cv_ehfun)
4506 #define eh_data  (cv_mem->cv_eh_data)
4507
4508 void CVProcessError(CVodeMem cv_mem, 
4509                     int error_code, const char *module, const char *fname, 
4510                     const char *msgfmt, ...)
4511 {
4512   va_list ap;
4513   char msg[256];
4514
4515   /* Initialize the argument pointer variable 
4516      (msgfmt is the last required argument to CVProcessError) */
4517
4518   va_start(ap, msgfmt);
4519
4520   if (cv_mem == NULL) {    /* We write to stderr */
4521
4522 #ifndef NO_FPRINTF_OUTPUT
4523     fprintf(stderr, "\n[%s ERROR]  %s\n  ", module, fname);
4524     fprintf(stderr, msgfmt);
4525     fprintf(stderr, "\n\n");
4526 #endif
4527
4528   } else {                 /* We can call ehfun */
4529
4530     /* Compose the message */
4531
4532     vsprintf(msg, msgfmt, ap);
4533
4534     /* Call ehfun */
4535
4536     ehfun(error_code, module, fname, msg, eh_data);
4537
4538   }
4539
4540   /* Finalize argument processing */
4541   
4542   va_end(ap);
4543
4544   return;
4545
4546 }
4547
4548 /* CVErrHandler is the default error handling function.
4549    It sends the error message to the stream pointed to by cv_errfp */
4550
4551 #define errfp    (cv_mem->cv_errfp)
4552
4553 void CVErrHandler(int error_code, const char *module,
4554                   const char *function, char *msg, void *data)
4555 {
4556   CVodeMem cv_mem;
4557   char err_type[10];
4558
4559   /* data points to cv_mem here */
4560
4561   cv_mem = (CVodeMem) data;
4562
4563   if (error_code == CV_WARNING)
4564     sprintf(err_type,"WARNING");
4565   else
4566     sprintf(err_type,"ERROR");
4567
4568 #ifndef NO_FPRINTF_OUTPUT
4569   if (errfp!=NULL) {
4570     fprintf(errfp,"\n[%s %s]  %s\n",module,err_type,function);
4571     fprintf(errfp,"  %s\n\n",msg);
4572   }
4573 #endif
4574
4575   return;
4576 }
4577
4578 /* SUNDIALS EXTENSION */
4579 /* ALL NEXT LINES ADDED FOR EXTENSION */
4580
4581 /*-----------------------------------------------------------------*/
4582
4583 /* 
4584 * CVRcheck1
4585 *
4586 * This routine completes the initialization of rootfinding memory
4587 * information, and checks whether g has a zero both at and very near
4588 * the initial point of the IVP.
4589 *
4590 * This routine returns an int equal to:
4591 *  INITROOT = -1 if a close pair of zeros was found, and
4592 *  CV_SUCCESS     =  0 otherwise.
4593 */
4594
4595 static int CVRcheck1Ext(CVodeMem cv_mem)
4596 {
4597         int i, retval;
4598         booleantype zroot;
4599
4600         for (i = 0; i < nrtfn; i++) iroots[i] = 0;
4601         tlo = tn;
4602         ttol = (ABS(tn) + ABS(h))*uround*HUN;
4603
4604         /* Evaluate g at initial t and check for zero values. */
4605         retval = gfun(tlo, zn[0], glo, g_data);
4606         nge = 1;
4607         if (retval != 0) return(CV_RTFUNC_FAIL);
4608
4609         zroot = FALSE;
4610         for (i = 0; i < nrtfn; i++) {
4611                 if (ABS(glo[i]) == ZERO) 
4612                         iroots[i] =MASKED; /* arbitrary choice*/
4613                 else 
4614                         iroots[i] =0;
4615         }
4616         return(CV_SUCCESS);
4617
4618         /* Some g_i is zero at t0; look at g at t0+(small increment). */
4619
4620
4621
4622 /*
4623 * CVRcheck2
4624 *
4625 * This routine checks for exact zeros of g at the last root found,
4626 * if the last return was a root.  It then checks for a close
4627 * pair of zeros (an error condition), and for a new root at a
4628 * nearby point.  The left endpoint (tlo) of the search interval
4629 * is adjusted if necessary to assure that all g_i are nonzero
4630 * there, before returning to do a root search in the interval.
4631 *
4632 * On entry, tlo = tretlast is the last value of tret returned by
4633 * CVode.  This may be the previous tn, the previous tout value, or
4634 * the last root location.
4635 *
4636 * This routine returns an int equal to:
4637 *      CLOSERT = -2 if a close pair of zeros was found,
4638 *      RTFOUND =  1 if a new zero of g was found near tlo, or
4639 *      CV_SUCCESS    =  0 otherwise.
4640 */
4641
4642 static int CVRcheck2Ext(CVodeMem cv_mem)
4643 {
4644         int i, retval;
4645
4646         if (irfnd == 0) return(CV_SUCCESS);
4647
4648         (void) CVodeGetDky(cv_mem, tlo, 0, y);
4649         retval = gfun(tlo, y, glo, g_data);
4650         nge++;
4651         if (retval != 0) return(CV_RTFUNC_FAIL);
4652
4653         for (i = 0; i < nrtfn; i++) {
4654                 if (ABS(glo[i]) == ZERO) 
4655                         iroots[i] =MASKED; /* arbitrary choice*/
4656                 else 
4657                         iroots[i] =0;
4658         }
4659         return(CV_SUCCESS);
4660
4661 }
4662
4663 /*
4664 * CVRcheck3
4665 *
4666 * This routine interfaces to CVRootfind to look for a root of g
4667 * between tlo and either tn or tout, whichever comes first.
4668 * Only roots beyond tlo in the direction of integration are sought.
4669 *
4670 * This routine returns an int equal to:
4671 *      RTFOUND =  1 if a root of g was found, or
4672 *      CV_SUCCESS    =  0 otherwise.
4673 */
4674
4675 static int CVRcheck3Ext(CVodeMem cv_mem)
4676 {
4677         int i, retval, ier;
4678         /* Set thi = tn or tout, whichever comes first; set y = y(thi). */
4679         if (taskc == CV_ONE_STEP) {
4680                 thi = tn;
4681                 N_VScale(ONE, zn[0], y);
4682         }
4683         if (taskc == CV_NORMAL) {
4684                 if ( (toutc - tn)*h >= ZERO) {
4685                         thi = tn; 
4686                         N_VScale(ONE, zn[0], y);
4687                 } else {
4688                         thi = toutc;
4689                         (void) CVodeGetDky(cv_mem, thi, 0, y);
4690                 }
4691         }
4692
4693         /* Set ghi = g(thi) and call CVRootfind to search (tlo,thi) for roots. */
4694         retval = gfun(thi, y, ghi, g_data);
4695         nge++;
4696         if (retval != 0) return(CV_RTFUNC_FAIL);
4697
4698         ttol = (ABS(tn) + ABS(h))*uround*HUN;
4699         ier = CVRootfind(cv_mem);
4700         tlo = trout;
4701         for (i = 0; i < nrtfn; i++) glo[i] = grout[i];
4702
4703         /* If no root found, return CV_SUCCESS. */  
4704         if (ier == CV_SUCCESS) return(CV_SUCCESS);
4705
4706         /* If a root was found, interpolate to get y(trout) and return.  */
4707         (void) CVodeGetDky(cv_mem, trout, 0, y);
4708         // return(RTFOUND);
4709
4710         if (ier == RTFOUND)
4711                 return(RTFOUND);
4712         else
4713                 return(ZERODETACHING);
4714 }
4715
4716 /*
4717 * CVRootFind
4718 *
4719 * This routine solves for a root of g(t) between tlo and thi, if
4720 * one exists.  Only roots of odd multiplicity (i.e. with a change
4721 * of sign in one of the g_i), or exact zeros, are found.
4722 * Here the sign of tlo - thi is arbitrary, but if multiple roots
4723 * are found, the one closest to tlo is returned.
4724 *
4725 * The method used is the Illinois algorithm, a modified secant method.
4726 * Reference: Kathie L. Hiebert and Lawrence F. Shampine, Implicitly
4727 * Defined Output Points for Solutions of ODEs, Sandia National
4728 * Laboratory Report SAND80-0180, February 1980.
4729 *
4730 * This routine uses the following parameters for communication:
4731 *
4732 * nrtfn    = number of functions g_i, or number of components of
4733 *            the vector-valued function g(t).  Input only.
4734 *
4735 * gfun     = user-defined function for g(t).  Its form is
4736 *            (void) gfun(t, y, gt, g_data)
4737 *
4738 * nge      = cumulative counter for gfun calls.
4739 *
4740 * ttol     = a convergence tolerance for trout.  Input only.
4741 *            When a root at trout is found, it is located only to
4742 *            within a tolerance of ttol.  Typically, ttol should
4743 *            be set to a value on the order of
4744 *               100 * UROUND * max (ABS(tlo), ABS(thi))
4745 *            where UROUND is the unit roundoff of the machine.
4746 *
4747 * tlo, thi = endpoints of the interval in which roots are sought.
4748 *            On input, and must be distinct, but tlo - thi may
4749 *            be of either sign.  The direction of integration is
4750 *            assumed to be from tlo to thi.  On return, tlo and thi
4751 *            are the endpoints of the final relevant interval.
4752 *
4753 * glo, ghi = arrays of length nrtfn containing the vectors g(tlo)
4754 *            and g(thi) respectively.  Input and output.  On input,
4755 *            none of the glo[i] should be zero.
4756 *
4757 * trout    = root location, if a root was found, or thi if not.
4758 *            Output only.  If a root was found other than an exact
4759 *            zero of g, trout is the endpoint thi of the final
4760 *            interval bracketing the root, with size at most ttol.
4761 *
4762 * grout    = array of length nrtfn containing g(trout) on return.
4763 *
4764 * iroots   = int array of length nrtfn with root information.
4765 *            Output only.  If a root was found, iroots indicates
4766 *            which components g_i have a root at trout.  For
4767 *            i = 0, ..., nrtfn-1, iroots[i] = 1 if g_i has a root
4768 *            and iroots[i] = 0 otherwise.
4769 *
4770 * This routine returns an int equal to:
4771 *      RTFOUND =  1 if a root of g was found, or
4772 *      CV_SUCCESS    =  0 otherwise.
4773 */
4774
4775 static int CVRootfindExt(CVodeMem cv_mem)
4776 {
4777         realtype alpha, tmid, gfrac, maxfrac, fracint, fracsub;
4778         int i, retval, imax, side, sideprev;
4779         int istuck,iunstuck,imaxold;
4780
4781         booleantype zroot, umroot, sgnchg;
4782
4783         imax = -1;
4784         istuck=-1;
4785         iunstuck=-1;
4786         maxfrac = ZERO;
4787
4788         /* First check for change in sign in ghi or for a zero in ghi. */
4789         zroot = FALSE;
4790
4791         for (i = 0;  i < nrtfn; i++) {
4792                 if ((ABS(ghi[i])==ZERO)&& (iroots[i]!=MASKED))  istuck=i;
4793                 if ((ABS(ghi[i])> ZERO)&& (iroots[i]==MASKED))  iunstuck=i;
4794                 if ((ABS(ghi[i])> ZERO)&& (glo[i]*ghi[i] <= ZERO)) {
4795                         gfrac = ABS(ghi[i]/(ghi[i] - glo[i]));
4796                         if (gfrac > maxfrac) { /* finding the very first root*/
4797                                 maxfrac = gfrac;
4798                                 imax = i;
4799                         }      
4800                 }
4801         }
4802
4803         if (imax>=0)
4804                 sgnchg=TRUE;
4805         else if (istuck>=0) {
4806                 sgnchg=TRUE;
4807                 imax=istuck;
4808         }else  if (iunstuck>=0) {
4809                 sgnchg=TRUE;
4810                 imax=iunstuck;
4811         }else
4812                 sgnchg = FALSE;
4813
4814         if (!sgnchg) {
4815                 trout = thi;
4816                 for (i = 0; i < nrtfn; i++) grout[i] = ghi[i];
4817                 return(CV_SUCCESS);
4818         }
4819
4820         /* Initialize alpha to avoid compiler warning */
4821         alpha = ONE;
4822
4823         /* A sign change was found.  Loop to locate nearest root. */
4824
4825         side = 0;  sideprev = -1;
4826         for(;;) {                                    /* Looping point */
4827
4828                 /* Set weight alpha.
4829                 On the first two passes, set alpha = 1.  Thereafter, reset alpha
4830                 according to the side (low vs high) of the subinterval in which
4831                 the sign change was found in the previous two passes.
4832                 If the sides were opposite, set alpha = 1.
4833                 If the sides were the same, then double alpha (if high side),
4834                 or halve alpha (if low side).
4835                 The next guess tmid is the secant method value if alpha = 1, but
4836                 is closer to tlo if alpha < 1, and closer to thi if alpha > 1.    */
4837
4838                 if (sideprev == side) {
4839                         alpha = (side == 2) ? alpha*TWO : alpha*HALF;
4840                 } else {
4841                         alpha = ONE;
4842                 }
4843                 /* Set next root approximation tmid and get g(tmid).
4844                 If tmid is too close to tlo or thi, adjust it inward,
4845                 by a fractional distance that is between 0.1 and 0.5.  */
4846                 if ((ABS(ghi[imax])==ZERO)||(ABS(glo[imax])==ZERO)){
4847                         tmid=(tlo+alpha*thi)/(1+alpha);
4848                 }else{
4849                         tmid = thi - (thi - tlo)*ghi[imax]/(ghi[imax] - alpha*glo[imax]);
4850                 }
4851
4852                 if (tmid+1 ==tmid) {printf("tmid is nan\n\r ");exit(0);};
4853
4854                 if (ABS(tmid - tlo) < HALF*ttol) {
4855                         fracint = ABS(thi - tlo)/ttol;
4856                         fracsub = (fracint > FIVE) ? TENTH : HALF/fracint;
4857                         tmid = tlo + fracsub*(thi - tlo);
4858                 }
4859
4860                 if (ABS(thi - tmid) < HALF*ttol) {
4861                         fracint = ABS(thi - tlo)/ttol;
4862                         fracsub = (fracint > FIVE) ? TENTH : HALF/fracint;
4863                         tmid = thi - fracsub*(thi - tlo);
4864                 }
4865
4866                 (void) CVodeGetDky(cv_mem, tmid, 0, y);
4867                 retval = gfun(tmid, y, grout, g_data);
4868                 nge++;
4869                 if (retval != 0) return(CV_RTFUNC_FAIL);
4870
4871                 /* Check to see in which subinterval g changes sign, and reset imax.
4872                 Set side = 1 if sign change is on low side, or 2 if on high side.  */  
4873
4874                 /* First check for change in sign in ghi or for a zero in ghi. */
4875                 zroot = FALSE;
4876                 sideprev = side;
4877                 imaxold=imax;
4878                 imax = -1;
4879                 istuck=-1;iunstuck=-1;
4880                 maxfrac = ZERO;
4881                 for (i = 0;  i < nrtfn; i++) {
4882                         if ((ABS(grout[i])==ZERO)&& (iroots[i]!=MASKED))  istuck=i;
4883                         if ((ABS(grout[i])> ZERO)&& (iroots[i]==MASKED))  iunstuck=i;
4884                         if ((ABS(grout[i])> ZERO)&& (glo[i]*grout[i] <= ZERO)) {
4885                                 gfrac = ABS(grout[i]/(grout[i] - glo[i]));
4886                                 if (gfrac > maxfrac) { /* finding the very first root*/
4887                                         maxfrac = gfrac;
4888                                         imax = i;
4889                                 }      
4890                         }
4891                 }
4892
4893                 if (imax>=0)
4894                         sgnchg=TRUE;
4895                 else if (istuck>=0) {
4896                         sgnchg=TRUE;
4897                         imax=istuck;
4898                 }else  if (iunstuck>=0) {
4899                         sgnchg=TRUE;
4900                         imax=iunstuck;
4901                 }else{
4902                         sgnchg = FALSE;
4903                         imax=imaxold;
4904                 }
4905
4906                 if (sgnchg) {
4907                         /* Sign change found in (tlo,tmid); replace thi with tmid. */
4908                         thi = tmid;
4909                         for (i = 0; i < nrtfn; i++) ghi[i] = grout[i];
4910                         side = 1;
4911                         /* Stop at root thi if converged; otherwise loop. */
4912                         if (ABS(thi - tlo) <= ttol) break;
4913                         continue;  /* Return to looping point. */
4914                 }
4915
4916                 /* here, either (ABS(thi - tlo) <= ttol) or NO SIGN CHANGE */
4917
4918                 /* No sign change in (tlo,tmid), and no zero at tmid.
4919                 Sign change must be in (tmid,thi).  Replace tlo with tmid. */
4920                 tlo = tmid;
4921                 for (i = 0; i < nrtfn; i++) glo[i] = grout[i];
4922                 side = 2;
4923                 /* Stop at root thi if converged; otherwise loop back. */
4924                 if (ABS(thi - tlo) <= ttol) break;
4925
4926         } /* End of root-search loop */
4927
4928         /* Reset trout and grout, set iroots, and return RTFOUND. */
4929         zroot = FALSE;
4930         umroot = FALSE;
4931         trout = thi;
4932         for (i = 0; i < nrtfn; i++) {
4933                 grout[i] = ghi[i];
4934                 if (iroots[i]==MASKED){
4935                         if (ABS(ghi[i]) != ZERO){ 
4936                                 iroots[i] = (ghi[i]> ZERO) ? 2 : -2;
4937                                 umroot=TRUE;
4938                         }else{
4939                                 iroots[i]=0;
4940                         }
4941                 }else{
4942                         if (ABS(ghi[i])== ZERO){ 
4943                                 iroots[i] = (glo[i]> ZERO) ? -1 : 1;
4944                                 zroot = TRUE;
4945                         }else{
4946                                 if (glo[i]*ghi[i] < ZERO){
4947                                         iroots[i] = (ghi[i]>glo[i]) ? 1 : -1;
4948                                         zroot = TRUE;
4949                                 }else
4950                                         iroots[i]=0;
4951                         }
4952                 }    
4953         }
4954         if (zroot) {
4955                 for (i = 0; i < nrtfn; i++) {
4956                         if ((iroots[i]==2)|| (iroots[i]==-2))  iroots[i]=0;
4957                 }
4958                 return(RTFOUND);
4959         }
4960         if (umroot) return(ZERODETACHING);
4961         return(CV_SUCCESS);
4962 }
4963
4964
4965 static int CVRcheck1(CVodeMem cv_mem)
4966 {
4967         /* SUNDIALS EXTENSION */
4968         if (is_sundials_with_extension())
4969         {
4970                 return CVRcheck1Ext(cv_mem);
4971         }
4972         else
4973         {
4974                 return CVRcheck1Std(cv_mem);
4975         }
4976 }
4977
4978 static int CVRcheck2(CVodeMem cv_mem)
4979 {
4980         /* SUNDIALS EXTENSION */
4981         if (is_sundials_with_extension())
4982         {
4983                 return CVRcheck2Ext(cv_mem);
4984         }
4985         else
4986         {
4987                 return CVRcheck2Std(cv_mem);
4988         }
4989 }
4990
4991 static int CVRcheck3(CVodeMem cv_mem)
4992 {
4993         /* SUNDIALS EXTENSION */
4994         if (is_sundials_with_extension())
4995         {
4996                 return CVRcheck3Ext(cv_mem);
4997         }
4998         else
4999         {
5000                 return CVRcheck3Std(cv_mem);
5001         }
5002 }
5003
5004 static int CVRootfind(CVodeMem cv_mem)
5005 {
5006         /* SUNDIALS EXTENSION */
5007         if (is_sundials_with_extension())
5008         {
5009                 return CVRootfindExt(cv_mem);
5010         }
5011         else
5012         {
5013                 return CVRootfindStd(cv_mem);
5014         }
5015 }
5016