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