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