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