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