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