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