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