1 /*
2 * -----------------------------------------------------------------
3 * Programmer(s): Radu Serban @ 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 CVODEA adjoint integrator.
16 * -----------------------------------------------------------------
17 */
18
19 /*
20 * =================================================================
21 * IMPORTED HEADER FILES
22 * =================================================================
23 */
24
25 #include <stdio.h>
26 #include <stdlib.h>
27
28 #include "cvodes_impl.h"
29
30 #include <sundials/sundials_math.h>
31 #include <sundials/sundials_types.h>
32
33 /*
34 * =================================================================
35 * CVODEA PRIVATE CONSTANTS
36 * =================================================================
37 */
38
39 #define ZERO RCONST(0.0) /* real 0.0 */
40 #define ONE RCONST(1.0) /* real 1.0 */
41 #define TWO RCONST(2.0) /* real 2.0 */
42 #define HUNDRED RCONST(100.0) /* real 100.0 */
43 #define FUZZ_FACTOR RCONST(1000000.0) /* fuzz factor for IMget */
44
45 /*
46 * =================================================================
47 * PRIVATE FUNCTION PROTOTYPES
48 * =================================================================
49 */
50
51 static CkpntMem CVAckpntInit(CVodeMem cv_mem);
52 static CkpntMem CVAckpntNew(CVodeMem cv_mem);
53 static void CVAckpntDelete(CkpntMem *ck_memPtr);
54
55 static void CVAbckpbDelete(CVodeBMem *cvB_memPtr);
56
57 static int CVAdataStore(CVodeMem cv_mem, CkpntMem ck_mem);
58 static int CVAckpntGet(CVodeMem cv_mem, CkpntMem ck_mem);
59
60 static int CVAfindIndex(CVodeMem cv_mem, realtype t,
61 long int *indx, booleantype *newpoint);
62
63 static booleantype CVAhermiteMalloc(CVodeMem cv_mem);
64 static void CVAhermiteFree(CVodeMem cv_mem);
65 static int CVAhermiteGetY(CVodeMem cv_mem, realtype t, N_Vector y, N_Vector *yS);
66 static int CVAhermiteStorePnt(CVodeMem cv_mem, DtpntMem d);
67
68 static booleantype CVApolynomialMalloc(CVodeMem cv_mem);
69 static void CVApolynomialFree(CVodeMem cv_mem);
70 static int CVApolynomialGetY(CVodeMem cv_mem, realtype t, N_Vector y, N_Vector *yS);
71 static int CVApolynomialStorePnt(CVodeMem cv_mem, DtpntMem d);
72
73 /* Wrappers */
74
75 static int CVArhs(realtype t, N_Vector yB,
76 N_Vector yBdot, void *cvode_mem);
77
78 static int CVArhsQ(realtype t, N_Vector yB,
79 N_Vector qBdot, void *cvode_mem);
80
81 /*
82 * =================================================================
83 * EXPORTED FUNCTIONS IMPLEMENTATION
84 * =================================================================
85 */
86
87 /*
88 * CVodeAdjInit
89 *
90 * This routine initializes ASA and allocates space for the adjoint
91 * memory structure.
92 */
93
CVodeAdjInit(void * cvode_mem,long int steps,int interp)94 int CVodeAdjInit(void *cvode_mem, long int steps, int interp)
95 {
96 CVadjMem ca_mem;
97 CVodeMem cv_mem;
98 long int i, ii;
99
100 /* ---------------
101 * Check arguments
102 * --------------- */
103
104 if (cvode_mem == NULL) {
105 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeAdjInit", MSGCV_NO_MEM);
106 return(CV_MEM_NULL);
107 }
108 cv_mem = (CVodeMem)cvode_mem;
109
110 if (steps <= 0) {
111 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeAdjInit", MSGCV_BAD_STEPS);
112 return(CV_ILL_INPUT);
113 }
114
115 if ( (interp != CV_HERMITE) && (interp != CV_POLYNOMIAL) ) {
116 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeAdjInit", MSGCV_BAD_INTERP);
117 return(CV_ILL_INPUT);
118 }
119
120 /* ----------------------------
121 * Allocate CVODEA memory block
122 * ---------------------------- */
123
124 ca_mem = NULL;
125 ca_mem = (CVadjMem) malloc(sizeof(struct CVadjMemRec));
126 if (ca_mem == NULL) {
127 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeAdjInit", MSGCV_MEM_FAIL);
128 return(CV_MEM_FAIL);
129 }
130
131 /* Attach ca_mem to CVodeMem structure */
132
133 cv_mem->cv_adj_mem = ca_mem;
134
135 /* ------------------------------
136 * Initialization of check points
137 * ------------------------------ */
138
139 /* Set Check Points linked list to NULL */
140 ca_mem->ck_mem = NULL;
141
142 /* Initialize nckpnts to ZERO */
143 ca_mem->ca_nckpnts = 0;
144
145 /* No interpolation data is available */
146 ca_mem->ca_ckpntData = NULL;
147
148 /* ------------------------------------
149 * Initialization of interpolation data
150 * ------------------------------------ */
151
152 /* Interpolation type */
153
154 ca_mem->ca_IMtype = interp;
155
156 /* Number of steps between check points */
157
158 ca_mem->ca_nsteps = steps;
159
160 /* Last index used in CVAfindIndex, initailize to invalid value */
161 ca_mem->ca_ilast = -1;
162
163 /* Allocate space for the array of Data Point structures */
164
165 ca_mem->dt_mem = NULL;
166 ca_mem->dt_mem = (DtpntMem *) malloc((steps+1)*sizeof(struct DtpntMemRec *));
167 if (ca_mem->dt_mem == NULL) {
168 free(ca_mem); ca_mem = NULL;
169 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeAdjInit", MSGCV_MEM_FAIL);
170 return(CV_MEM_FAIL);
171 }
172
173 for (i=0; i<=steps; i++) {
174 ca_mem->dt_mem[i] = NULL;
175 ca_mem->dt_mem[i] = (DtpntMem) malloc(sizeof(struct DtpntMemRec));
176 if (ca_mem->dt_mem[i] == NULL) {
177 for(ii=0; ii<i; ii++) {free(ca_mem->dt_mem[ii]); ca_mem->dt_mem[ii] = NULL;}
178 free(ca_mem->dt_mem); ca_mem->dt_mem = NULL;
179 free(ca_mem); ca_mem = NULL;
180 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeAdjInit", MSGCV_MEM_FAIL);
181 return(CV_MEM_FAIL);
182 }
183 }
184
185 /* Attach functions for the appropriate interpolation module */
186
187 switch(interp) {
188
189 case CV_HERMITE:
190
191 ca_mem->ca_IMmalloc = CVAhermiteMalloc;
192 ca_mem->ca_IMfree = CVAhermiteFree;
193 ca_mem->ca_IMget = CVAhermiteGetY;
194 ca_mem->ca_IMstore = CVAhermiteStorePnt;
195
196 break;
197
198 case CV_POLYNOMIAL:
199
200 ca_mem->ca_IMmalloc = CVApolynomialMalloc;
201 ca_mem->ca_IMfree = CVApolynomialFree;
202 ca_mem->ca_IMget = CVApolynomialGetY;
203 ca_mem->ca_IMstore = CVApolynomialStorePnt;
204
205 break;
206
207 }
208
209 /* The interpolation module has not been initialized yet */
210
211 ca_mem->ca_IMmallocDone = SUNFALSE;
212
213 /* By default we will store but not interpolate sensitivities
214 * - IMstoreSensi will be set in CVodeF to SUNFALSE if FSA is not enabled
215 * or if the user can force this through CVodeSetAdjNoSensi
216 * - IMinterpSensi will be set in CVodeB to SUNTRUE if IMstoreSensi is
217 * SUNTRUE and if at least one backward problem requires sensitivities */
218
219 ca_mem->ca_IMstoreSensi = SUNTRUE;
220 ca_mem->ca_IMinterpSensi = SUNFALSE;
221
222 /* ------------------------------------
223 * Initialize list of backward problems
224 * ------------------------------------ */
225
226 ca_mem->cvB_mem = NULL;
227 ca_mem->ca_bckpbCrt = NULL;
228 ca_mem->ca_nbckpbs = 0;
229
230 /* --------------------------------
231 * CVodeF and CVodeB not called yet
232 * -------------------------------- */
233
234 ca_mem->ca_firstCVodeFcall = SUNTRUE;
235 ca_mem->ca_tstopCVodeFcall = SUNFALSE;
236
237 ca_mem->ca_firstCVodeBcall = SUNTRUE;
238
239 ca_mem->ca_rootret = SUNFALSE;
240
241 /* ---------------------------------------------
242 * ASA initialized and allocated
243 * --------------------------------------------- */
244
245 cv_mem->cv_adj = SUNTRUE;
246 cv_mem->cv_adjMallocDone = SUNTRUE;
247
248 return(CV_SUCCESS);
249 }
250
251 /* CVodeAdjReInit
252 *
253 * This routine reinitializes the CVODEA memory structure assuming that the
254 * the number of steps between check points and the type of interpolation
255 * remain unchanged.
256 * The list of check points (and associated memory) is deleted.
257 * The list of backward problems is kept (however, new backward problems can
258 * be added to this list by calling CVodeCreateB).
259 * The CVODES memory for the forward and backward problems can be reinitialized
260 * separately by calling CVodeReInit and CVodeReInitB, respectively.
261 * NOTE: if a completely new list of backward problems is also needed, then
262 * simply free the adjoint memory (by calling CVodeAdjFree) and reinitialize
263 * ASA with CVodeAdjInit.
264 */
265
CVodeAdjReInit(void * cvode_mem)266 int CVodeAdjReInit(void *cvode_mem)
267 {
268 CVadjMem ca_mem;
269 CVodeMem cv_mem;
270
271 /* Check cvode_mem */
272 if (cvode_mem == NULL) {
273 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeAdjReInit", MSGCV_NO_MEM);
274 return(CV_MEM_NULL);
275 }
276 cv_mem = (CVodeMem) cvode_mem;
277
278 /* Was ASA initialized? */
279 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
280 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeAdjReInit", MSGCV_NO_ADJ);
281 return(CV_NO_ADJ);
282 }
283
284 ca_mem = cv_mem->cv_adj_mem;
285
286 /* Free current list of Check Points */
287
288 while (ca_mem->ck_mem != NULL) CVAckpntDelete(&(ca_mem->ck_mem));
289
290 /* Initialization of check points */
291
292 ca_mem->ck_mem = NULL;
293 ca_mem->ca_nckpnts = 0;
294 ca_mem->ca_ckpntData = NULL;
295
296 /* CVodeF and CVodeB not called yet */
297
298 ca_mem->ca_firstCVodeFcall = SUNTRUE;
299 ca_mem->ca_tstopCVodeFcall = SUNFALSE;
300 ca_mem->ca_firstCVodeBcall = SUNTRUE;
301
302 return(CV_SUCCESS);
303 }
304
305 /*
306 * CVodeAdjFree
307 *
308 * This routine frees the memory allocated by CVodeAdjInit.
309 */
310
CVodeAdjFree(void * cvode_mem)311 void CVodeAdjFree(void *cvode_mem)
312 {
313 CVodeMem cv_mem;
314 CVadjMem ca_mem;
315 long int i;
316
317 if (cvode_mem == NULL) return;
318 cv_mem = (CVodeMem) cvode_mem;
319
320 if (cv_mem->cv_adjMallocDone) {
321
322 ca_mem = cv_mem->cv_adj_mem;
323
324 /* Delete check points one by one */
325 while (ca_mem->ck_mem != NULL) CVAckpntDelete(&(ca_mem->ck_mem));
326
327 /* Free vectors at all data points */
328 if (ca_mem->ca_IMmallocDone) {
329 ca_mem->ca_IMfree(cv_mem);
330 }
331 for(i=0; i<=ca_mem->ca_nsteps; i++) {
332 free(ca_mem->dt_mem[i]);
333 ca_mem->dt_mem[i] = NULL;
334 }
335 free(ca_mem->dt_mem);
336 ca_mem->dt_mem = NULL;
337
338 /* Delete backward problems one by one */
339 while (ca_mem->cvB_mem != NULL) CVAbckpbDelete(&(ca_mem->cvB_mem));
340
341 /* Free CVODEA memory */
342 free(ca_mem);
343 cv_mem->cv_adj_mem = NULL;
344
345 }
346
347 }
348
349 /*
350 * CVodeF
351 *
352 * This routine integrates to tout and returns solution into yout.
353 * In the same time, it stores check point data every 'steps' steps.
354 *
355 * CVodeF can be called repeatedly by the user.
356 *
357 * ncheckPtr points to the number of check points stored so far.
358 */
359
CVodeF(void * cvode_mem,realtype tout,N_Vector yout,realtype * tret,int itask,int * ncheckPtr)360 int CVodeF(void *cvode_mem, realtype tout, N_Vector yout,
361 realtype *tret, int itask, int *ncheckPtr)
362 {
363 CVadjMem ca_mem;
364 CVodeMem cv_mem;
365 CkpntMem tmp;
366 DtpntMem *dt_mem;
367 long int nstloc;
368 int flag, i;
369 booleantype allocOK, earlyret;
370 realtype ttest;
371
372 /* Check if cvode_mem exists */
373 if (cvode_mem == NULL) {
374 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeF", MSGCV_NO_MEM);
375 return(CV_MEM_NULL);
376 }
377 cv_mem = (CVodeMem) cvode_mem;
378
379 /* Was ASA initialized? */
380 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
381 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeF", MSGCV_NO_ADJ);
382 return(CV_NO_ADJ);
383 }
384
385 ca_mem = cv_mem->cv_adj_mem;
386
387 /* Check for yout != NULL */
388 if (yout == NULL) {
389 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeF", MSGCV_YOUT_NULL);
390 return(CV_ILL_INPUT);
391 }
392
393 /* Check for tret != NULL */
394 if (tret == NULL) {
395 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeF", MSGCV_TRET_NULL);
396 return(CV_ILL_INPUT);
397 }
398
399 /* Check for valid itask */
400 if ( (itask != CV_NORMAL) && (itask != CV_ONE_STEP) ) {
401 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeF", MSGCV_BAD_ITASK);
402 return(CV_ILL_INPUT);
403 }
404
405 /* All error checking done */
406
407 dt_mem = ca_mem->dt_mem;
408
409 /* If tstop is enabled, store some info */
410 if (cv_mem->cv_tstopset) {
411 ca_mem->ca_tstopCVodeFcall = SUNTRUE;
412 ca_mem->ca_tstopCVodeF = cv_mem->cv_tstop;
413 }
414
415 /* On the first step:
416 * - set tinitial
417 * - initialize list of check points
418 * - if needed, initialize the interpolation module
419 * - load dt_mem[0]
420 * On subsequent steps, test if taking a new step is necessary.
421 */
422 if ( ca_mem->ca_firstCVodeFcall ) {
423
424 ca_mem->ca_tinitial = cv_mem->cv_tn;
425
426 ca_mem->ck_mem = CVAckpntInit(cv_mem);
427 if (ca_mem->ck_mem == NULL) {
428 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeF", MSGCV_MEM_FAIL);
429 return(CV_MEM_FAIL);
430 }
431
432 if ( !ca_mem->ca_IMmallocDone ) {
433
434 /* Do we need to store sensitivities? */
435 if (!cv_mem->cv_sensi) ca_mem->ca_IMstoreSensi = SUNFALSE;
436
437 /* Allocate space for interpolation data */
438 allocOK = ca_mem->ca_IMmalloc(cv_mem);
439 if (!allocOK) {
440 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeF", MSGCV_MEM_FAIL);
441 return(CV_MEM_FAIL);
442 }
443
444 /* Rename zn and, if needed, znS for use in interpolation */
445 for (i=0;i<L_MAX;i++) ca_mem->ca_Y[i] = cv_mem->cv_zn[i];
446 if (ca_mem->ca_IMstoreSensi) {
447 for (i=0;i<L_MAX;i++) ca_mem->ca_YS[i] = cv_mem->cv_znS[i];
448 }
449
450 ca_mem->ca_IMmallocDone = SUNTRUE;
451
452 }
453
454 dt_mem[0]->t = ca_mem->ck_mem->ck_t0;
455 ca_mem->ca_IMstore(cv_mem, dt_mem[0]);
456
457 ca_mem->ca_firstCVodeFcall = SUNFALSE;
458
459 } else if ( itask == CV_NORMAL ) {
460
461 /* When in normal mode, check if tout was passed or if a previous root was
462 not reported and return an interpolated solution. No changes to ck_mem
463 or dt_mem are needed. */
464
465 /* flag to signal if an early return is needed */
466 earlyret = SUNFALSE;
467
468 /* if a root needs to be reported compare tout to troot otherwise compare
469 to the current time tn */
470 ttest = (ca_mem->ca_rootret) ? ca_mem->ca_troot : cv_mem->cv_tn;
471
472 if ((ttest - tout)*cv_mem->cv_h >= ZERO) {
473 /* ttest is after tout, interpolate to tout */
474 *tret = tout;
475 flag = CVodeGetDky(cv_mem, tout, 0, yout);
476 earlyret = SUNTRUE;
477 } else if (ca_mem->ca_rootret) {
478 /* tout is after troot, interpolate to troot */
479 *tret = ca_mem->ca_troot;
480 flag = CVodeGetDky(cv_mem, ca_mem->ca_troot, 0, yout);
481 flag = CV_ROOT_RETURN;
482 ca_mem->ca_rootret = SUNFALSE;
483 earlyret = SUNTRUE;
484 }
485
486 /* return if necessary */
487 if (earlyret) {
488 *ncheckPtr = ca_mem->ca_nckpnts;
489 ca_mem->ca_IMnewData = SUNTRUE;
490 ca_mem->ca_ckpntData = ca_mem->ck_mem;
491 ca_mem->ca_np = cv_mem->cv_nst % ca_mem->ca_nsteps + 1;
492 return(flag);
493 }
494
495 }
496
497 /* Integrate to tout (in CV_ONE_STEP mode) while loading check points */
498 nstloc = 0;
499 for(;;) {
500
501 /* Check for too many steps */
502
503 if ( (cv_mem->cv_mxstep>0) && (nstloc >= cv_mem->cv_mxstep) ) {
504 cvProcessError(cv_mem, CV_TOO_MUCH_WORK, "CVODEA", "CVodeF",
505 MSGCV_MAX_STEPS, cv_mem->cv_tn);
506 flag = CV_TOO_MUCH_WORK;
507 break;
508 }
509
510 /* Perform one step of the integration */
511
512 flag = CVode(cv_mem, tout, yout, tret, CV_ONE_STEP);
513 if (flag < 0) break;
514
515 nstloc++;
516
517 /* Test if a new check point is needed */
518
519 if ( cv_mem->cv_nst % ca_mem->ca_nsteps == 0 ) {
520
521 ca_mem->ck_mem->ck_t1 = cv_mem->cv_tn;
522
523 /* Create a new check point, load it, and append it to the list */
524 tmp = CVAckpntNew(cv_mem);
525 if (tmp == NULL) {
526 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeF", MSGCV_MEM_FAIL);
527 flag = CV_MEM_FAIL;
528 break;
529 }
530 tmp->ck_next = ca_mem->ck_mem;
531 ca_mem->ck_mem = tmp;
532 ca_mem->ca_nckpnts++;
533 cv_mem->cv_forceSetup = SUNTRUE;
534
535 /* Reset i=0 and load dt_mem[0] */
536 dt_mem[0]->t = ca_mem->ck_mem->ck_t0;
537 ca_mem->ca_IMstore(cv_mem, dt_mem[0]);
538
539 } else {
540
541 /* Load next point in dt_mem */
542 dt_mem[cv_mem->cv_nst % ca_mem->ca_nsteps]->t = cv_mem->cv_tn;
543 ca_mem->ca_IMstore(cv_mem, dt_mem[cv_mem->cv_nst % ca_mem->ca_nsteps]);
544
545 }
546
547 /* Set t1 field of the current ckeck point structure
548 for the case in which there will be no future
549 check points */
550 ca_mem->ck_mem->ck_t1 = cv_mem->cv_tn;
551
552 /* tfinal is now set to tn */
553 ca_mem->ca_tfinal = cv_mem->cv_tn;
554
555 /* Return if in CV_ONE_STEP mode */
556 if (itask == CV_ONE_STEP) break;
557
558 /* CV_NORMAL_STEP returns */
559
560 /* Return if tout reached */
561 if ( (*tret - tout)*cv_mem->cv_h >= ZERO ) {
562
563 /* If this was a root return, save the root time to return later */
564 if (flag == CV_ROOT_RETURN) {
565 ca_mem->ca_rootret = SUNTRUE;
566 ca_mem->ca_troot = *tret;
567 }
568
569 /* Get solution value at tout to return now */
570 *tret = tout;
571 flag = CVodeGetDky(cv_mem, tout, 0, yout);
572
573 /* Reset tretlast in cv_mem so that CVodeGetQuad and CVodeGetSens
574 * evaluate quadratures and/or sensitivities at the proper time */
575 cv_mem->cv_tretlast = tout;
576
577 break;
578 }
579
580 /* Return if tstop or a root was found */
581 if ( (flag == CV_TSTOP_RETURN) || (flag == CV_ROOT_RETURN) ) break;
582
583 } /* end of for(;;) */
584
585 /* Get ncheck from ca_mem */
586 *ncheckPtr = ca_mem->ca_nckpnts;
587
588 /* Data is available for the last interval */
589 ca_mem->ca_IMnewData = SUNTRUE;
590 ca_mem->ca_ckpntData = ca_mem->ck_mem;
591 ca_mem->ca_np = cv_mem->cv_nst % ca_mem->ca_nsteps + 1;
592
593 return(flag);
594 }
595
596
597
598 /*
599 * =================================================================
600 * FUNCTIONS FOR BACKWARD PROBLEMS
601 * =================================================================
602 */
603
604
CVodeCreateB(void * cvode_mem,int lmmB,int * which)605 int CVodeCreateB(void *cvode_mem, int lmmB, int *which)
606 {
607 CVodeMem cv_mem;
608 CVadjMem ca_mem;
609 CVodeBMem new_cvB_mem;
610 void *cvodeB_mem;
611
612 /* Check if cvode_mem exists */
613 if (cvode_mem == NULL) {
614 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeCreateB", MSGCV_NO_MEM);
615 return(CV_MEM_NULL);
616 }
617 cv_mem = (CVodeMem) cvode_mem;
618
619 /* Was ASA initialized? */
620 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
621 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeCreateB", MSGCV_NO_ADJ);
622 return(CV_NO_ADJ);
623 }
624 ca_mem = cv_mem->cv_adj_mem;
625
626 /* Allocate space for new CVodeBMem object */
627
628 new_cvB_mem = NULL;
629 new_cvB_mem = (CVodeBMem) malloc(sizeof(struct CVodeBMemRec));
630 if (new_cvB_mem == NULL) {
631 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeCreateB", MSGCV_MEM_FAIL);
632 return(CV_MEM_FAIL);
633 }
634
635 /* Create and set a new CVODES object for the backward problem */
636
637 cvodeB_mem = CVodeCreate(lmmB);
638 if (cvodeB_mem == NULL) {
639 cvProcessError(cv_mem, CV_MEM_FAIL, "CVODEA", "CVodeCreateB", MSGCV_MEM_FAIL);
640 return(CV_MEM_FAIL);
641 }
642
643 CVodeSetUserData(cvodeB_mem, cvode_mem);
644
645 CVodeSetMaxHnilWarns(cvodeB_mem, -1);
646
647 CVodeSetErrHandlerFn(cvodeB_mem, cv_mem->cv_ehfun, cv_mem->cv_eh_data);
648 CVodeSetErrFile(cvodeB_mem, cv_mem->cv_errfp);
649
650 /* Set/initialize fields in the new CVodeBMem object, new_cvB_mem */
651
652 new_cvB_mem->cv_index = ca_mem->ca_nbckpbs;
653
654 new_cvB_mem->cv_mem = (CVodeMem) cvodeB_mem;
655
656 new_cvB_mem->cv_f = NULL;
657 new_cvB_mem->cv_fs = NULL;
658
659 new_cvB_mem->cv_fQ = NULL;
660 new_cvB_mem->cv_fQs = NULL;
661
662 new_cvB_mem->cv_user_data = NULL;
663
664 new_cvB_mem->cv_lmem = NULL;
665 new_cvB_mem->cv_lfree = NULL;
666 new_cvB_mem->cv_pmem = NULL;
667 new_cvB_mem->cv_pfree = NULL;
668
669 new_cvB_mem->cv_y = NULL;
670
671 new_cvB_mem->cv_f_withSensi = SUNFALSE;
672 new_cvB_mem->cv_fQ_withSensi = SUNFALSE;
673
674 /* Attach the new object to the linked list cvB_mem */
675
676 new_cvB_mem->cv_next = ca_mem->cvB_mem;
677 ca_mem->cvB_mem = new_cvB_mem;
678
679 /* Return the index of the newly created CVodeBMem object.
680 * This must be passed to CVodeInitB and to other ***B
681 * functions to set optional inputs for this backward problem */
682
683 *which = ca_mem->ca_nbckpbs;
684
685 ca_mem->ca_nbckpbs++;
686
687 return(CV_SUCCESS);
688 }
689
CVodeInitB(void * cvode_mem,int which,CVRhsFnB fB,realtype tB0,N_Vector yB0)690 int CVodeInitB(void *cvode_mem, int which,
691 CVRhsFnB fB,
692 realtype tB0, N_Vector yB0)
693 {
694 CVodeMem cv_mem;
695 CVadjMem ca_mem;
696 CVodeBMem cvB_mem;
697 void *cvodeB_mem;
698 int flag;
699
700 /* Check if cvode_mem exists */
701
702 if (cvode_mem == NULL) {
703 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeInitB", MSGCV_NO_MEM);
704 return(CV_MEM_NULL);
705 }
706 cv_mem = (CVodeMem) cvode_mem;
707
708 /* Was ASA initialized? */
709
710 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
711 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeInitB", MSGCV_NO_ADJ);
712 return(CV_NO_ADJ);
713 }
714 ca_mem = cv_mem->cv_adj_mem;
715
716 /* Check the value of which */
717
718 if ( which >= ca_mem->ca_nbckpbs ) {
719 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeInitB", MSGCV_BAD_WHICH);
720 return(CV_ILL_INPUT);
721 }
722
723 /* Find the CVodeBMem entry in the linked list corresponding to which */
724
725 cvB_mem = ca_mem->cvB_mem;
726 while (cvB_mem != NULL) {
727 if ( which == cvB_mem->cv_index ) break;
728 cvB_mem = cvB_mem->cv_next;
729 }
730
731 cvodeB_mem = (void *) (cvB_mem->cv_mem);
732
733 /* Allocate and set the CVODES object */
734
735 flag = CVodeInit(cvodeB_mem, CVArhs, tB0, yB0);
736
737 if (flag != CV_SUCCESS) return(flag);
738
739 /* Copy fB function in cvB_mem */
740
741 cvB_mem->cv_f_withSensi = SUNFALSE;
742 cvB_mem->cv_f = fB;
743
744 /* Allocate space and initialize the y Nvector in cvB_mem */
745
746 cvB_mem->cv_t0 = tB0;
747 cvB_mem->cv_y = N_VClone(yB0);
748 N_VScale(ONE, yB0, cvB_mem->cv_y);
749
750 return(CV_SUCCESS);
751 }
752
CVodeInitBS(void * cvode_mem,int which,CVRhsFnBS fBs,realtype tB0,N_Vector yB0)753 int CVodeInitBS(void *cvode_mem, int which,
754 CVRhsFnBS fBs,
755 realtype tB0, N_Vector yB0)
756 {
757 CVodeMem cv_mem;
758 CVadjMem ca_mem;
759 CVodeBMem cvB_mem;
760 void *cvodeB_mem;
761 int flag;
762
763 /* Check if cvode_mem exists */
764
765 if (cvode_mem == NULL) {
766 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeInitBS", MSGCV_NO_MEM);
767 return(CV_MEM_NULL);
768 }
769 cv_mem = (CVodeMem) cvode_mem;
770
771 /* Was ASA initialized? */
772
773 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
774 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeInitBS", MSGCV_NO_ADJ);
775 return(CV_NO_ADJ);
776 }
777 ca_mem = cv_mem->cv_adj_mem;
778
779 /* Check the value of which */
780
781 if ( which >= ca_mem->ca_nbckpbs ) {
782 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeInitBS", MSGCV_BAD_WHICH);
783 return(CV_ILL_INPUT);
784 }
785
786 /* Find the CVodeBMem entry in the linked list corresponding to which */
787
788 cvB_mem = ca_mem->cvB_mem;
789 while (cvB_mem != NULL) {
790 if ( which == cvB_mem->cv_index ) break;
791 cvB_mem = cvB_mem->cv_next;
792 }
793
794 cvodeB_mem = (void *) (cvB_mem->cv_mem);
795
796 /* Allocate and set the CVODES object */
797
798 flag = CVodeInit(cvodeB_mem, CVArhs, tB0, yB0);
799
800 if (flag != CV_SUCCESS) return(flag);
801
802 /* Copy fBs function in cvB_mem */
803
804 cvB_mem->cv_f_withSensi = SUNTRUE;
805 cvB_mem->cv_fs = fBs;
806
807 /* Allocate space and initialize the y Nvector in cvB_mem */
808
809 cvB_mem->cv_t0 = tB0;
810 cvB_mem->cv_y = N_VClone(yB0);
811 N_VScale(ONE, yB0, cvB_mem->cv_y);
812
813 return(CV_SUCCESS);
814 }
815
816
CVodeReInitB(void * cvode_mem,int which,realtype tB0,N_Vector yB0)817 int CVodeReInitB(void *cvode_mem, int which,
818 realtype tB0, N_Vector yB0)
819 {
820 CVodeMem cv_mem;
821 CVadjMem ca_mem;
822 CVodeBMem cvB_mem;
823 void *cvodeB_mem;
824 int flag;
825
826 /* Check if cvode_mem exists */
827 if (cvode_mem == NULL) {
828 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeReInitB", MSGCV_NO_MEM);
829 return(CV_MEM_NULL);
830 }
831 cv_mem = (CVodeMem) cvode_mem;
832
833 /* Was ASA initialized? */
834 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
835 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeReInitB", MSGCV_NO_ADJ);
836 return(CV_NO_ADJ);
837 }
838 ca_mem = cv_mem->cv_adj_mem;
839
840 /* Check the value of which */
841 if ( which >= ca_mem->ca_nbckpbs ) {
842 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeReInitB", MSGCV_BAD_WHICH);
843 return(CV_ILL_INPUT);
844 }
845
846 /* Find the CVodeBMem entry in the linked list corresponding to which */
847 cvB_mem = ca_mem->cvB_mem;
848 while (cvB_mem != NULL) {
849 if ( which == cvB_mem->cv_index ) break;
850 cvB_mem = cvB_mem->cv_next;
851 }
852
853 cvodeB_mem = (void *) (cvB_mem->cv_mem);
854
855 /* Reinitialize CVODES object */
856
857 flag = CVodeReInit(cvodeB_mem, tB0, yB0);
858
859 return(flag);
860 }
861
862
CVodeSStolerancesB(void * cvode_mem,int which,realtype reltolB,realtype abstolB)863 int CVodeSStolerancesB(void *cvode_mem, int which, realtype reltolB, realtype abstolB)
864 {
865 CVodeMem cv_mem;
866 CVadjMem ca_mem;
867 CVodeBMem cvB_mem;
868 void *cvodeB_mem;
869 int flag;
870
871 /* Check if cvode_mem exists */
872
873 if (cvode_mem == NULL) {
874 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeSStolerancesB", MSGCV_NO_MEM);
875 return(CV_MEM_NULL);
876 }
877 cv_mem = (CVodeMem) cvode_mem;
878
879 /* Was ASA initialized? */
880
881 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
882 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeSStolerancesB", MSGCV_NO_ADJ);
883 return(CV_NO_ADJ);
884 }
885 ca_mem = cv_mem->cv_adj_mem;
886
887 /* Check the value of which */
888
889 if ( which >= ca_mem->ca_nbckpbs ) {
890 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeSStolerancesB", MSGCV_BAD_WHICH);
891 return(CV_ILL_INPUT);
892 }
893
894 /* Find the CVodeBMem entry in the linked list corresponding to which */
895
896 cvB_mem = ca_mem->cvB_mem;
897 while (cvB_mem != NULL) {
898 if ( which == cvB_mem->cv_index ) break;
899 cvB_mem = cvB_mem->cv_next;
900 }
901
902 cvodeB_mem = (void *) (cvB_mem->cv_mem);
903
904 /* Set tolerances */
905
906 flag = CVodeSStolerances(cvodeB_mem, reltolB, abstolB);
907
908 return(flag);
909 }
910
911
CVodeSVtolerancesB(void * cvode_mem,int which,realtype reltolB,N_Vector abstolB)912 int CVodeSVtolerancesB(void *cvode_mem, int which, realtype reltolB, N_Vector abstolB)
913 {
914 CVodeMem cv_mem;
915 CVadjMem ca_mem;
916 CVodeBMem cvB_mem;
917 void *cvodeB_mem;
918 int flag;
919
920 /* Check if cvode_mem exists */
921
922 if (cvode_mem == NULL) {
923 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeSVtolerancesB", MSGCV_NO_MEM);
924 return(CV_MEM_NULL);
925 }
926 cv_mem = (CVodeMem) cvode_mem;
927
928 /* Was ASA initialized? */
929
930 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
931 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeSVtolerancesB", MSGCV_NO_ADJ);
932 return(CV_NO_ADJ);
933 }
934 ca_mem = cv_mem->cv_adj_mem;
935
936 /* Check the value of which */
937
938 if ( which >= ca_mem->ca_nbckpbs ) {
939 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeSVtolerancesB", MSGCV_BAD_WHICH);
940 return(CV_ILL_INPUT);
941 }
942
943 /* Find the CVodeBMem entry in the linked list corresponding to which */
944
945 cvB_mem = ca_mem->cvB_mem;
946 while (cvB_mem != NULL) {
947 if ( which == cvB_mem->cv_index ) break;
948 cvB_mem = cvB_mem->cv_next;
949 }
950
951 cvodeB_mem = (void *) (cvB_mem->cv_mem);
952
953 /* Set tolerances */
954
955 flag = CVodeSVtolerances(cvodeB_mem, reltolB, abstolB);
956
957 return(flag);
958 }
959
960
CVodeQuadInitB(void * cvode_mem,int which,CVQuadRhsFnB fQB,N_Vector yQB0)961 int CVodeQuadInitB(void *cvode_mem, int which,
962 CVQuadRhsFnB fQB, N_Vector yQB0)
963 {
964 CVodeMem cv_mem;
965 CVadjMem ca_mem;
966 CVodeBMem cvB_mem;
967 void *cvodeB_mem;
968 int flag;
969
970 /* Check if cvode_mem exists */
971 if (cvode_mem == NULL) {
972 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeQuadInitB", MSGCV_NO_MEM);
973 return(CV_MEM_NULL);
974 }
975 cv_mem = (CVodeMem) cvode_mem;
976
977 /* Was ASA initialized? */
978 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
979 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeQuadInitB", MSGCV_NO_ADJ);
980 return(CV_NO_ADJ);
981 }
982 ca_mem = cv_mem->cv_adj_mem;
983
984 /* Check which */
985 if ( which >= ca_mem->ca_nbckpbs ) {
986 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeQuadInitB", MSGCV_BAD_WHICH);
987 return(CV_ILL_INPUT);
988 }
989
990 /* Find the CVodeBMem entry in the linked list corresponding to which */
991 cvB_mem = ca_mem->cvB_mem;
992 while (cvB_mem != NULL) {
993 if ( which == cvB_mem->cv_index ) break;
994 cvB_mem = cvB_mem->cv_next;
995 }
996
997 cvodeB_mem = (void *) (cvB_mem->cv_mem);
998
999 flag = CVodeQuadInit(cvodeB_mem, CVArhsQ, yQB0);
1000 if (flag != CV_SUCCESS) return(flag);
1001
1002 cvB_mem->cv_fQ_withSensi = SUNFALSE;
1003 cvB_mem->cv_fQ = fQB;
1004
1005 return(CV_SUCCESS);
1006 }
1007
CVodeQuadInitBS(void * cvode_mem,int which,CVQuadRhsFnBS fQBs,N_Vector yQB0)1008 int CVodeQuadInitBS(void *cvode_mem, int which,
1009 CVQuadRhsFnBS fQBs, N_Vector yQB0)
1010 {
1011 CVodeMem cv_mem;
1012 CVadjMem ca_mem;
1013 CVodeBMem cvB_mem;
1014 void *cvodeB_mem;
1015 int flag;
1016
1017 /* Check if cvode_mem exists */
1018 if (cvode_mem == NULL) {
1019 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeQuadInitBS", MSGCV_NO_MEM);
1020 return(CV_MEM_NULL);
1021 }
1022 cv_mem = (CVodeMem) cvode_mem;
1023
1024 /* Was ASA initialized? */
1025 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
1026 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeQuadInitBS", MSGCV_NO_ADJ);
1027 return(CV_NO_ADJ);
1028 }
1029 ca_mem = cv_mem->cv_adj_mem;
1030
1031 /* Check which */
1032 if ( which >= ca_mem->ca_nbckpbs ) {
1033 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeQuadInitBS", MSGCV_BAD_WHICH);
1034 return(CV_ILL_INPUT);
1035 }
1036
1037 /* Find the CVodeBMem entry in the linked list corresponding to which */
1038 cvB_mem = ca_mem->cvB_mem;
1039 while (cvB_mem != NULL) {
1040 if ( which == cvB_mem->cv_index ) break;
1041 cvB_mem = cvB_mem->cv_next;
1042 }
1043
1044 cvodeB_mem = (void *) (cvB_mem->cv_mem);
1045
1046 flag = CVodeQuadInit(cvodeB_mem, CVArhsQ, yQB0);
1047 if (flag != CV_SUCCESS) return(flag);
1048
1049 cvB_mem->cv_fQ_withSensi = SUNTRUE;
1050 cvB_mem->cv_fQs = fQBs;
1051
1052 return(CV_SUCCESS);
1053 }
1054
CVodeQuadReInitB(void * cvode_mem,int which,N_Vector yQB0)1055 int CVodeQuadReInitB(void *cvode_mem, int which, N_Vector yQB0)
1056 {
1057 CVodeMem cv_mem;
1058 CVadjMem ca_mem;
1059 CVodeBMem cvB_mem;
1060 void *cvodeB_mem;
1061 int flag;
1062
1063 /* Check if cvode_mem exists */
1064 if (cvode_mem == NULL) {
1065 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeQuadReInitB", MSGCV_NO_MEM);
1066 return(CV_MEM_NULL);
1067 }
1068 cv_mem = (CVodeMem) cvode_mem;
1069
1070 /* Was ASA initialized? */
1071 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
1072 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeQuadReInitB", MSGCV_NO_ADJ);
1073 return(CV_NO_ADJ);
1074 }
1075 ca_mem = cv_mem->cv_adj_mem;
1076
1077 /* Check the value of which */
1078 if ( which >= ca_mem->ca_nbckpbs ) {
1079 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeQuadReInitB", MSGCV_BAD_WHICH);
1080 return(CV_ILL_INPUT);
1081 }
1082
1083 /* Find the CVodeBMem entry in the linked list corresponding to which */
1084 cvB_mem = ca_mem->cvB_mem;
1085 while (cvB_mem != NULL) {
1086 if ( which == cvB_mem->cv_index ) break;
1087 cvB_mem = cvB_mem->cv_next;
1088 }
1089
1090 cvodeB_mem = (void *) (cvB_mem->cv_mem);
1091
1092 flag = CVodeQuadReInit(cvodeB_mem, yQB0);
1093 if (flag != CV_SUCCESS) return(flag);
1094
1095 return(CV_SUCCESS);
1096 }
1097
CVodeQuadSStolerancesB(void * cvode_mem,int which,realtype reltolQB,realtype abstolQB)1098 int CVodeQuadSStolerancesB(void *cvode_mem, int which, realtype reltolQB, realtype abstolQB)
1099 {
1100 CVodeMem cv_mem;
1101 CVadjMem ca_mem;
1102 CVodeBMem cvB_mem;
1103 void *cvodeB_mem;
1104 int flag;
1105
1106 /* Check if cvode_mem exists */
1107 if (cvode_mem == NULL) {
1108 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeQuadSStolerancesB", MSGCV_NO_MEM);
1109 return(CV_MEM_NULL);
1110 }
1111 cv_mem = (CVodeMem) cvode_mem;
1112
1113 /* Was ASA initialized? */
1114 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
1115 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeQuadSStolerancesB", MSGCV_NO_ADJ);
1116 return(CV_NO_ADJ);
1117 }
1118 ca_mem = cv_mem->cv_adj_mem;
1119
1120 /* Check which */
1121 if ( which >= ca_mem->ca_nbckpbs ) {
1122 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeQuadSStolerancesB", MSGCV_BAD_WHICH);
1123 return(CV_ILL_INPUT);
1124 }
1125
1126 /* Find the CVodeBMem entry in the linked list corresponding to which */
1127 cvB_mem = ca_mem->cvB_mem;
1128 while (cvB_mem != NULL) {
1129 if ( which == cvB_mem->cv_index ) break;
1130 cvB_mem = cvB_mem->cv_next;
1131 }
1132
1133 cvodeB_mem = (void *) (cvB_mem->cv_mem);
1134
1135 flag = CVodeQuadSStolerances(cvodeB_mem, reltolQB, abstolQB);
1136
1137 return(flag);
1138 }
1139
CVodeQuadSVtolerancesB(void * cvode_mem,int which,realtype reltolQB,N_Vector abstolQB)1140 int CVodeQuadSVtolerancesB(void *cvode_mem, int which, realtype reltolQB, N_Vector abstolQB)
1141 {
1142 CVodeMem cv_mem;
1143 CVadjMem ca_mem;
1144 CVodeBMem cvB_mem;
1145 void *cvodeB_mem;
1146 int flag;
1147
1148 /* Check if cvode_mem exists */
1149 if (cvode_mem == NULL) {
1150 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeQuadSStolerancesB", MSGCV_NO_MEM);
1151 return(CV_MEM_NULL);
1152 }
1153 cv_mem = (CVodeMem) cvode_mem;
1154
1155 /* Was ASA initialized? */
1156 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
1157 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeQuadSStolerancesB", MSGCV_NO_ADJ);
1158 return(CV_NO_ADJ);
1159 }
1160 ca_mem = cv_mem->cv_adj_mem;
1161
1162 /* Check which */
1163 if ( which >= ca_mem->ca_nbckpbs ) {
1164 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeQuadSStolerancesB", MSGCV_BAD_WHICH);
1165 return(CV_ILL_INPUT);
1166 }
1167
1168 /* Find the CVodeBMem entry in the linked list corresponding to which */
1169 cvB_mem = ca_mem->cvB_mem;
1170 while (cvB_mem != NULL) {
1171 if ( which == cvB_mem->cv_index ) break;
1172 cvB_mem = cvB_mem->cv_next;
1173 }
1174
1175 cvodeB_mem = (void *) (cvB_mem->cv_mem);
1176
1177 flag = CVodeQuadSVtolerances(cvodeB_mem, reltolQB, abstolQB);
1178
1179 return(flag);
1180 }
1181
1182 /*
1183 * CVodeB
1184 *
1185 * This routine performs the backward integration towards tBout
1186 * of all backward problems that were defined.
1187 * When necessary, it performs a forward integration between two
1188 * consecutive check points to update interpolation data.
1189 *
1190 * On a successful return, CVodeB returns CV_SUCCESS.
1191 *
1192 * NOTE that CVodeB DOES NOT return the solution for the backward
1193 * problem(s). Use CVodeGetB to extract the solution at tBret
1194 * for any given backward problem.
1195 *
1196 * If there are multiple backward problems and multiple check points,
1197 * CVodeB may not succeed in getting all problems to take one step
1198 * when called in ONE_STEP mode.
1199 */
1200
CVodeB(void * cvode_mem,realtype tBout,int itaskB)1201 int CVodeB(void *cvode_mem, realtype tBout, int itaskB)
1202 {
1203 CVodeMem cv_mem;
1204 CVadjMem ca_mem;
1205 CVodeBMem cvB_mem, tmp_cvB_mem;
1206 CkpntMem ck_mem;
1207 int sign, flag=0;
1208 realtype tfuzz, tBret, tBn;
1209 booleantype gotCheckpoint, isActive, reachedTBout;
1210
1211 /* Check if cvode_mem exists */
1212
1213 if (cvode_mem == NULL) {
1214 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeB", MSGCV_NO_MEM);
1215 return(CV_MEM_NULL);
1216 }
1217 cv_mem = (CVodeMem) cvode_mem;
1218
1219 /* Was ASA initialized? */
1220
1221 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
1222 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeB", MSGCV_NO_ADJ);
1223 return(CV_NO_ADJ);
1224 }
1225 ca_mem = cv_mem->cv_adj_mem;
1226
1227 /* Check if any backward problem has been defined */
1228
1229 if ( ca_mem->ca_nbckpbs == 0 ) {
1230 cvProcessError(cv_mem, CV_NO_BCK, "CVODEA", "CVodeB", MSGCV_NO_BCK);
1231 return(CV_NO_BCK);
1232 }
1233 cvB_mem = ca_mem->cvB_mem;
1234
1235 /* Check whether CVodeF has been called */
1236
1237 if ( ca_mem->ca_firstCVodeFcall ) {
1238 cvProcessError(cv_mem, CV_NO_FWD, "CVODEA", "CVodeB", MSGCV_NO_FWD);
1239 return(CV_NO_FWD);
1240 }
1241 sign = (ca_mem->ca_tfinal - ca_mem->ca_tinitial > ZERO) ? 1 : -1;
1242
1243 /* If this is the first call, loop over all backward problems and
1244 * - check that tB0 is valid
1245 * - check that tBout is ahead of tB0 in the backward direction
1246 * - check whether we need to interpolate forward sensitivities
1247 */
1248
1249 if ( ca_mem->ca_firstCVodeBcall ) {
1250
1251 tmp_cvB_mem = cvB_mem;
1252
1253 while(tmp_cvB_mem != NULL) {
1254
1255 tBn = tmp_cvB_mem->cv_mem->cv_tn;
1256
1257 if ( (sign*(tBn-ca_mem->ca_tinitial) < ZERO) || (sign*(ca_mem->ca_tfinal-tBn) < ZERO) ) {
1258 cvProcessError(cv_mem, CV_BAD_TB0, "CVODEA", "CVodeB", MSGCV_BAD_TB0,
1259 tmp_cvB_mem->cv_index);
1260 return(CV_BAD_TB0);
1261 }
1262
1263 if (sign*(tBn-tBout) <= ZERO) {
1264 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeB", MSGCV_BAD_TBOUT,
1265 tmp_cvB_mem->cv_index);
1266 return(CV_ILL_INPUT);
1267 }
1268
1269 if ( tmp_cvB_mem->cv_f_withSensi || tmp_cvB_mem->cv_fQ_withSensi )
1270 ca_mem->ca_IMinterpSensi = SUNTRUE;
1271
1272 tmp_cvB_mem = tmp_cvB_mem->cv_next;
1273
1274 }
1275
1276 if ( ca_mem->ca_IMinterpSensi && !ca_mem->ca_IMstoreSensi) {
1277 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeB", MSGCV_BAD_SENSI);
1278 return(CV_ILL_INPUT);
1279 }
1280
1281 ca_mem->ca_firstCVodeBcall = SUNFALSE;
1282 }
1283
1284 /* Check if itaskB is legal */
1285
1286 if ( (itaskB != CV_NORMAL) && (itaskB != CV_ONE_STEP) ) {
1287 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeB", MSGCV_BAD_ITASKB);
1288 return(CV_ILL_INPUT);
1289 }
1290
1291 /* Check if tBout is legal */
1292
1293 if ( (sign*(tBout-ca_mem->ca_tinitial) < ZERO) || (sign*(ca_mem->ca_tfinal-tBout) < ZERO) ) {
1294 tfuzz = HUNDRED*cv_mem->cv_uround*(SUNRabs(ca_mem->ca_tinitial) + SUNRabs(ca_mem->ca_tfinal));
1295 if ( (sign*(tBout-ca_mem->ca_tinitial) < ZERO) && (SUNRabs(tBout-ca_mem->ca_tinitial) < tfuzz) ) {
1296 tBout = ca_mem->ca_tinitial;
1297 } else {
1298 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeB", MSGCV_BAD_TBOUT);
1299 return(CV_ILL_INPUT);
1300 }
1301 }
1302
1303 /* Loop through the check points and stop as soon as a backward
1304 * problem has its tn value behind the current check point's t0_
1305 * value (in the backward direction) */
1306
1307 ck_mem = ca_mem->ck_mem;
1308
1309 gotCheckpoint = SUNFALSE;
1310
1311 for(;;) {
1312
1313 tmp_cvB_mem = cvB_mem;
1314 while(tmp_cvB_mem != NULL) {
1315 tBn = tmp_cvB_mem->cv_mem->cv_tn;
1316
1317 if ( sign*(tBn-ck_mem->ck_t0) > ZERO ) {
1318 gotCheckpoint = SUNTRUE;
1319 break;
1320 }
1321
1322 if ( (itaskB==CV_NORMAL) && (tBn == ck_mem->ck_t0) && (sign*(tBout-ck_mem->ck_t0) >= ZERO) ) {
1323 gotCheckpoint = SUNTRUE;
1324 break;
1325 }
1326
1327 tmp_cvB_mem = tmp_cvB_mem->cv_next;
1328 }
1329
1330 if (gotCheckpoint) break;
1331
1332 if (ck_mem->ck_next == NULL) break;
1333
1334 ck_mem = ck_mem->ck_next;
1335 }
1336
1337 /* Starting with the current check point from above, loop over check points
1338 while propagating backward problems */
1339
1340 for(;;) {
1341
1342 /* Store interpolation data if not available.
1343 This is the 2nd forward integration pass */
1344
1345 if (ck_mem != ca_mem->ca_ckpntData) {
1346 flag = CVAdataStore(cv_mem, ck_mem);
1347 if (flag != CV_SUCCESS) break;
1348 }
1349
1350 /* Loop through all backward problems and, if needed,
1351 * propagate their solution towards tBout */
1352
1353 tmp_cvB_mem = cvB_mem;
1354 while (tmp_cvB_mem != NULL) {
1355
1356 /* Decide if current backward problem is "active" in this check point */
1357
1358 isActive = SUNTRUE;
1359
1360 tBn = tmp_cvB_mem->cv_mem->cv_tn;
1361
1362 if ( (tBn == ck_mem->ck_t0) && (sign*(tBout-ck_mem->ck_t0) < ZERO ) ) isActive = SUNFALSE;
1363 if ( (tBn == ck_mem->ck_t0) && (itaskB==CV_ONE_STEP) ) isActive = SUNFALSE;
1364
1365 if ( sign * (tBn - ck_mem->ck_t0) < ZERO ) isActive = SUNFALSE;
1366
1367 if ( isActive ) {
1368
1369 /* Store the address of current backward problem memory
1370 * in ca_mem to be used in the wrapper functions */
1371 ca_mem->ca_bckpbCrt = tmp_cvB_mem;
1372
1373 /* Integrate current backward problem */
1374 CVodeSetStopTime(tmp_cvB_mem->cv_mem, ck_mem->ck_t0);
1375 flag = CVode(tmp_cvB_mem->cv_mem, tBout, tmp_cvB_mem->cv_y, &tBret, itaskB);
1376
1377 /* Set the time at which we will report solution and/or quadratures */
1378 tmp_cvB_mem->cv_tout = tBret;
1379
1380 /* If an error occurred, exit while loop */
1381 if (flag < 0) break;
1382
1383 } else {
1384 flag = CV_SUCCESS;
1385 tmp_cvB_mem->cv_tout = tBn;
1386 }
1387
1388 /* Move to next backward problem */
1389
1390 tmp_cvB_mem = tmp_cvB_mem->cv_next;
1391 }
1392
1393 /* If an error occurred, return now */
1394
1395 if (flag <0) {
1396 cvProcessError(cv_mem, flag, "CVODEA", "CVodeB", MSGCV_BACK_ERROR,
1397 tmp_cvB_mem->cv_index);
1398 return(flag);
1399 }
1400
1401 /* If in CV_ONE_STEP mode, return now (flag = CV_SUCCESS) */
1402
1403 if (itaskB == CV_ONE_STEP) break;
1404
1405 /* If all backward problems have succesfully reached tBout, return now */
1406
1407 reachedTBout = SUNTRUE;
1408
1409 tmp_cvB_mem = cvB_mem;
1410 while(tmp_cvB_mem != NULL) {
1411 if ( sign*(tmp_cvB_mem->cv_tout - tBout) > ZERO ) {
1412 reachedTBout = SUNFALSE;
1413 break;
1414 }
1415 tmp_cvB_mem = tmp_cvB_mem->cv_next;
1416 }
1417
1418 if ( reachedTBout ) break;
1419
1420 /* Move check point in linked list to next one */
1421
1422 ck_mem = ck_mem->ck_next;
1423
1424 }
1425
1426 return(flag);
1427 }
1428
1429
CVodeGetB(void * cvode_mem,int which,realtype * tret,N_Vector yB)1430 int CVodeGetB(void *cvode_mem, int which, realtype *tret, N_Vector yB)
1431 {
1432 CVodeMem cv_mem;
1433 CVadjMem ca_mem;
1434 CVodeBMem cvB_mem;
1435
1436 /* Check if cvode_mem exists */
1437 if (cvode_mem == NULL) {
1438 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeGetB", MSGCV_NO_MEM);
1439 return(CV_MEM_NULL);
1440 }
1441 cv_mem = (CVodeMem) cvode_mem;
1442
1443 /* Was ASA initialized? */
1444 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
1445 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeGetB", MSGCV_NO_ADJ);
1446 return(CV_NO_ADJ);
1447 }
1448
1449 ca_mem = cv_mem->cv_adj_mem;
1450
1451 /* Check the value of which */
1452 if ( which >= ca_mem->ca_nbckpbs ) {
1453 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeGetB", MSGCV_BAD_WHICH);
1454 return(CV_ILL_INPUT);
1455 }
1456
1457 /* Find the CVodeBMem entry in the linked list corresponding to which */
1458 cvB_mem = ca_mem->cvB_mem;
1459 while (cvB_mem != NULL) {
1460 if ( which == cvB_mem->cv_index ) break;
1461 cvB_mem = cvB_mem->cv_next;
1462 }
1463
1464 N_VScale(ONE, cvB_mem->cv_y, yB);
1465 *tret = cvB_mem->cv_tout;
1466
1467 return(CV_SUCCESS);
1468 }
1469
1470
1471 /*
1472 * CVodeGetQuadB
1473 */
1474
CVodeGetQuadB(void * cvode_mem,int which,realtype * tret,N_Vector qB)1475 int CVodeGetQuadB(void *cvode_mem, int which, realtype *tret, N_Vector qB)
1476 {
1477 CVodeMem cv_mem;
1478 CVadjMem ca_mem;
1479 CVodeBMem cvB_mem;
1480 void *cvodeB_mem;
1481 long int nstB;
1482 int flag;
1483
1484 /* Check if cvode_mem exists */
1485 if (cvode_mem == NULL) {
1486 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeGetQuadB", MSGCV_NO_MEM);
1487 return(CV_MEM_NULL);
1488 }
1489 cv_mem = (CVodeMem) cvode_mem;
1490
1491 /* Was ASA initialized? */
1492 if (cv_mem->cv_adjMallocDone == SUNFALSE) {
1493 cvProcessError(cv_mem, CV_NO_ADJ, "CVODEA", "CVodeGetQuadB", MSGCV_NO_ADJ);
1494 return(CV_NO_ADJ);
1495 }
1496
1497 ca_mem = cv_mem->cv_adj_mem;
1498
1499 /* Check the value of which */
1500 if ( which >= ca_mem->ca_nbckpbs ) {
1501 cvProcessError(cv_mem, CV_ILL_INPUT, "CVODEA", "CVodeGetQuadB", MSGCV_BAD_WHICH);
1502 return(CV_ILL_INPUT);
1503 }
1504
1505 /* Find the CVodeBMem entry in the linked list corresponding to which */
1506 cvB_mem = ca_mem->cvB_mem;
1507 while (cvB_mem != NULL) {
1508 if ( which == cvB_mem->cv_index ) break;
1509 cvB_mem = cvB_mem->cv_next;
1510 }
1511
1512 cvodeB_mem = (void *) (cvB_mem->cv_mem);
1513
1514 /* If the integration for this backward problem has not started yet,
1515 * simply return the current value of qB (i.e. the final conditions) */
1516
1517 flag = CVodeGetNumSteps(cvodeB_mem, &nstB);
1518
1519 if (nstB == 0) {
1520 N_VScale(ONE, cvB_mem->cv_mem->cv_znQ[0], qB);
1521 *tret = cvB_mem->cv_tout;
1522 } else {
1523 flag = CVodeGetQuad(cvodeB_mem, tret, qB);
1524 }
1525
1526 return(flag);
1527 }
1528
1529
1530 /*
1531 * =================================================================
1532 * PRIVATE FUNCTIONS FOR CHECK POINTS
1533 * =================================================================
1534 */
1535
1536 /*
1537 * CVAckpntInit
1538 *
1539 * This routine initializes the check point linked list with
1540 * information from the initial time.
1541 */
1542
CVAckpntInit(CVodeMem cv_mem)1543 static CkpntMem CVAckpntInit(CVodeMem cv_mem)
1544 {
1545 CkpntMem ck_mem;
1546 int is;
1547
1548 /* Allocate space for ckdata */
1549 ck_mem = NULL;
1550 ck_mem = (CkpntMem) malloc(sizeof(struct CkpntMemRec));
1551 if (ck_mem == NULL) return(NULL);
1552
1553 ck_mem->ck_zn[0] = N_VClone(cv_mem->cv_tempv);
1554 if (ck_mem->ck_zn[0] == NULL) {
1555 free(ck_mem); ck_mem = NULL;
1556 return(NULL);
1557 }
1558
1559 ck_mem->ck_zn[1] = N_VClone(cv_mem->cv_tempv);
1560 if (ck_mem->ck_zn[1] == NULL) {
1561 N_VDestroy(ck_mem->ck_zn[0]);
1562 free(ck_mem); ck_mem = NULL;
1563 return(NULL);
1564 }
1565
1566 /* ck_mem->ck_zn[qmax] was not allocated */
1567 ck_mem->ck_zqm = 0;
1568
1569 /* Load ckdata from cv_mem */
1570 N_VScale(ONE, cv_mem->cv_zn[0], ck_mem->ck_zn[0]);
1571 ck_mem->ck_t0 = cv_mem->cv_tn;
1572 ck_mem->ck_nst = 0;
1573 ck_mem->ck_q = 1;
1574 ck_mem->ck_h = 0.0;
1575
1576 /* Do we need to carry quadratures */
1577 ck_mem->ck_quadr = cv_mem->cv_quadr && cv_mem->cv_errconQ;
1578
1579 if (ck_mem->ck_quadr) {
1580
1581 ck_mem->ck_znQ[0] = N_VClone(cv_mem->cv_tempvQ);
1582 if (ck_mem->ck_znQ[0] == NULL) {
1583 N_VDestroy(ck_mem->ck_zn[0]);
1584 N_VDestroy(ck_mem->ck_zn[1]);
1585 free(ck_mem); ck_mem = NULL;
1586 return(NULL);
1587 }
1588
1589 N_VScale(ONE, cv_mem->cv_znQ[0], ck_mem->ck_znQ[0]);
1590
1591 }
1592
1593 /* Do we need to carry sensitivities? */
1594 ck_mem->ck_sensi = cv_mem->cv_sensi;
1595
1596 if (ck_mem->ck_sensi) {
1597
1598 ck_mem->ck_Ns = cv_mem->cv_Ns;
1599
1600 ck_mem->ck_znS[0] = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempv);
1601 if (ck_mem->ck_znS[0] == NULL) {
1602 N_VDestroy(ck_mem->ck_zn[0]);
1603 N_VDestroy(ck_mem->ck_zn[1]);
1604 if (ck_mem->ck_quadr) N_VDestroy(ck_mem->ck_znQ[0]);
1605 free(ck_mem); ck_mem = NULL;
1606 return(NULL);
1607 }
1608
1609 for (is=0; is<cv_mem->cv_Ns; is++)
1610 cv_mem->cv_cvals[is] = ONE;
1611
1612 (void) N_VScaleVectorArray(cv_mem->cv_Ns, cv_mem->cv_cvals,
1613 cv_mem->cv_znS[0], ck_mem->ck_znS[0]);
1614 }
1615
1616 /* Do we need to carry quadrature sensitivities? */
1617 ck_mem->ck_quadr_sensi = cv_mem->cv_quadr_sensi && cv_mem->cv_errconQS;
1618
1619 if (ck_mem->ck_quadr_sensi) {
1620 ck_mem->ck_znQS[0] = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempvQ);
1621 if (ck_mem->ck_znQS[0] == NULL) {
1622 N_VDestroy(ck_mem->ck_zn[0]);
1623 N_VDestroy(ck_mem->ck_zn[1]);
1624 if (ck_mem->ck_quadr) N_VDestroy(ck_mem->ck_znQ[0]);
1625 N_VDestroyVectorArray(ck_mem->ck_znS[0], cv_mem->cv_Ns);
1626 free(ck_mem); ck_mem = NULL;
1627 return(NULL);
1628 }
1629
1630 for (is=0; is<cv_mem->cv_Ns; is++)
1631 cv_mem->cv_cvals[is] = ONE;
1632
1633 (void) N_VScaleVectorArray(cv_mem->cv_Ns, cv_mem->cv_cvals,
1634 cv_mem->cv_znQS[0], ck_mem->ck_znQS[0]);
1635 }
1636
1637 /* Next in list */
1638 ck_mem->ck_next = NULL;
1639
1640 return(ck_mem);
1641 }
1642
1643 /*
1644 * CVAckpntNew
1645 *
1646 * This routine allocates space for a new check point and sets
1647 * its data from current values in cv_mem.
1648 */
1649
CVAckpntNew(CVodeMem cv_mem)1650 static CkpntMem CVAckpntNew(CVodeMem cv_mem)
1651 {
1652 CkpntMem ck_mem;
1653 int j, jj, is, qmax;
1654
1655 /* Allocate space for ckdata */
1656 ck_mem = NULL;
1657 ck_mem = (CkpntMem) malloc(sizeof(struct CkpntMemRec));
1658 if (ck_mem == NULL) return(NULL);
1659
1660 /* Set cv_next to NULL */
1661 ck_mem->ck_next = NULL;
1662
1663 /* Test if we need to allocate space for the last zn.
1664 * NOTE: zn(qmax) may be needed for a hot restart, if an order
1665 * increase is deemed necessary at the first step after a check point */
1666 qmax = cv_mem->cv_qmax;
1667 ck_mem->ck_zqm = (cv_mem->cv_q < qmax) ? qmax : 0;
1668
1669 for (j=0; j<=cv_mem->cv_q; j++) {
1670 ck_mem->ck_zn[j] = N_VClone(cv_mem->cv_tempv);
1671 if (ck_mem->ck_zn[j] == NULL) {
1672 for (jj=0; jj<j; jj++) N_VDestroy(ck_mem->ck_zn[jj]);
1673 free(ck_mem); ck_mem = NULL;
1674 return(NULL);
1675 }
1676 }
1677
1678 if (cv_mem->cv_q < qmax) {
1679 ck_mem->ck_zn[qmax] = N_VClone(cv_mem->cv_tempv);
1680 if (ck_mem->ck_zn[qmax] == NULL) {
1681 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_zn[jj]);
1682 free(ck_mem); ck_mem = NULL;
1683 return(NULL);
1684 }
1685 }
1686
1687 /* Test if we need to carry quadratures */
1688 ck_mem->ck_quadr = cv_mem->cv_quadr && cv_mem->cv_errconQ;
1689
1690 if (ck_mem->ck_quadr) {
1691
1692 for (j=0; j<=cv_mem->cv_q; j++) {
1693 ck_mem->ck_znQ[j] = N_VClone(cv_mem->cv_tempvQ);
1694 if(ck_mem->ck_znQ[j] == NULL) {
1695 for (jj=0; jj<j; jj++) N_VDestroy(ck_mem->ck_znQ[jj]);
1696 if (cv_mem->cv_q < qmax) N_VDestroy(ck_mem->ck_zn[qmax]);
1697 for (jj=0; jj<=cv_mem->cv_q; j++) N_VDestroy(ck_mem->ck_zn[jj]);
1698 free(ck_mem); ck_mem = NULL;
1699 return(NULL);
1700 }
1701 }
1702
1703 if (cv_mem->cv_q < qmax) {
1704 ck_mem->ck_znQ[qmax] = N_VClone(cv_mem->cv_tempvQ);
1705 if (ck_mem->ck_znQ[qmax] == NULL) {
1706 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_znQ[jj]);
1707 N_VDestroy(ck_mem->ck_zn[qmax]);
1708 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_zn[jj]);
1709 free(ck_mem); ck_mem = NULL;
1710 return(NULL);
1711 }
1712 }
1713
1714 }
1715
1716 /* Test if we need to carry sensitivities */
1717 ck_mem->ck_sensi = cv_mem->cv_sensi;
1718
1719 if (ck_mem->ck_sensi) {
1720
1721 ck_mem->ck_Ns = cv_mem->cv_Ns;
1722
1723 for (j=0; j<=cv_mem->cv_q; j++) {
1724 ck_mem->ck_znS[j] = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempv);
1725 if (ck_mem->ck_znS[j] == NULL) {
1726 for (jj=0; jj<j; jj++) N_VDestroyVectorArray(ck_mem->ck_znS[jj], cv_mem->cv_Ns);
1727 if (ck_mem->ck_quadr) {
1728 if (cv_mem->cv_q < qmax) N_VDestroy(ck_mem->ck_znQ[qmax]);
1729 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_znQ[jj]);
1730 }
1731 if (cv_mem->cv_q < qmax) N_VDestroy(ck_mem->ck_zn[qmax]);
1732 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_zn[jj]);
1733 free(ck_mem); ck_mem = NULL;
1734 return(NULL);
1735 }
1736 }
1737
1738 if ( cv_mem->cv_q < qmax) {
1739 ck_mem->ck_znS[qmax] = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempv);
1740 if (ck_mem->ck_znS[qmax] == NULL) {
1741 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroyVectorArray(ck_mem->ck_znS[jj], cv_mem->cv_Ns);
1742 if (ck_mem->ck_quadr) {
1743 N_VDestroy(ck_mem->ck_znQ[qmax]);
1744 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_znQ[jj]);
1745 }
1746 N_VDestroy(ck_mem->ck_zn[qmax]);
1747 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_zn[jj]);
1748 free(ck_mem); ck_mem = NULL;
1749 return(NULL);
1750 }
1751 }
1752
1753 }
1754
1755 /* Test if we need to carry quadrature sensitivities */
1756 ck_mem->ck_quadr_sensi = cv_mem->cv_quadr_sensi && cv_mem->cv_errconQS;
1757
1758 if (ck_mem->ck_quadr_sensi) {
1759
1760 for (j=0; j<=cv_mem->cv_q; j++) {
1761 ck_mem->ck_znQS[j] = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempvQ);
1762 if (ck_mem->ck_znQS[j] == NULL) {
1763 for (jj=0; jj<j; jj++) N_VDestroyVectorArray(ck_mem->ck_znQS[jj], cv_mem->cv_Ns);
1764 if (cv_mem->cv_q < qmax) N_VDestroyVectorArray(ck_mem->ck_znS[qmax], cv_mem->cv_Ns);
1765 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroyVectorArray(ck_mem->ck_znS[jj], cv_mem->cv_Ns);
1766 if (ck_mem->ck_quadr) {
1767 if (cv_mem->cv_q < qmax) N_VDestroy(ck_mem->ck_znQ[qmax]);
1768 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_znQ[jj]);
1769 }
1770 if (cv_mem->cv_q < qmax) N_VDestroy(ck_mem->ck_zn[qmax]);
1771 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_zn[jj]);
1772 free(ck_mem); ck_mem = NULL;
1773 return(NULL);
1774 }
1775 }
1776
1777 if ( cv_mem->cv_q < qmax) {
1778 ck_mem->ck_znQS[qmax] = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempvQ);
1779 if (ck_mem->ck_znQS[qmax] == NULL) {
1780 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroyVectorArray(ck_mem->ck_znQS[jj], cv_mem->cv_Ns);
1781 N_VDestroyVectorArray(ck_mem->ck_znS[qmax], cv_mem->cv_Ns);
1782 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroyVectorArray(ck_mem->ck_znS[jj], cv_mem->cv_Ns);
1783 if (ck_mem->ck_quadr) {
1784 N_VDestroy(ck_mem->ck_znQ[qmax]);
1785 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_zn[jj]);
1786 }
1787 N_VDestroy(ck_mem->ck_zn[qmax]);
1788 for (jj=0; jj<=cv_mem->cv_q; jj++) N_VDestroy(ck_mem->ck_zn[jj]);
1789 free(ck_mem); ck_mem = NULL;
1790 return(NULL);
1791 }
1792 }
1793
1794 }
1795
1796 /* Load check point data from cv_mem */
1797
1798 for (j=0; j<=cv_mem->cv_q; j++)
1799 cv_mem->cv_cvals[j] = ONE;
1800
1801 (void) N_VScaleVectorArray(cv_mem->cv_q+1, cv_mem->cv_cvals,
1802 cv_mem->cv_zn, ck_mem->ck_zn);
1803
1804 if ( cv_mem->cv_q < qmax )
1805 N_VScale(ONE, cv_mem->cv_zn[qmax], ck_mem->ck_zn[qmax]);
1806
1807 if (ck_mem->ck_quadr) {
1808 for (j=0; j<=cv_mem->cv_q; j++)
1809 cv_mem->cv_cvals[j] = ONE;
1810
1811 (void) N_VScaleVectorArray(cv_mem->cv_q+1, cv_mem->cv_cvals,
1812 cv_mem->cv_znQ, ck_mem->ck_znQ);
1813
1814 if ( cv_mem->cv_q < qmax )
1815 N_VScale(ONE, cv_mem->cv_znQ[qmax], ck_mem->ck_znQ[qmax]);
1816 }
1817
1818 if (ck_mem->ck_sensi) {
1819 for (j=0; j<=cv_mem->cv_q; j++) {
1820 for (is=0; is<cv_mem->cv_Ns; is++) {
1821 cv_mem->cv_cvals[j*cv_mem->cv_Ns+is] = ONE;
1822 cv_mem->cv_Xvecs[j*cv_mem->cv_Ns+is] = cv_mem->cv_znS[j][is];
1823 cv_mem->cv_Zvecs[j*cv_mem->cv_Ns+is] = ck_mem->ck_znS[j][is];
1824 }
1825 }
1826
1827 (void) N_VScaleVectorArray(cv_mem->cv_Ns*(cv_mem->cv_q+1),
1828 cv_mem->cv_cvals,
1829 cv_mem->cv_Xvecs, cv_mem->cv_Zvecs);
1830
1831 if ( cv_mem->cv_q < qmax ) {
1832 for (is=0; is<cv_mem->cv_Ns; is++)
1833 cv_mem->cv_cvals[is] = ONE;
1834
1835 (void) N_VScaleVectorArray(cv_mem->cv_Ns, cv_mem->cv_cvals,
1836 cv_mem->cv_znS[qmax], ck_mem->ck_znS[qmax]);
1837 }
1838 }
1839
1840 if (ck_mem->ck_quadr_sensi) {
1841 for (j=0; j<=cv_mem->cv_q; j++) {
1842 for (is=0; is<cv_mem->cv_Ns; is++) {
1843 cv_mem->cv_cvals[j*cv_mem->cv_Ns+is] = ONE;
1844 cv_mem->cv_Xvecs[j*cv_mem->cv_Ns+is] = cv_mem->cv_znQS[j][is];
1845 cv_mem->cv_Zvecs[j*cv_mem->cv_Ns+is] = ck_mem->ck_znQS[j][is];
1846 }
1847 }
1848
1849 (void) N_VScaleVectorArray(cv_mem->cv_Ns, cv_mem->cv_cvals,
1850 cv_mem->cv_Xvecs, cv_mem->cv_Zvecs);
1851
1852 if ( cv_mem->cv_q < qmax ) {
1853 for (is=0; is<cv_mem->cv_Ns; is++)
1854 cv_mem->cv_cvals[is] = ONE;
1855
1856 (void) N_VScaleVectorArray(cv_mem->cv_Ns, cv_mem->cv_cvals,
1857 cv_mem->cv_znQS[qmax], ck_mem->ck_znQS[qmax]);
1858 }
1859 }
1860
1861 for (j=0; j<=L_MAX; j++) ck_mem->ck_tau[j] = cv_mem->cv_tau[j];
1862 for (j=0; j<=NUM_TESTS; j++) ck_mem->ck_tq[j] = cv_mem->cv_tq[j];
1863 for (j=0; j<=cv_mem->cv_q; j++) ck_mem->ck_l[j] = cv_mem->cv_l[j];
1864 ck_mem->ck_nst = cv_mem->cv_nst;
1865 ck_mem->ck_tretlast = cv_mem->cv_tretlast;
1866 ck_mem->ck_q = cv_mem->cv_q;
1867 ck_mem->ck_qprime = cv_mem->cv_qprime;
1868 ck_mem->ck_qwait = cv_mem->cv_qwait;
1869 ck_mem->ck_L = cv_mem->cv_L;
1870 ck_mem->ck_gammap = cv_mem->cv_gammap;
1871 ck_mem->ck_h = cv_mem->cv_h;
1872 ck_mem->ck_hprime = cv_mem->cv_hprime;
1873 ck_mem->ck_hscale = cv_mem->cv_hscale;
1874 ck_mem->ck_eta = cv_mem->cv_eta;
1875 ck_mem->ck_etamax = cv_mem->cv_etamax;
1876 ck_mem->ck_t0 = cv_mem->cv_tn;
1877 ck_mem->ck_saved_tq5 = cv_mem->cv_saved_tq5;
1878
1879 return(ck_mem);
1880 }
1881
1882 /*
1883 * CVAckpntDelete
1884 *
1885 * This routine deletes the first check point in list and returns
1886 * the new list head
1887 */
1888
CVAckpntDelete(CkpntMem * ck_memPtr)1889 static void CVAckpntDelete(CkpntMem *ck_memPtr)
1890 {
1891 CkpntMem tmp;
1892 int j;
1893
1894 if (*ck_memPtr == NULL) return;
1895
1896 /* store head of list */
1897 tmp = *ck_memPtr;
1898
1899 /* move head of list */
1900 *ck_memPtr = (*ck_memPtr)->ck_next;
1901
1902 /* free N_Vectors in tmp */
1903 for (j=0;j<=tmp->ck_q;j++) N_VDestroy(tmp->ck_zn[j]);
1904 if (tmp->ck_zqm != 0) N_VDestroy(tmp->ck_zn[tmp->ck_zqm]);
1905
1906 /* free N_Vectors for quadratures in tmp
1907 * Note that at the check point at t_initial, only znQ_[0]
1908 * was allocated */
1909 if (tmp->ck_quadr) {
1910
1911 if (tmp->ck_next != NULL) {
1912 for (j=0;j<=tmp->ck_q;j++) N_VDestroy(tmp->ck_znQ[j]);
1913 if (tmp->ck_zqm != 0) N_VDestroy(tmp->ck_znQ[tmp->ck_zqm]);
1914 } else {
1915 N_VDestroy(tmp->ck_znQ[0]);
1916 }
1917
1918 }
1919
1920 /* free N_Vectors for sensitivities in tmp
1921 * Note that at the check point at t_initial, only znS_[0]
1922 * was allocated */
1923 if (tmp->ck_sensi) {
1924
1925 if (tmp->ck_next != NULL) {
1926 for (j=0;j<=tmp->ck_q;j++) N_VDestroyVectorArray(tmp->ck_znS[j], tmp->ck_Ns);
1927 if (tmp->ck_zqm != 0) N_VDestroyVectorArray(tmp->ck_znS[tmp->ck_zqm], tmp->ck_Ns);
1928 } else {
1929 N_VDestroyVectorArray(tmp->ck_znS[0], tmp->ck_Ns);
1930 }
1931
1932 }
1933
1934 /* free N_Vectors for quadrature sensitivities in tmp
1935 * Note that at the check point at t_initial, only znQS_[0]
1936 * was allocated */
1937 if (tmp->ck_quadr_sensi) {
1938
1939 if (tmp->ck_next != NULL) {
1940 for (j=0;j<=tmp->ck_q;j++) N_VDestroyVectorArray(tmp->ck_znQS[j], tmp->ck_Ns);
1941 if (tmp->ck_zqm != 0) N_VDestroyVectorArray(tmp->ck_znQS[tmp->ck_zqm], tmp->ck_Ns);
1942 } else {
1943 N_VDestroyVectorArray(tmp->ck_znQS[0], tmp->ck_Ns);
1944 }
1945
1946 }
1947
1948 free(tmp); tmp = NULL;
1949
1950 }
1951
1952 /*
1953 * =================================================================
1954 * PRIVATE FUNCTIONS FOR BACKWARD PROBLEMS
1955 * =================================================================
1956 */
1957
CVAbckpbDelete(CVodeBMem * cvB_memPtr)1958 static void CVAbckpbDelete(CVodeBMem *cvB_memPtr)
1959 {
1960 CVodeBMem tmp;
1961 void *cvode_mem;
1962
1963 if (*cvB_memPtr != NULL) {
1964
1965 /* Save head of the list */
1966 tmp = *cvB_memPtr;
1967
1968 /* Move head of the list */
1969 *cvB_memPtr = (*cvB_memPtr)->cv_next;
1970
1971 /* Free CVODES memory in tmp */
1972 cvode_mem = (void *)(tmp->cv_mem);
1973 CVodeFree(&cvode_mem);
1974
1975 /* Free linear solver memory */
1976 if (tmp->cv_lfree != NULL) tmp->cv_lfree(tmp);
1977
1978 /* Free preconditioner memory */
1979 if (tmp->cv_pfree != NULL) tmp->cv_pfree(tmp);
1980
1981 /* Free workspace Nvector */
1982 N_VDestroy(tmp->cv_y);
1983
1984 free(tmp); tmp = NULL;
1985
1986 }
1987
1988 }
1989
1990 /*
1991 * =================================================================
1992 * PRIVATE FUNCTIONS FOR INTERPOLATION
1993 * =================================================================
1994 */
1995
1996 /*
1997 * CVAdataStore
1998 *
1999 * This routine integrates the forward model starting at the check
2000 * point ck_mem and stores y and yprime at all intermediate steps.
2001 *
2002 * Return values:
2003 * CV_SUCCESS
2004 * CV_REIFWD_FAIL
2005 * CV_FWD_FAIL
2006 */
2007
CVAdataStore(CVodeMem cv_mem,CkpntMem ck_mem)2008 static int CVAdataStore(CVodeMem cv_mem, CkpntMem ck_mem)
2009 {
2010 CVadjMem ca_mem;
2011 DtpntMem *dt_mem;
2012 realtype t;
2013 long int i;
2014 int flag, sign;
2015
2016 ca_mem = cv_mem->cv_adj_mem;
2017 dt_mem = ca_mem->dt_mem;
2018
2019 /* Initialize cv_mem with data from ck_mem */
2020 flag = CVAckpntGet(cv_mem, ck_mem);
2021 if (flag != CV_SUCCESS)
2022 return(CV_REIFWD_FAIL);
2023
2024 /* Set first structure in dt_mem[0] */
2025 dt_mem[0]->t = ck_mem->ck_t0;
2026 ca_mem->ca_IMstore(cv_mem, dt_mem[0]);
2027
2028 /* Decide whether TSTOP must be activated */
2029 if (ca_mem->ca_tstopCVodeFcall) {
2030 CVodeSetStopTime(cv_mem, ca_mem->ca_tstopCVodeF);
2031 }
2032
2033 sign = (ca_mem->ca_tfinal - ca_mem->ca_tinitial > ZERO) ? 1 : -1;
2034
2035
2036 /* Run CVode to set following structures in dt_mem[i] */
2037 i = 1;
2038 do {
2039
2040 flag = CVode(cv_mem, ck_mem->ck_t1, ca_mem->ca_ytmp, &t, CV_ONE_STEP);
2041 if (flag < 0) return(CV_FWD_FAIL);
2042
2043 dt_mem[i]->t = t;
2044 ca_mem->ca_IMstore(cv_mem, dt_mem[i]);
2045 i++;
2046
2047 } while ( sign*(ck_mem->ck_t1 - t) > ZERO );
2048
2049
2050 ca_mem->ca_IMnewData = SUNTRUE; /* New data is now available */
2051 ca_mem->ca_ckpntData = ck_mem; /* starting at this check point */
2052 ca_mem->ca_np = i; /* and we have this many points */
2053
2054 return(CV_SUCCESS);
2055 }
2056
2057 /*
2058 * CVAckpntGet
2059 *
2060 * This routine prepares CVODES for a hot restart from
2061 * the check point ck_mem
2062 */
2063
CVAckpntGet(CVodeMem cv_mem,CkpntMem ck_mem)2064 static int CVAckpntGet(CVodeMem cv_mem, CkpntMem ck_mem)
2065 {
2066 int flag, j, is, qmax, retval;
2067
2068 if (ck_mem->ck_next == NULL) {
2069
2070 /* In this case, we just call the reinitialization routine,
2071 * but make sure we use the same initial stepsize as on
2072 * the first run. */
2073
2074 CVodeSetInitStep(cv_mem, cv_mem->cv_h0u);
2075
2076 flag = CVodeReInit(cv_mem, ck_mem->ck_t0, ck_mem->ck_zn[0]);
2077 if (flag != CV_SUCCESS) return(flag);
2078
2079 if (ck_mem->ck_quadr) {
2080 flag = CVodeQuadReInit(cv_mem, ck_mem->ck_znQ[0]);
2081 if (flag != CV_SUCCESS) return(flag);
2082 }
2083
2084 if (ck_mem->ck_sensi) {
2085 flag = CVodeSensReInit(cv_mem, cv_mem->cv_ism, ck_mem->ck_znS[0]);
2086 if (flag != CV_SUCCESS) return(flag);
2087 }
2088
2089 if (ck_mem->ck_quadr_sensi) {
2090 flag = CVodeQuadSensReInit(cv_mem, ck_mem->ck_znQS[0]);
2091 if (flag != CV_SUCCESS) return(flag);
2092 }
2093
2094 } else {
2095
2096 qmax = cv_mem->cv_qmax;
2097
2098 /* Copy parameters from check point data structure */
2099
2100 cv_mem->cv_nst = ck_mem->ck_nst;
2101 cv_mem->cv_tretlast = ck_mem->ck_tretlast;
2102 cv_mem->cv_q = ck_mem->ck_q;
2103 cv_mem->cv_qprime = ck_mem->ck_qprime;
2104 cv_mem->cv_qwait = ck_mem->ck_qwait;
2105 cv_mem->cv_L = ck_mem->ck_L;
2106 cv_mem->cv_gammap = ck_mem->ck_gammap;
2107 cv_mem->cv_h = ck_mem->ck_h;
2108 cv_mem->cv_hprime = ck_mem->ck_hprime;
2109 cv_mem->cv_hscale = ck_mem->ck_hscale;
2110 cv_mem->cv_eta = ck_mem->ck_eta;
2111 cv_mem->cv_etamax = ck_mem->ck_etamax;
2112 cv_mem->cv_tn = ck_mem->ck_t0;
2113 cv_mem->cv_saved_tq5 = ck_mem->ck_saved_tq5;
2114
2115 /* Copy the arrays from check point data structure */
2116
2117 for (j=0; j<=cv_mem->cv_q; j++)
2118 cv_mem->cv_cvals[j] = ONE;
2119
2120 retval = N_VScaleVectorArray(cv_mem->cv_q+1, cv_mem->cv_cvals,
2121 ck_mem->ck_zn, cv_mem->cv_zn);
2122 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2123
2124 if ( cv_mem->cv_q < qmax )
2125 N_VScale(ONE, ck_mem->ck_zn[qmax], cv_mem->cv_zn[qmax]);
2126
2127 if (ck_mem->ck_quadr) {
2128 for (j=0; j<=cv_mem->cv_q; j++)
2129 cv_mem->cv_cvals[j] = ONE;
2130
2131 retval = N_VScaleVectorArray(cv_mem->cv_q+1, cv_mem->cv_cvals,
2132 ck_mem->ck_znQ, cv_mem->cv_znQ);
2133 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2134
2135 if ( cv_mem->cv_q < qmax )
2136 N_VScale(ONE, ck_mem->ck_znQ[qmax], cv_mem->cv_znQ[qmax]);
2137 }
2138
2139 if (ck_mem->ck_sensi) {
2140 for (j=0; j<=cv_mem->cv_q; j++) {
2141 for (is=0; is<cv_mem->cv_Ns; is++) {
2142 cv_mem->cv_cvals[j*cv_mem->cv_Ns+is] = ONE;
2143 cv_mem->cv_Xvecs[j*cv_mem->cv_Ns+is] = ck_mem->ck_znS[j][is];
2144 cv_mem->cv_Zvecs[j*cv_mem->cv_Ns+is] = cv_mem->cv_znS[j][is];
2145 }
2146 }
2147
2148 retval = N_VScaleVectorArray(cv_mem->cv_Ns*(cv_mem->cv_q+1),
2149 cv_mem->cv_cvals,
2150 cv_mem->cv_Xvecs, cv_mem->cv_Zvecs);
2151 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2152
2153 if ( cv_mem->cv_q < qmax ) {
2154 for (is=0; is<cv_mem->cv_Ns; is++)
2155 cv_mem->cv_cvals[is] = ONE;
2156
2157 retval = N_VScaleVectorArray(cv_mem->cv_Ns, cv_mem->cv_cvals,
2158 ck_mem->ck_znS[qmax], cv_mem->cv_znS[qmax]);
2159 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2160 }
2161 }
2162
2163 if (ck_mem->ck_quadr_sensi) {
2164 for (j=0; j<=cv_mem->cv_q; j++) {
2165 for (is=0; is<cv_mem->cv_Ns; is++) {
2166 cv_mem->cv_cvals[j*cv_mem->cv_Ns+is] = ONE;
2167 cv_mem->cv_Xvecs[j*cv_mem->cv_Ns+is] = ck_mem->ck_znQS[j][is];
2168 cv_mem->cv_Zvecs[j*cv_mem->cv_Ns+is] = cv_mem->cv_znQS[j][is];
2169 }
2170 }
2171
2172 retval = N_VScaleVectorArray(cv_mem->cv_Ns*(cv_mem->cv_q+1),
2173 cv_mem->cv_cvals,
2174 cv_mem->cv_Xvecs, cv_mem->cv_Zvecs);
2175 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2176
2177 if ( cv_mem->cv_q < qmax ) {
2178 for (is=0; is<cv_mem->cv_Ns; is++)
2179 cv_mem->cv_cvals[is] = ONE;
2180
2181 retval = N_VScaleVectorArray(cv_mem->cv_Ns, cv_mem->cv_cvals,
2182 ck_mem->ck_znQS[qmax], cv_mem->cv_znQS[qmax]);
2183 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2184 }
2185 }
2186
2187 for (j=0; j<=L_MAX; j++) cv_mem->cv_tau[j] = ck_mem->ck_tau[j];
2188 for (j=0; j<=NUM_TESTS; j++) cv_mem->cv_tq[j] = ck_mem->ck_tq[j];
2189 for (j=0; j<=cv_mem->cv_q; j++) cv_mem->cv_l[j] = ck_mem->ck_l[j];
2190
2191 /* Force a call to setup */
2192
2193 cv_mem->cv_forceSetup = SUNTRUE;
2194
2195 }
2196
2197 return(CV_SUCCESS);
2198 }
2199
2200 /*
2201 * -----------------------------------------------------------------
2202 * Functions for interpolation
2203 * -----------------------------------------------------------------
2204 */
2205
2206 /*
2207 * CVAfindIndex
2208 *
2209 * Finds the index in the array of data point strctures such that
2210 * dt_mem[indx-1].t <= t < dt_mem[indx].t
2211 * If indx is changed from the previous invocation, then newpoint = SUNTRUE
2212 *
2213 * If t is beyond the leftmost limit, but close enough, indx=0.
2214 *
2215 * Returns CV_SUCCESS if successful and CV_GETY_BADT if unable to
2216 * find indx (t is too far beyond limits).
2217 */
2218
CVAfindIndex(CVodeMem cv_mem,realtype t,long int * indx,booleantype * newpoint)2219 static int CVAfindIndex(CVodeMem cv_mem, realtype t,
2220 long int *indx, booleantype *newpoint)
2221 {
2222 CVadjMem ca_mem;
2223 DtpntMem *dt_mem;
2224 int sign;
2225 booleantype to_left, to_right;
2226
2227 ca_mem = cv_mem->cv_adj_mem;
2228 dt_mem = ca_mem->dt_mem;
2229
2230 *newpoint = SUNFALSE;
2231
2232 /* Find the direction of integration */
2233 sign = (ca_mem->ca_tfinal - ca_mem->ca_tinitial > ZERO) ? 1 : -1;
2234
2235 /* If this is the first time we use new data */
2236 if (ca_mem->ca_IMnewData) {
2237 ca_mem->ca_ilast = ca_mem->ca_np-1;
2238 *newpoint = SUNTRUE;
2239 ca_mem->ca_IMnewData = SUNFALSE;
2240 }
2241
2242 /* Search for indx starting from ilast */
2243 to_left = ( sign*(t - dt_mem[ca_mem->ca_ilast-1]->t) < ZERO);
2244 to_right = ( sign*(t - dt_mem[ca_mem->ca_ilast]->t) > ZERO);
2245
2246 if ( to_left ) {
2247 /* look for a new indx to the left */
2248
2249 *newpoint = SUNTRUE;
2250
2251 *indx = ca_mem->ca_ilast;
2252 for(;;) {
2253 if ( *indx == 0 ) break;
2254 if ( sign*(t - dt_mem[*indx-1]->t) <= ZERO ) (*indx)--;
2255 else break;
2256 }
2257
2258 if ( *indx == 0 )
2259 ca_mem->ca_ilast = 1;
2260 else
2261 ca_mem->ca_ilast = *indx;
2262
2263 if ( *indx == 0 ) {
2264 /* t is beyond leftmost limit. Is it too far? */
2265 if ( SUNRabs(t - dt_mem[0]->t) > FUZZ_FACTOR * cv_mem->cv_uround ) {
2266 return(CV_GETY_BADT);
2267 }
2268 }
2269
2270 } else if ( to_right ) {
2271 /* look for a new indx to the right */
2272
2273 *newpoint = SUNTRUE;
2274
2275 *indx = ca_mem->ca_ilast;
2276 for(;;) {
2277 if ( sign*(t - dt_mem[*indx]->t) > ZERO) (*indx)++;
2278 else break;
2279 }
2280
2281 ca_mem->ca_ilast = *indx;
2282
2283
2284 } else {
2285 /* ilast is still OK */
2286
2287 *indx = ca_mem->ca_ilast;
2288
2289 }
2290
2291 return(CV_SUCCESS);
2292
2293
2294 }
2295
2296 /*
2297 * CVodeGetAdjY
2298 *
2299 * This routine returns the interpolated forward solution at time t.
2300 * The user must allocate space for y.
2301 */
2302
CVodeGetAdjY(void * cvode_mem,realtype t,N_Vector y)2303 int CVodeGetAdjY(void *cvode_mem, realtype t, N_Vector y)
2304 {
2305 CVodeMem cv_mem;
2306 CVadjMem ca_mem;
2307 int flag;
2308
2309 if (cvode_mem == NULL) {
2310 cvProcessError(NULL, CV_MEM_NULL, "CVODEA", "CVodeGetAdjY", MSGCV_NO_MEM);
2311 return(CV_MEM_NULL);
2312 }
2313 cv_mem = (CVodeMem) cvode_mem;
2314
2315 ca_mem = cv_mem->cv_adj_mem;
2316
2317 flag = ca_mem->ca_IMget(cv_mem, t, y, NULL);
2318
2319 return(flag);
2320 }
2321
2322 /*
2323 * -----------------------------------------------------------------
2324 * Functions specific to cubic Hermite interpolation
2325 * -----------------------------------------------------------------
2326 */
2327
2328 /*
2329 * CVAhermiteMalloc
2330 *
2331 * This routine allocates memory for storing information at all
2332 * intermediate points between two consecutive check points.
2333 * This data is then used to interpolate the forward solution
2334 * at any other time.
2335 */
2336
CVAhermiteMalloc(CVodeMem cv_mem)2337 static booleantype CVAhermiteMalloc(CVodeMem cv_mem)
2338 {
2339 CVadjMem ca_mem;
2340 DtpntMem *dt_mem;
2341 HermiteDataMem content;
2342 long int i, ii=0;
2343 booleantype allocOK;
2344
2345 allocOK = SUNTRUE;
2346
2347 ca_mem = cv_mem->cv_adj_mem;
2348
2349 /* Allocate space for the vectors ytmp and yStmp */
2350
2351 ca_mem->ca_ytmp = N_VClone(cv_mem->cv_tempv);
2352 if (ca_mem->ca_ytmp == NULL) {
2353 return(SUNFALSE);
2354 }
2355
2356 if (ca_mem->ca_IMstoreSensi) {
2357 ca_mem->ca_yStmp = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempv);
2358 if (ca_mem->ca_yStmp == NULL) {
2359 N_VDestroy(ca_mem->ca_ytmp);
2360 return(SUNFALSE);
2361 }
2362 }
2363
2364 /* Allocate space for the content field of the dt structures */
2365
2366 dt_mem = ca_mem->dt_mem;
2367
2368 for (i=0; i<=ca_mem->ca_nsteps; i++) {
2369
2370 content = NULL;
2371 content = (HermiteDataMem) malloc(sizeof(struct HermiteDataMemRec));
2372 if (content == NULL) {
2373 ii = i;
2374 allocOK = SUNFALSE;
2375 break;
2376 }
2377
2378 content->y = N_VClone(cv_mem->cv_tempv);
2379 if (content->y == NULL) {
2380 free(content); content = NULL;
2381 ii = i;
2382 allocOK = SUNFALSE;
2383 break;
2384 }
2385
2386 content->yd = N_VClone(cv_mem->cv_tempv);
2387 if (content->yd == NULL) {
2388 N_VDestroy(content->y);
2389 free(content); content = NULL;
2390 ii = i;
2391 allocOK = SUNFALSE;
2392 break;
2393 }
2394
2395 if (ca_mem->ca_IMstoreSensi) {
2396
2397 content->yS = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempv);
2398 if (content->yS == NULL) {
2399 N_VDestroy(content->y);
2400 N_VDestroy(content->yd);
2401 free(content); content = NULL;
2402 ii = i;
2403 allocOK = SUNFALSE;
2404 break;
2405 }
2406
2407 content->ySd = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempv);
2408 if (content->ySd == NULL) {
2409 N_VDestroy(content->y);
2410 N_VDestroy(content->yd);
2411 N_VDestroyVectorArray(content->yS, cv_mem->cv_Ns);
2412 free(content); content = NULL;
2413 ii = i;
2414 allocOK = SUNFALSE;
2415 break;
2416 }
2417
2418 }
2419
2420 dt_mem[i]->content = content;
2421
2422 }
2423
2424 /* If an error occurred, deallocate and return */
2425
2426 if (!allocOK) {
2427
2428 N_VDestroy(ca_mem->ca_ytmp);
2429
2430 if (ca_mem->ca_IMstoreSensi) {
2431 N_VDestroyVectorArray(ca_mem->ca_yStmp, cv_mem->cv_Ns);
2432 }
2433
2434 for (i=0; i<ii; i++) {
2435 content = (HermiteDataMem) (dt_mem[i]->content);
2436 N_VDestroy(content->y);
2437 N_VDestroy(content->yd);
2438 if (ca_mem->ca_IMstoreSensi) {
2439 N_VDestroyVectorArray(content->yS, cv_mem->cv_Ns);
2440 N_VDestroyVectorArray(content->ySd, cv_mem->cv_Ns);
2441 }
2442 free(dt_mem[i]->content); dt_mem[i]->content = NULL;
2443 }
2444
2445 }
2446
2447 return(allocOK);
2448 }
2449
2450 /*
2451 * CVAhermiteFree
2452 *
2453 * This routine frees the memory allocated for data storage.
2454 */
2455
CVAhermiteFree(CVodeMem cv_mem)2456 static void CVAhermiteFree(CVodeMem cv_mem)
2457 {
2458 CVadjMem ca_mem;
2459 DtpntMem *dt_mem;
2460 HermiteDataMem content;
2461 long int i;
2462
2463 ca_mem = cv_mem->cv_adj_mem;
2464
2465 N_VDestroy(ca_mem->ca_ytmp);
2466
2467 if (ca_mem->ca_IMstoreSensi) {
2468 N_VDestroyVectorArray(ca_mem->ca_yStmp, cv_mem->cv_Ns);
2469 }
2470
2471 dt_mem = ca_mem->dt_mem;
2472
2473 for (i=0; i<=ca_mem->ca_nsteps; i++) {
2474 content = (HermiteDataMem) (dt_mem[i]->content);
2475 N_VDestroy(content->y);
2476 N_VDestroy(content->yd);
2477 if (ca_mem->ca_IMstoreSensi) {
2478 N_VDestroyVectorArray(content->yS, cv_mem->cv_Ns);
2479 N_VDestroyVectorArray(content->ySd, cv_mem->cv_Ns);
2480 }
2481 free(dt_mem[i]->content); dt_mem[i]->content = NULL;
2482 }
2483 }
2484
2485 /*
2486 * CVAhermiteStorePnt ( -> IMstore )
2487 *
2488 * This routine stores a new point (y,yd) in the structure d for use
2489 * in the cubic Hermite interpolation.
2490 * Note that the time is already stored.
2491 */
2492
CVAhermiteStorePnt(CVodeMem cv_mem,DtpntMem d)2493 static int CVAhermiteStorePnt(CVodeMem cv_mem, DtpntMem d)
2494 {
2495 CVadjMem ca_mem;
2496 HermiteDataMem content;
2497 int is, retval;
2498
2499 ca_mem = cv_mem->cv_adj_mem;
2500
2501 content = (HermiteDataMem) d->content;
2502
2503 /* Load solution */
2504
2505 N_VScale(ONE, cv_mem->cv_zn[0], content->y);
2506
2507 if (ca_mem->ca_IMstoreSensi) {
2508 for (is=0; is<cv_mem->cv_Ns; is++)
2509 cv_mem->cv_cvals[is] = ONE;
2510
2511 retval = N_VScaleVectorArray(cv_mem->cv_Ns, cv_mem->cv_cvals,
2512 cv_mem->cv_znS[0], content->yS);
2513 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2514 }
2515
2516 /* Load derivative */
2517
2518 if (cv_mem->cv_nst == 0) {
2519
2520 /* retval = */ cv_mem->cv_f(cv_mem->cv_tn, content->y, content->yd, cv_mem->cv_user_data);
2521
2522 if (ca_mem->ca_IMstoreSensi) {
2523 /* retval = */ cvSensRhsWrapper(cv_mem, cv_mem->cv_tn, content->y, content->yd,
2524 content->yS, content->ySd,
2525 cv_mem->cv_tempv, cv_mem->cv_ftemp);
2526 }
2527
2528 } else {
2529
2530 N_VScale(ONE/cv_mem->cv_h, cv_mem->cv_zn[1], content->yd);
2531
2532 if (ca_mem->ca_IMstoreSensi) {
2533 for (is=0; is<cv_mem->cv_Ns; is++)
2534 cv_mem->cv_cvals[is] = ONE/cv_mem->cv_h;
2535
2536 retval = N_VScaleVectorArray(cv_mem->cv_Ns, cv_mem->cv_cvals,
2537 cv_mem->cv_znS[1], content->ySd);
2538 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2539 }
2540
2541 }
2542
2543 return(0);
2544 }
2545
2546 /*
2547 * CVAhermiteGetY ( -> IMget )
2548 *
2549 * This routine uses cubic piece-wise Hermite interpolation for
2550 * the forward solution vector.
2551 * It is typically called by the wrapper routines before calling
2552 * user provided routines (fB, djacB, bjacB, jtimesB, psolB) but
2553 * can be directly called by the user through CVodeGetAdjY
2554 */
2555
CVAhermiteGetY(CVodeMem cv_mem,realtype t,N_Vector y,N_Vector * yS)2556 static int CVAhermiteGetY(CVodeMem cv_mem, realtype t,
2557 N_Vector y, N_Vector *yS)
2558 {
2559 CVadjMem ca_mem;
2560 DtpntMem *dt_mem;
2561 HermiteDataMem content0, content1;
2562
2563 realtype t0, t1, delta;
2564 realtype factor1, factor2, factor3;
2565
2566 N_Vector y0, yd0, y1, yd1;
2567 N_Vector *yS0=NULL, *ySd0=NULL, *yS1, *ySd1;
2568
2569 int flag, is, NS;
2570 long int indx;
2571 booleantype newpoint;
2572
2573 /* local variables for fused vector oerations */
2574 int retval;
2575 realtype cvals[4];
2576 N_Vector Xvecs[4];
2577 N_Vector* XXvecs[4];
2578
2579 ca_mem = cv_mem->cv_adj_mem;
2580 dt_mem = ca_mem->dt_mem;
2581
2582 /* Local value of Ns */
2583
2584 NS = (ca_mem->ca_IMinterpSensi && (yS != NULL)) ? cv_mem->cv_Ns : 0;
2585
2586 /* Get the index in dt_mem */
2587
2588 flag = CVAfindIndex(cv_mem, t, &indx, &newpoint);
2589 if (flag != CV_SUCCESS) return(flag);
2590
2591 /* If we are beyond the left limit but close enough,
2592 then return y at the left limit. */
2593
2594 if (indx == 0) {
2595 content0 = (HermiteDataMem) (dt_mem[0]->content);
2596 N_VScale(ONE, content0->y, y);
2597
2598 if (NS > 0) {
2599 for (is=0; is<NS; is++)
2600 cv_mem->cv_cvals[is] = ONE;
2601
2602 retval = N_VScaleVectorArray(NS, cv_mem->cv_cvals,
2603 content0->yS, yS);
2604 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2605 }
2606
2607 return(CV_SUCCESS);
2608 }
2609
2610 /* Extract stuff from the appropriate data points */
2611
2612 t0 = dt_mem[indx-1]->t;
2613 t1 = dt_mem[indx]->t;
2614 delta = t1 - t0;
2615
2616 content0 = (HermiteDataMem) (dt_mem[indx-1]->content);
2617 y0 = content0->y;
2618 yd0 = content0->yd;
2619 if (ca_mem->ca_IMinterpSensi) {
2620 yS0 = content0->yS;
2621 ySd0 = content0->ySd;
2622 }
2623
2624 if (newpoint) {
2625
2626 /* Recompute Y0 and Y1 */
2627
2628 content1 = (HermiteDataMem) (dt_mem[indx]->content);
2629
2630 y1 = content1->y;
2631 yd1 = content1->yd;
2632
2633 /* Y1 = delta (yd1 + yd0) - 2 (y1 - y0) */
2634 cvals[0] = -TWO; Xvecs[0] = y1;
2635 cvals[1] = TWO; Xvecs[1] = y0;
2636 cvals[2] = delta; Xvecs[2] = yd1;
2637 cvals[3] = delta; Xvecs[3] = yd0;
2638
2639 retval = N_VLinearCombination(4, cvals, Xvecs, ca_mem->ca_Y[1]);
2640 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2641
2642 /* Y0 = y1 - y0 - delta * yd0 */
2643 cvals[0] = ONE; Xvecs[0] = y1;
2644 cvals[1] = -ONE; Xvecs[1] = y0;
2645 cvals[2] = -delta; Xvecs[2] = yd0;
2646
2647 retval = N_VLinearCombination(3, cvals, Xvecs, ca_mem->ca_Y[0]);
2648 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2649
2650 /* Recompute YS0 and YS1, if needed */
2651
2652 if (NS > 0) {
2653
2654 yS1 = content1->yS;
2655 ySd1 = content1->ySd;
2656
2657 /* YS1 = delta (ySd1 + ySd0) - 2 (yS1 - yS0) */
2658 cvals[0] = -TWO; XXvecs[0] = yS1;
2659 cvals[1] = TWO; XXvecs[1] = yS0;
2660 cvals[2] = delta; XXvecs[2] = ySd1;
2661 cvals[3] = delta; XXvecs[3] = ySd0;
2662
2663 retval = N_VLinearCombinationVectorArray(NS, 4, cvals, XXvecs, ca_mem->ca_YS[1]);
2664 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2665
2666 /* YS0 = yS1 - yS0 - delta * ySd0 */
2667 cvals[0] = ONE; XXvecs[0] = yS1;
2668 cvals[1] = -ONE; XXvecs[1] = yS0;
2669 cvals[2] = -delta; XXvecs[2] = ySd0;
2670
2671 retval = N_VLinearCombinationVectorArray(NS, 3, cvals, XXvecs, ca_mem->ca_YS[0]);
2672 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2673
2674 }
2675
2676 }
2677
2678 /* Perform the actual interpolation. */
2679
2680 factor1 = t - t0;
2681
2682 factor2 = factor1/delta;
2683 factor2 = factor2*factor2;
2684
2685 factor3 = factor2*(t-t1)/delta;
2686
2687 cvals[0] = ONE;
2688 cvals[1] = factor1;
2689 cvals[2] = factor2;
2690 cvals[3] = factor3;
2691
2692 /* y = y0 + factor1 yd0 + factor2 * Y[0] + factor3 Y[1] */
2693 Xvecs[0] = y0;
2694 Xvecs[1] = yd0;
2695 Xvecs[2] = ca_mem->ca_Y[0];
2696 Xvecs[3] = ca_mem->ca_Y[1];
2697
2698 retval = N_VLinearCombination(4, cvals, Xvecs, y);
2699 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2700
2701 /* yS = yS0 + factor1 ySd0 + factor2 * YS[0] + factor3 YS[1], if needed */
2702 if (NS > 0) {
2703
2704 XXvecs[0] = yS0;
2705 XXvecs[1] = ySd0;
2706 XXvecs[2] = ca_mem->ca_YS[0];
2707 XXvecs[3] = ca_mem->ca_YS[1];
2708
2709 retval = N_VLinearCombinationVectorArray(NS, 4, cvals, XXvecs, yS);
2710 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2711
2712 }
2713
2714 return(CV_SUCCESS);
2715 }
2716
2717 /*
2718 * -----------------------------------------------------------------
2719 * Functions specific to Polynomial interpolation
2720 * -----------------------------------------------------------------
2721 */
2722
2723 /*
2724 * CVApolynomialMalloc
2725 *
2726 * This routine allocates memory for storing information at all
2727 * intermediate points between two consecutive check points.
2728 * This data is then used to interpolate the forward solution
2729 * at any other time.
2730 */
2731
CVApolynomialMalloc(CVodeMem cv_mem)2732 static booleantype CVApolynomialMalloc(CVodeMem cv_mem)
2733 {
2734 CVadjMem ca_mem;
2735 DtpntMem *dt_mem;
2736 PolynomialDataMem content;
2737 long int i, ii=0;
2738 booleantype allocOK;
2739
2740 allocOK = SUNTRUE;
2741
2742 ca_mem = cv_mem->cv_adj_mem;
2743
2744 /* Allocate space for the vectors ytmp and yStmp */
2745
2746 ca_mem->ca_ytmp = N_VClone(cv_mem->cv_tempv);
2747 if (ca_mem->ca_ytmp == NULL) {
2748 return(SUNFALSE);
2749 }
2750
2751 if (ca_mem->ca_IMstoreSensi) {
2752 ca_mem->ca_yStmp = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempv);
2753 if (ca_mem->ca_yStmp == NULL) {
2754 N_VDestroy(ca_mem->ca_ytmp);
2755 return(SUNFALSE);
2756 }
2757 }
2758
2759 /* Allocate space for the content field of the dt structures */
2760
2761 dt_mem = ca_mem->dt_mem;
2762
2763 for (i=0; i<=ca_mem->ca_nsteps; i++) {
2764
2765 content = NULL;
2766 content = (PolynomialDataMem) malloc(sizeof(struct PolynomialDataMemRec));
2767 if (content == NULL) {
2768 ii = i;
2769 allocOK = SUNFALSE;
2770 break;
2771 }
2772
2773 content->y = N_VClone(cv_mem->cv_tempv);
2774 if (content->y == NULL) {
2775 free(content); content = NULL;
2776 ii = i;
2777 allocOK = SUNFALSE;
2778 break;
2779 }
2780
2781 if (ca_mem->ca_IMstoreSensi) {
2782
2783 content->yS = N_VCloneVectorArray(cv_mem->cv_Ns, cv_mem->cv_tempv);
2784 if (content->yS == NULL) {
2785 N_VDestroy(content->y);
2786 free(content); content = NULL;
2787 ii = i;
2788 allocOK = SUNFALSE;
2789 break;
2790 }
2791
2792 }
2793
2794 dt_mem[i]->content = content;
2795
2796 }
2797
2798 /* If an error occurred, deallocate and return */
2799
2800 if (!allocOK) {
2801
2802 N_VDestroy(ca_mem->ca_ytmp);
2803
2804 if (ca_mem->ca_IMstoreSensi) {
2805 N_VDestroyVectorArray(ca_mem->ca_yStmp, cv_mem->cv_Ns);
2806 }
2807
2808 for (i=0; i<ii; i++) {
2809 content = (PolynomialDataMem) (dt_mem[i]->content);
2810 N_VDestroy(content->y);
2811 if (ca_mem->ca_IMstoreSensi) {
2812 N_VDestroyVectorArray(content->yS, cv_mem->cv_Ns);
2813 }
2814 free(dt_mem[i]->content); dt_mem[i]->content = NULL;
2815 }
2816
2817 }
2818
2819 return(allocOK);
2820
2821 }
2822
2823 /*
2824 * CVApolynomialFree
2825 *
2826 * This routine frees the memeory allocated for data storage.
2827 */
2828
CVApolynomialFree(CVodeMem cv_mem)2829 static void CVApolynomialFree(CVodeMem cv_mem)
2830 {
2831 CVadjMem ca_mem;
2832 DtpntMem *dt_mem;
2833 PolynomialDataMem content;
2834 long int i;
2835
2836 ca_mem = cv_mem->cv_adj_mem;
2837
2838 N_VDestroy(ca_mem->ca_ytmp);
2839
2840 if (ca_mem->ca_IMstoreSensi) {
2841 N_VDestroyVectorArray(ca_mem->ca_yStmp, cv_mem->cv_Ns);
2842 }
2843
2844 dt_mem = ca_mem->dt_mem;
2845
2846 for (i=0; i<=ca_mem->ca_nsteps; i++) {
2847 content = (PolynomialDataMem) (dt_mem[i]->content);
2848 N_VDestroy(content->y);
2849 if (ca_mem->ca_IMstoreSensi) {
2850 N_VDestroyVectorArray(content->yS, cv_mem->cv_Ns);
2851 }
2852 free(dt_mem[i]->content); dt_mem[i]->content = NULL;
2853 }
2854 }
2855
2856 /*
2857 * CVApolynomialStorePnt ( -> IMstore )
2858 *
2859 * This routine stores a new point y in the structure d for use
2860 * in the Polynomial interpolation.
2861 * Note that the time is already stored.
2862 */
2863
CVApolynomialStorePnt(CVodeMem cv_mem,DtpntMem d)2864 static int CVApolynomialStorePnt(CVodeMem cv_mem, DtpntMem d)
2865 {
2866 CVadjMem ca_mem;
2867 PolynomialDataMem content;
2868 int is, retval;
2869
2870 ca_mem = cv_mem->cv_adj_mem;
2871
2872 content = (PolynomialDataMem) d->content;
2873
2874 N_VScale(ONE, cv_mem->cv_zn[0], content->y);
2875
2876 if (ca_mem->ca_IMstoreSensi) {
2877 for (is=0; is<cv_mem->cv_Ns; is++)
2878 cv_mem->cv_cvals[is] = ONE;
2879 retval = N_VScaleVectorArray(cv_mem->cv_Ns, cv_mem->cv_cvals,
2880 cv_mem->cv_znS[0], content->yS);
2881 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2882 }
2883
2884 content->order = cv_mem->cv_qu;
2885
2886 return(0);
2887 }
2888
2889 /*
2890 * CVApolynomialGetY ( -> IMget )
2891 *
2892 * This routine uses polynomial interpolation for the forward solution vector.
2893 * It is typically called by the wrapper routines before calling
2894 * user provided routines (fB, djacB, bjacB, jtimesB, psolB)) but
2895 * can be directly called by the user through CVodeGetAdjY.
2896 */
2897
CVApolynomialGetY(CVodeMem cv_mem,realtype t,N_Vector y,N_Vector * yS)2898 static int CVApolynomialGetY(CVodeMem cv_mem, realtype t,
2899 N_Vector y, N_Vector *yS)
2900 {
2901 CVadjMem ca_mem;
2902 DtpntMem *dt_mem;
2903 PolynomialDataMem content;
2904
2905 int flag, dir, order, i, j, is, NS, retval;
2906 long int indx, base;
2907 booleantype newpoint;
2908 realtype dt, factor;
2909
2910 ca_mem = cv_mem->cv_adj_mem;
2911 dt_mem = ca_mem->dt_mem;
2912
2913 /* Local value of Ns */
2914
2915 NS = (ca_mem->ca_IMinterpSensi && (yS != NULL)) ? cv_mem->cv_Ns : 0;
2916
2917 /* Get the index in dt_mem */
2918
2919 flag = CVAfindIndex(cv_mem, t, &indx, &newpoint);
2920 if (flag != CV_SUCCESS) return(flag);
2921
2922 /* If we are beyond the left limit but close enough,
2923 then return y at the left limit. */
2924
2925 if (indx == 0) {
2926 content = (PolynomialDataMem) (dt_mem[0]->content);
2927 N_VScale(ONE, content->y, y);
2928
2929 if (NS > 0) {
2930 for (is=0; is<NS; is++)
2931 cv_mem->cv_cvals[is] = ONE;
2932 retval = N_VScaleVectorArray(NS, cv_mem->cv_cvals, content->yS, yS);
2933 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2934 }
2935
2936 return(CV_SUCCESS);
2937 }
2938
2939 /* Scaling factor */
2940
2941 dt = SUNRabs(dt_mem[indx]->t - dt_mem[indx-1]->t);
2942
2943 /* Find the direction of the forward integration */
2944
2945 dir = (ca_mem->ca_tfinal - ca_mem->ca_tinitial > ZERO) ? 1 : -1;
2946
2947 /* Establish the base point depending on the integration direction.
2948 Modify the base if there are not enough points for the current order */
2949
2950 if (dir == 1) {
2951 base = indx;
2952 content = (PolynomialDataMem) (dt_mem[base]->content);
2953 order = content->order;
2954 if(indx < order) base += order-indx;
2955 } else {
2956 base = indx-1;
2957 content = (PolynomialDataMem) (dt_mem[base]->content);
2958 order = content->order;
2959 if (ca_mem->ca_np-indx > order) base -= indx+order-ca_mem->ca_np;
2960 }
2961
2962 /* Recompute Y (divided differences for Newton polynomial) if needed */
2963
2964 if (newpoint) {
2965
2966 /* Store 0-th order DD */
2967 if (dir == 1) {
2968 for(j=0;j<=order;j++) {
2969 ca_mem->ca_T[j] = dt_mem[base-j]->t;
2970 content = (PolynomialDataMem) (dt_mem[base-j]->content);
2971 N_VScale(ONE, content->y, ca_mem->ca_Y[j]);
2972
2973 if (NS > 0) {
2974 for (is=0; is<NS; is++)
2975 cv_mem->cv_cvals[is] = ONE;
2976 retval = N_VScaleVectorArray(NS, cv_mem->cv_cvals,
2977 content->yS, ca_mem->ca_YS[j]);
2978 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2979 }
2980 }
2981 } else {
2982 for(j=0;j<=order;j++) {
2983 ca_mem->ca_T[j] = dt_mem[base-1+j]->t;
2984 content = (PolynomialDataMem) (dt_mem[base-1+j]->content);
2985 N_VScale(ONE, content->y, ca_mem->ca_Y[j]);
2986 if (NS > 0) {
2987 for (is=0; is<NS; is++)
2988 cv_mem->cv_cvals[is] = ONE;
2989 retval = N_VScaleVectorArray(NS, cv_mem->cv_cvals,
2990 content->yS, ca_mem->ca_YS[j]);
2991 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
2992 }
2993 }
2994 }
2995
2996 /* Compute higher-order DD */
2997 for(i=1;i<=order;i++) {
2998 for(j=order;j>=i;j--) {
2999 factor = dt/(ca_mem->ca_T[j]-ca_mem->ca_T[j-i]);
3000 N_VLinearSum(factor, ca_mem->ca_Y[j], -factor, ca_mem->ca_Y[j-1], ca_mem->ca_Y[j]);
3001
3002 if (NS > 0) {
3003 retval = N_VLinearSumVectorArray(NS,
3004 factor, ca_mem->ca_YS[j],
3005 -factor, ca_mem->ca_YS[j-1],
3006 ca_mem->ca_YS[j]);
3007 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
3008 }
3009 }
3010 }
3011 }
3012
3013 /* Perform the actual interpolation using nested multiplications */
3014
3015 cv_mem->cv_cvals[0] = ONE;
3016 for (i=0; i<order; i++)
3017 cv_mem->cv_cvals[i+1] = cv_mem->cv_cvals[i] * (t-ca_mem->ca_T[i]) / dt;
3018
3019 retval = N_VLinearCombination(order+1, cv_mem->cv_cvals, ca_mem->ca_Y, y);
3020 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
3021
3022 if (NS > 0) {
3023 retval = N_VLinearCombinationVectorArray(NS, order+1, cv_mem->cv_cvals, ca_mem->ca_YS, yS);
3024 if (retval != CV_SUCCESS) return (CV_VECTOROP_ERR);
3025 }
3026
3027 return(CV_SUCCESS);
3028
3029 }
3030
3031 /*
3032 * =================================================================
3033 * WRAPPERS FOR ADJOINT SYSTEM
3034 * =================================================================
3035 */
3036 /*
3037 * CVArhs
3038 *
3039 * This routine interfaces to the CVRhsFnB (or CVRhsFnBS) routine
3040 * provided by the user.
3041 */
3042
CVArhs(realtype t,N_Vector yB,N_Vector yBdot,void * cvode_mem)3043 static int CVArhs(realtype t, N_Vector yB,
3044 N_Vector yBdot, void *cvode_mem)
3045 {
3046 CVodeMem cv_mem;
3047 CVadjMem ca_mem;
3048 CVodeBMem cvB_mem;
3049 int flag, retval;
3050
3051 cv_mem = (CVodeMem) cvode_mem;
3052
3053 ca_mem = cv_mem->cv_adj_mem;
3054
3055 cvB_mem = ca_mem->ca_bckpbCrt;
3056
3057 /* Get forward solution from interpolation */
3058
3059 if (ca_mem->ca_IMinterpSensi)
3060 flag = ca_mem->ca_IMget(cv_mem, t, ca_mem->ca_ytmp, ca_mem->ca_yStmp);
3061 else
3062 flag = ca_mem->ca_IMget(cv_mem, t, ca_mem->ca_ytmp, NULL);
3063
3064 if (flag != CV_SUCCESS) {
3065 cvProcessError(cv_mem, -1, "CVODEA", "CVArhs", MSGCV_BAD_TINTERP, t);
3066 return(-1);
3067 }
3068
3069 /* Call the user's RHS function */
3070
3071 if (cvB_mem->cv_f_withSensi)
3072 retval = (cvB_mem->cv_fs)(t, ca_mem->ca_ytmp, ca_mem->ca_yStmp, yB, yBdot, cvB_mem->cv_user_data);
3073 else
3074 retval = (cvB_mem->cv_f)(t, ca_mem->ca_ytmp, yB, yBdot, cvB_mem->cv_user_data);
3075
3076 return(retval);
3077 }
3078
3079 /*
3080 * CVArhsQ
3081 *
3082 * This routine interfaces to the CVQuadRhsFnB (or CVQuadRhsFnBS) routine
3083 * provided by the user.
3084 */
3085
CVArhsQ(realtype t,N_Vector yB,N_Vector qBdot,void * cvode_mem)3086 static int CVArhsQ(realtype t, N_Vector yB,
3087 N_Vector qBdot, void *cvode_mem)
3088 {
3089 CVodeMem cv_mem;
3090 CVadjMem ca_mem;
3091 CVodeBMem cvB_mem;
3092 /* int flag; */
3093 int retval;
3094
3095 cv_mem = (CVodeMem) cvode_mem;
3096
3097 ca_mem = cv_mem->cv_adj_mem;
3098
3099 cvB_mem = ca_mem->ca_bckpbCrt;
3100
3101 /* Get forward solution from interpolation */
3102
3103 if (ca_mem->ca_IMinterpSensi)
3104 /* flag = */ ca_mem->ca_IMget(cv_mem, t, ca_mem->ca_ytmp, ca_mem->ca_yStmp);
3105 else
3106 /* flag = */ ca_mem->ca_IMget(cv_mem, t, ca_mem->ca_ytmp, NULL);
3107
3108 /* Call the user's RHS function */
3109
3110 if (cvB_mem->cv_fQ_withSensi)
3111 retval = (cvB_mem->cv_fQs)(t, ca_mem->ca_ytmp, ca_mem->ca_yStmp, yB, qBdot, cvB_mem->cv_user_data);
3112 else
3113 retval = (cvB_mem->cv_fQ)(t, ca_mem->ca_ytmp, yB, qBdot, cvB_mem->cv_user_data);
3114
3115 return(retval);
3116 }
3117