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