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 IDAA adjoint integrator.
16 * -----------------------------------------------------------------
17 */
18
19 /*=================================================================*/
20 /* Import Header Files */
21 /*=================================================================*/
22
23 #include <stdio.h>
24 #include <stdlib.h>
25
26 #include "idas_impl.h"
27 #include <sundials/sundials_math.h>
28
29 /*=================================================================*/
30 /* IDAA Private Constants */
31 /*=================================================================*/
32
33 #define ZERO RCONST(0.0) /* real 0.0 */
34 #define ONE RCONST(1.0) /* real 1.0 */
35 #define TWO RCONST(2.0) /* real 2.0 */
36 #define HUNDRED RCONST(100.0) /* real 100.0 */
37 #define FUZZ_FACTOR RCONST(1000000.0) /* fuzz factor for IDAAgetY */
38
39
40 /*=================================================================*/
41 /* Private Functions Prototypes */
42 /*=================================================================*/
43
44 static CkpntMem IDAAckpntInit(IDAMem IDA_mem);
45 static CkpntMem IDAAckpntNew(IDAMem IDA_mem);
46 static void IDAAckpntCopyVectors(IDAMem IDA_mem, CkpntMem ck_mem);
47 static booleantype IDAAckpntAllocVectors(IDAMem IDA_mem, CkpntMem ck_mem);
48 static void IDAAckpntDelete(CkpntMem *ck_memPtr);
49
50 static void IDAAbckpbDelete(IDABMem *IDAB_memPtr);
51
52 static booleantype IDAAdataMalloc(IDAMem IDA_mem);
53 static void IDAAdataFree(IDAMem IDA_mem);
54 static int IDAAdataStore(IDAMem IDA_mem, CkpntMem ck_mem);
55
56 static int IDAAckpntGet(IDAMem IDA_mem, CkpntMem ck_mem);
57
58 static booleantype IDAAhermiteMalloc(IDAMem IDA_mem);
59 static void IDAAhermiteFree(IDAMem IDA_mem);
60 static int IDAAhermiteStorePnt(IDAMem IDA_mem, DtpntMem d);
61 static int IDAAhermiteGetY(IDAMem IDA_mem, realtype t,
62 N_Vector yy, N_Vector yp,
63 N_Vector *yyS, N_Vector *ypS);
64
65 static booleantype IDAApolynomialMalloc(IDAMem IDA_mem);
66 static void IDAApolynomialFree(IDAMem IDA_mem);
67 static int IDAApolynomialStorePnt(IDAMem IDA_mem, DtpntMem d);
68 static int IDAApolynomialGetY(IDAMem IDA_mem, realtype t,
69 N_Vector yy, N_Vector yp,
70 N_Vector *yyS, N_Vector *ypS);
71
72 static int IDAAfindIndex(IDAMem ida_mem, realtype t,
73 long int *indx, booleantype *newpoint);
74
75 static int IDAAres(realtype tt,
76 N_Vector yyB, N_Vector ypB,
77 N_Vector resvalB, void *ida_mem);
78
79 static int IDAArhsQ(realtype tt,
80 N_Vector yyB, N_Vector ypB,
81 N_Vector rrQB, void *ida_mem);
82
83 static int IDAAGettnSolutionYp(IDAMem IDA_mem, N_Vector yp);
84 static int IDAAGettnSolutionYpS(IDAMem IDA_mem, N_Vector *ypS);
85
86 extern int IDAGetSolution(void *ida_mem, realtype t, N_Vector yret, N_Vector ypret);
87
88
89 /*=================================================================*/
90 /* Exported Functions */
91 /*=================================================================*/
92
93 /*
94 * IDAAdjInit
95 *
96 * This routine allocates space for the global IDAA memory
97 * structure.
98 */
99
100
IDAAdjInit(void * ida_mem,long int steps,int interp)101 int IDAAdjInit(void *ida_mem, long int steps, int interp)
102 {
103 IDAadjMem IDAADJ_mem;
104 IDAMem IDA_mem;
105
106 /* Check arguments */
107
108 if (ida_mem == NULL) {
109 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAAdjInit", MSGAM_NULL_IDAMEM);
110 return(IDA_MEM_NULL);
111 }
112 IDA_mem = (IDAMem)ida_mem;
113
114 if (steps <= 0) {
115 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAAdjInit", MSGAM_BAD_STEPS);
116 return(IDA_ILL_INPUT);
117 }
118
119 if ( (interp != IDA_HERMITE) && (interp != IDA_POLYNOMIAL) ) {
120 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAAdjInit", MSGAM_BAD_INTERP);
121 return(IDA_ILL_INPUT);
122 }
123
124 /* Allocate memory block for IDAadjMem. */
125 IDAADJ_mem = (IDAadjMem) malloc(sizeof(struct IDAadjMemRec));
126 if (IDAADJ_mem == NULL) {
127 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDAA", "IDAAdjInit", MSGAM_MEM_FAIL);
128 return(IDA_MEM_FAIL);
129 }
130
131 /* Attach IDAS memory for forward runs */
132 IDA_mem->ida_adj_mem = IDAADJ_mem;
133
134 /* Initialization of check points. */
135 IDAADJ_mem->ck_mem = NULL;
136 IDAADJ_mem->ia_nckpnts = 0;
137 IDAADJ_mem->ia_ckpntData = NULL;
138
139
140 /* Initialization of interpolation data. */
141 IDAADJ_mem->ia_interpType = interp;
142 IDAADJ_mem->ia_nsteps = steps;
143
144 /* Last index used in IDAAfindIndex, initailize to invalid value */
145 IDAADJ_mem->ia_ilast = -1;
146
147 /* Allocate space for the array of Data Point structures. */
148 if (IDAAdataMalloc(IDA_mem) == SUNFALSE) {
149 free(IDAADJ_mem); IDAADJ_mem = NULL;
150 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDAA", "IDAAdjInit", MSGAM_MEM_FAIL);
151 return(IDA_MEM_FAIL);
152 }
153
154 /* Attach functions for the appropriate interpolation module */
155 switch(interp) {
156
157 case IDA_HERMITE:
158 IDAADJ_mem->ia_malloc = IDAAhermiteMalloc;
159 IDAADJ_mem->ia_free = IDAAhermiteFree;
160 IDAADJ_mem->ia_getY = IDAAhermiteGetY;
161 IDAADJ_mem->ia_storePnt = IDAAhermiteStorePnt;
162 break;
163
164 case IDA_POLYNOMIAL:
165
166 IDAADJ_mem->ia_malloc = IDAApolynomialMalloc;
167 IDAADJ_mem->ia_free = IDAApolynomialFree;
168 IDAADJ_mem->ia_getY = IDAApolynomialGetY;
169 IDAADJ_mem->ia_storePnt = IDAApolynomialStorePnt;
170 break;
171 }
172
173 /* The interpolation module has not been initialized yet */
174 IDAADJ_mem->ia_mallocDone = SUNFALSE;
175
176 /* By default we will store but not interpolate sensitivities
177 * - storeSensi will be set in IDASolveF to SUNFALSE if FSA is not enabled
178 * or if the user forced this through IDAAdjSetNoSensi
179 * - interpSensi will be set in IDASolveB to SUNTRUE if storeSensi is SUNTRUE
180 * and if at least one backward problem requires sensitivities
181 * - noInterp will be set in IDACalcICB to SUNTRUE before the call to
182 * IDACalcIC and SUNFALSE after.*/
183
184 IDAADJ_mem->ia_storeSensi = SUNTRUE;
185 IDAADJ_mem->ia_interpSensi = SUNFALSE;
186 IDAADJ_mem->ia_noInterp = SUNFALSE;
187
188 /* Initialize backward problems. */
189 IDAADJ_mem->IDAB_mem = NULL;
190 IDAADJ_mem->ia_bckpbCrt = NULL;
191 IDAADJ_mem->ia_nbckpbs = 0;
192
193 /* IDASolveF and IDASolveB not called yet. */
194 IDAADJ_mem->ia_firstIDAFcall = SUNTRUE;
195 IDAADJ_mem->ia_tstopIDAFcall = SUNFALSE;
196
197 IDAADJ_mem->ia_firstIDABcall = SUNTRUE;
198
199 IDAADJ_mem->ia_rootret = SUNFALSE;
200
201 /* Adjoint module initialized and allocated. */
202 IDA_mem->ida_adj = SUNTRUE;
203 IDA_mem->ida_adjMallocDone = SUNTRUE;
204
205 return(IDA_SUCCESS);
206 }
207
208 /*
209 * IDAAdjReInit
210 *
211 * IDAAdjReInit reinitializes the IDAS memory structure for ASA
212 */
213
IDAAdjReInit(void * ida_mem)214 int IDAAdjReInit(void *ida_mem)
215 {
216 IDAadjMem IDAADJ_mem;
217 IDAMem IDA_mem;
218
219 /* Check arguments */
220
221 if (ida_mem == NULL) {
222 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAAdjReInit", MSGAM_NULL_IDAMEM);
223 return(IDA_MEM_NULL);
224 }
225 IDA_mem = (IDAMem)ida_mem;
226
227 /* Was ASA previously initialized? */
228 if(IDA_mem->ida_adjMallocDone == SUNFALSE) {
229 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAAdjReInit", MSGAM_NO_ADJ);
230 return(IDA_NO_ADJ);
231 }
232
233 IDAADJ_mem = IDA_mem->ida_adj_mem;
234
235 /* Free all stored checkpoints. */
236 while (IDAADJ_mem->ck_mem != NULL)
237 IDAAckpntDelete(&(IDAADJ_mem->ck_mem));
238
239 IDAADJ_mem->ck_mem = NULL;
240 IDAADJ_mem->ia_nckpnts = 0;
241 IDAADJ_mem->ia_ckpntData = NULL;
242
243 /* Flags for tracking the first calls to IDASolveF and IDASolveF. */
244 IDAADJ_mem->ia_firstIDAFcall = SUNTRUE;
245 IDAADJ_mem->ia_tstopIDAFcall = SUNFALSE;
246 IDAADJ_mem->ia_firstIDABcall = SUNTRUE;
247
248 return(IDA_SUCCESS);
249 }
250
251 /*
252 * IDAAdjFree
253 *
254 * IDAAdjFree routine frees the memory allocated by IDAAdjInit.
255 */
256
257
IDAAdjFree(void * ida_mem)258 void IDAAdjFree(void *ida_mem)
259 {
260 IDAMem IDA_mem;
261 IDAadjMem IDAADJ_mem;
262
263 if (ida_mem == NULL) return;
264 IDA_mem = (IDAMem) ida_mem;
265
266 if(IDA_mem->ida_adjMallocDone) {
267
268 /* Data for adjoint. */
269 IDAADJ_mem = IDA_mem->ida_adj_mem;
270
271 /* Delete check points one by one */
272 while (IDAADJ_mem->ck_mem != NULL) {
273 IDAAckpntDelete(&(IDAADJ_mem->ck_mem));
274 }
275
276 IDAAdataFree(IDA_mem);
277
278 /* Free all backward problems. */
279 while (IDAADJ_mem->IDAB_mem != NULL)
280 IDAAbckpbDelete( &(IDAADJ_mem->IDAB_mem) );
281
282 /* Free IDAA memory. */
283 free(IDAADJ_mem);
284
285 IDA_mem->ida_adj_mem = NULL;
286 }
287 }
288
289 /*
290 * =================================================================
291 * PRIVATE FUNCTIONS FOR BACKWARD PROBLEMS
292 * =================================================================
293 */
294
IDAAbckpbDelete(IDABMem * IDAB_memPtr)295 static void IDAAbckpbDelete(IDABMem *IDAB_memPtr)
296 {
297 IDABMem IDAB_mem = (*IDAB_memPtr);
298 void * ida_mem;
299
300 if (IDAB_mem == NULL) return;
301
302 /* Move head to the next element in list. */
303 *IDAB_memPtr = IDAB_mem->ida_next;
304
305 /* IDAB_mem is going to be deallocated. */
306
307 /* Free IDAS memory for this backward problem. */
308 ida_mem = (void *)IDAB_mem->IDA_mem;
309 IDAFree(&ida_mem);
310
311 /* Free linear solver memory. */
312 if (IDAB_mem->ida_lfree != NULL) IDAB_mem->ida_lfree(IDAB_mem);
313
314 /* Free preconditioner memory. */
315 if (IDAB_mem->ida_pfree != NULL) IDAB_mem->ida_pfree(IDAB_mem);
316
317 /* Free any workspace vectors. */
318 N_VDestroy(IDAB_mem->ida_yy);
319 N_VDestroy(IDAB_mem->ida_yp);
320
321 /* Free the node itself. */
322 free(IDAB_mem);
323 IDAB_mem = NULL;
324 }
325
326 /*=================================================================*/
327 /* Wrappers for IDAA */
328 /*=================================================================*/
329
330 /*
331 * IDASolveF
332 *
333 * This routine integrates to tout and returns solution into yout.
334 * In the same time, it stores check point data every 'steps' steps.
335 *
336 * IDASolveF can be called repeatedly by the user. The last tout
337 * will be used as the starting time for the backward integration.
338 *
339 * ncheckPtr points to the number of check points stored so far.
340 */
341
IDASolveF(void * ida_mem,realtype tout,realtype * tret,N_Vector yret,N_Vector ypret,int itask,int * ncheckPtr)342 int IDASolveF(void *ida_mem, realtype tout, realtype *tret,
343 N_Vector yret, N_Vector ypret, int itask, int *ncheckPtr)
344 {
345 IDAadjMem IDAADJ_mem;
346 IDAMem IDA_mem;
347 CkpntMem tmp;
348 DtpntMem *dt_mem;
349 long int nstloc;
350 int flag, i;
351 booleantype allocOK, earlyret;
352 realtype ttest;
353
354 /* Is the mem OK? */
355 if (ida_mem == NULL) {
356 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDASolveF", MSGAM_NULL_IDAMEM);
357 return(IDA_MEM_NULL);
358 }
359 IDA_mem = (IDAMem) ida_mem;
360
361 /* Is ASA initialized ? */
362 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
363 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDASolveF", MSGAM_NO_ADJ);
364 return(IDA_NO_ADJ);
365 }
366 IDAADJ_mem = IDA_mem->ida_adj_mem;
367
368 /* Check for yret != NULL */
369 if (yret == NULL) {
370 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDASolveF", MSG_YRET_NULL);
371 return(IDA_ILL_INPUT);
372 }
373
374 /* Check for ypret != NULL */
375 if (ypret == NULL) {
376 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDASolveF", MSG_YPRET_NULL);
377 return(IDA_ILL_INPUT);
378 }
379 /* Check for tret != NULL */
380 if (tret == NULL) {
381 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDASolveF", MSG_TRET_NULL);
382 return(IDA_ILL_INPUT);
383 }
384
385 /* Check for valid itask */
386 if ( (itask != IDA_NORMAL) && (itask != IDA_ONE_STEP) ) {
387 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDASolveF", MSG_BAD_ITASK);
388 return(IDA_ILL_INPUT);
389 }
390
391 /* All memory checks done, proceed ... */
392
393 dt_mem = IDAADJ_mem->dt_mem;
394
395 /* If tstop is enabled, store some info */
396 if (IDA_mem->ida_tstopset) {
397 IDAADJ_mem->ia_tstopIDAFcall = SUNTRUE;
398 IDAADJ_mem->ia_tstopIDAF = IDA_mem->ida_tstop;
399 }
400
401 /* On the first step:
402 * - set tinitial
403 * - initialize list of check points
404 * - if needed, initialize the interpolation module
405 * - load dt_mem[0]
406 * On subsequent steps, test if taking a new step is necessary.
407 */
408 if ( IDAADJ_mem->ia_firstIDAFcall ) {
409
410 IDAADJ_mem->ia_tinitial = IDA_mem->ida_tn;
411 IDAADJ_mem->ck_mem = IDAAckpntInit(IDA_mem);
412 if (IDAADJ_mem->ck_mem == NULL) {
413 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDAA", "IDASolveF", MSG_MEM_FAIL);
414 return(IDA_MEM_FAIL);
415 }
416
417 if (!IDAADJ_mem->ia_mallocDone) {
418 /* Do we need to store sensitivities? */
419 if (!IDA_mem->ida_sensi) IDAADJ_mem->ia_storeSensi = SUNFALSE;
420
421 /* Allocate space for interpolation data */
422 allocOK = IDAADJ_mem->ia_malloc(IDA_mem);
423 if (!allocOK) {
424 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDAA", "IDASolveF", MSG_MEM_FAIL);
425 return(IDA_MEM_FAIL);
426 }
427
428 /* Rename phi and, if needed, phiS for use in interpolation */
429 for (i=0;i<MXORDP1;i++) IDAADJ_mem->ia_Y[i] = IDA_mem->ida_phi[i];
430 if (IDAADJ_mem->ia_storeSensi) {
431 for (i=0;i<MXORDP1;i++)
432 IDAADJ_mem->ia_YS[i] = IDA_mem->ida_phiS[i];
433 }
434
435 IDAADJ_mem->ia_mallocDone = SUNTRUE;
436 }
437
438 dt_mem[0]->t = IDAADJ_mem->ck_mem->ck_t0;
439 IDAADJ_mem->ia_storePnt(IDA_mem, dt_mem[0]);
440
441 IDAADJ_mem->ia_firstIDAFcall = SUNFALSE;
442
443 } else if ( itask == IDA_NORMAL ) {
444
445 /* When in normal mode, check if tout was passed or if a previous root was
446 not reported and return an interpolated solution. No changes to ck_mem
447 or dt_mem are needed. */
448
449 /* flag to signal if an early return is needed */
450 earlyret = SUNFALSE;
451
452 /* if a root needs to be reported compare tout to troot otherwise compare
453 to the current time tn */
454 ttest = (IDAADJ_mem->ia_rootret) ? IDAADJ_mem->ia_troot : IDA_mem->ida_tn;
455
456 if ((ttest - tout)*IDA_mem->ida_hh >= ZERO) {
457 /* ttest is after tout, interpolate to tout */
458 *tret = tout;
459 flag = IDAGetSolution(IDA_mem, tout, yret, ypret);
460 earlyret = SUNTRUE;
461 } else if (IDAADJ_mem->ia_rootret) {
462 /* tout is after troot, interpolate to troot */
463 *tret = IDAADJ_mem->ia_troot;
464 flag = IDAGetSolution(IDA_mem, IDAADJ_mem->ia_troot, yret, ypret);
465 flag = IDA_ROOT_RETURN;
466 IDAADJ_mem->ia_rootret = SUNFALSE;
467 earlyret = SUNTRUE;
468 }
469
470 /* return if necessary */
471 if (earlyret) {
472 *ncheckPtr = IDAADJ_mem->ia_nckpnts;
473 IDAADJ_mem->ia_newData = SUNTRUE;
474 IDAADJ_mem->ia_ckpntData = IDAADJ_mem->ck_mem;
475 IDAADJ_mem->ia_np = IDA_mem->ida_nst % IDAADJ_mem->ia_nsteps + 1;
476 return(flag);
477 }
478
479 }
480
481 /* Integrate to tout (in IDA_ONE_STEP mode) while loading check points */
482 nstloc = 0;
483 for(;;) {
484
485 /* Check for too many steps */
486
487 if ( (IDA_mem->ida_mxstep>0) && (nstloc >= IDA_mem->ida_mxstep) ) {
488 IDAProcessError(IDA_mem, IDA_TOO_MUCH_WORK, "IDAA", "IDASolveF",
489 MSG_MAX_STEPS, IDA_mem->ida_tn);
490 flag = IDA_TOO_MUCH_WORK;
491 break;
492 }
493
494 /* Perform one step of the integration */
495
496 flag = IDASolve(IDA_mem, tout, tret, yret, ypret, IDA_ONE_STEP);
497 if (flag < 0) break;
498
499 nstloc++;
500
501 /* Test if a new check point is needed */
502
503 if ( IDA_mem->ida_nst % IDAADJ_mem->ia_nsteps == 0 ) {
504
505 IDAADJ_mem->ck_mem->ck_t1 = IDA_mem->ida_tn;
506
507 /* Create a new check point, load it, and append it to the list */
508 tmp = IDAAckpntNew(IDA_mem);
509 if (tmp == NULL) {
510 flag = IDA_MEM_FAIL;
511 break;
512 }
513
514 tmp->ck_next = IDAADJ_mem->ck_mem;
515 IDAADJ_mem->ck_mem = tmp;
516 IDAADJ_mem->ia_nckpnts++;
517
518 IDA_mem->ida_forceSetup = SUNTRUE;
519
520 /* Reset i=0 and load dt_mem[0] */
521 dt_mem[0]->t = IDAADJ_mem->ck_mem->ck_t0;
522 IDAADJ_mem->ia_storePnt(IDA_mem, dt_mem[0]);
523
524 } else {
525
526 /* Load next point in dt_mem */
527 dt_mem[IDA_mem->ida_nst%IDAADJ_mem->ia_nsteps]->t = IDA_mem->ida_tn;
528 IDAADJ_mem->ia_storePnt(IDA_mem, dt_mem[IDA_mem->ida_nst % IDAADJ_mem->ia_nsteps]);
529
530 }
531
532 /* Set t1 field of the current ckeck point structure
533 for the case in which there will be no future
534 check points */
535 IDAADJ_mem->ck_mem->ck_t1 = IDA_mem->ida_tn;
536
537 /* tfinal is now set to tn */
538 IDAADJ_mem->ia_tfinal = IDA_mem->ida_tn;
539
540 /* Return if in IDA_ONE_STEP mode */
541 if (itask == IDA_ONE_STEP) break;
542
543 /* IDA_NORMAL_STEP returns */
544
545 /* Return if tout reached */
546 if ( (*tret - tout)*IDA_mem->ida_hh >= ZERO ) {
547
548 /* If this was a root return, save the root time to return later */
549 if (flag == IDA_ROOT_RETURN) {
550 IDAADJ_mem->ia_rootret = SUNTRUE;
551 IDAADJ_mem->ia_troot = *tret;
552 }
553
554 /* Get solution value at tout to return now */
555 *tret = tout;
556 flag = IDAGetSolution(IDA_mem, tout, yret, ypret);
557
558 /* Reset tretlast in IDA_mem so that IDAGetQuad and IDAGetSens
559 * evaluate quadratures and/or sensitivities at the proper time */
560 IDA_mem->ida_tretlast = tout;
561
562 break;
563 }
564
565 /* Return if tstop or a root was found */
566 if ((flag == IDA_TSTOP_RETURN) || (flag == IDA_ROOT_RETURN)) break;
567
568 } /* end of for(;;) */
569
570 /* Get ncheck from IDAADJ_mem */
571 *ncheckPtr = IDAADJ_mem->ia_nckpnts;
572
573 /* Data is available for the last interval */
574 IDAADJ_mem->ia_newData = SUNTRUE;
575 IDAADJ_mem->ia_ckpntData = IDAADJ_mem->ck_mem;
576 IDAADJ_mem->ia_np = IDA_mem->ida_nst % IDAADJ_mem->ia_nsteps + 1;
577
578 return(flag);
579 }
580
581
582
583
584 /*
585 * =================================================================
586 * FUNCTIONS FOR BACKWARD PROBLEMS
587 * =================================================================
588 */
589
IDACreateB(void * ida_mem,int * which)590 int IDACreateB(void *ida_mem, int *which)
591 {
592 IDAMem IDA_mem;
593 void* ida_memB;
594 IDABMem new_IDAB_mem;
595 IDAadjMem IDAADJ_mem;
596
597 /* Is the mem OK? */
598 if (ida_mem == NULL) {
599 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDACreateB", MSGAM_NULL_IDAMEM);
600 return(IDA_MEM_NULL);
601 }
602 IDA_mem = (IDAMem) ida_mem;
603
604 /* Is ASA initialized ? */
605 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
606 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDACreateB", MSGAM_NO_ADJ);
607 return(IDA_NO_ADJ);
608 }
609 IDAADJ_mem = IDA_mem->ida_adj_mem;
610
611 /* Allocate a new IDABMem struct. */
612 new_IDAB_mem = (IDABMem) malloc( sizeof( struct IDABMemRec ) );
613 if (new_IDAB_mem == NULL) {
614 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDAA", "IDACreateB", MSG_MEM_FAIL);
615 return(IDA_MEM_FAIL);
616 }
617
618 /* Allocate the IDAMem struct needed by this backward problem. */
619 ida_memB = IDACreate();
620 if (ida_memB == NULL) {
621 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDAA", "IDACreateB", MSG_MEM_FAIL);
622 return(IDA_MEM_FAIL);
623 }
624
625 /* Save ida_mem in ida_memB as user data. */
626 IDASetUserData(ida_memB, ida_mem);
627
628 /* Set same error output and handler for ida_memB. */
629 IDASetErrHandlerFn(ida_memB, IDA_mem->ida_ehfun, IDA_mem->ida_eh_data);
630 IDASetErrFile(ida_memB, IDA_mem->ida_errfp);
631
632 /* Initialize fields in the IDABMem struct. */
633 new_IDAB_mem->ida_index = IDAADJ_mem->ia_nbckpbs;
634 new_IDAB_mem->IDA_mem = (IDAMem) ida_memB;
635
636 new_IDAB_mem->ida_res = NULL;
637 new_IDAB_mem->ida_resS = NULL;
638 new_IDAB_mem->ida_rhsQ = NULL;
639 new_IDAB_mem->ida_rhsQS = NULL;
640
641
642 new_IDAB_mem->ida_user_data = NULL;
643
644 new_IDAB_mem->ida_lmem = NULL;
645 new_IDAB_mem->ida_lfree = NULL;
646 new_IDAB_mem->ida_pmem = NULL;
647 new_IDAB_mem->ida_pfree = NULL;
648
649 new_IDAB_mem->ida_yy = NULL;
650 new_IDAB_mem->ida_yp = NULL;
651
652 new_IDAB_mem->ida_res_withSensi = SUNFALSE;
653 new_IDAB_mem->ida_rhsQ_withSensi = SUNFALSE;
654
655 /* Attach the new object to the beginning of the linked list IDAADJ_mem->IDAB_mem. */
656 new_IDAB_mem->ida_next = IDAADJ_mem->IDAB_mem;
657 IDAADJ_mem->IDAB_mem = new_IDAB_mem;
658
659 /* Return the assigned index. This id is used as identificator and has to be passed
660 to IDAInitB and other ***B functions that set the optional inputs for this
661 backward problem. */
662 *which = IDAADJ_mem->ia_nbckpbs;
663
664 /*Increase the counter of the backward problems stored. */
665 IDAADJ_mem->ia_nbckpbs++;
666
667 return(IDA_SUCCESS);
668
669 }
670
IDAInitB(void * ida_mem,int which,IDAResFnB resB,realtype tB0,N_Vector yyB0,N_Vector ypB0)671 int IDAInitB(void *ida_mem, int which, IDAResFnB resB,
672 realtype tB0, N_Vector yyB0, N_Vector ypB0)
673 {
674 IDAadjMem IDAADJ_mem;
675 IDAMem IDA_mem;
676 IDABMem IDAB_mem;
677 void * ida_memB;
678 int flag;
679
680 if (ida_mem == NULL) {
681 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAInitB", MSGAM_NULL_IDAMEM);
682 return(IDA_MEM_NULL);
683 }
684 IDA_mem = (IDAMem) ida_mem;
685
686 /* Is ASA initialized ? */
687 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
688 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAInitB", MSGAM_NO_ADJ);
689 return(IDA_NO_ADJ);
690 }
691 IDAADJ_mem = IDA_mem->ida_adj_mem;
692
693 /* Check the initial time for this backward problem against the adjoint data. */
694 if ( (tB0 < IDAADJ_mem->ia_tinitial) || (tB0 > IDAADJ_mem->ia_tfinal) ) {
695 IDAProcessError(IDA_mem, IDA_BAD_TB0, "IDAA", "IDAInitB", MSGAM_BAD_TB0);
696 return(IDA_BAD_TB0);
697 }
698
699 /* Check the value of which */
700 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
701 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAInitB", MSGAM_BAD_WHICH);
702 return(IDA_ILL_INPUT);
703 }
704
705 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
706 IDAB_mem = IDAADJ_mem->IDAB_mem;
707 while (IDAB_mem != NULL) {
708 if( which == IDAB_mem->ida_index ) break;
709 /* advance */
710 IDAB_mem = IDAB_mem->ida_next;
711 }
712
713 /* Get the IDAMem corresponding to this backward problem. */
714 ida_memB = (void*) IDAB_mem->IDA_mem;
715
716 /* Call the IDAInit for this backward problem. */
717 flag = IDAInit(ida_memB, IDAAres, tB0, yyB0, ypB0);
718 if (IDA_SUCCESS != flag) return(flag);
719
720 /* Copy residual function in IDAB_mem. */
721 IDAB_mem->ida_res = resB;
722 IDAB_mem->ida_res_withSensi = SUNFALSE;
723
724 /* Initialized the initial time field. */
725 IDAB_mem->ida_t0 = tB0;
726
727 /* Allocate and initialize space workspace vectors. */
728 IDAB_mem->ida_yy = N_VClone(yyB0);
729 IDAB_mem->ida_yp = N_VClone(yyB0);
730 N_VScale(ONE, yyB0, IDAB_mem->ida_yy);
731 N_VScale(ONE, ypB0, IDAB_mem->ida_yp);
732
733 return(flag);
734
735 }
736
IDAInitBS(void * ida_mem,int which,IDAResFnBS resS,realtype tB0,N_Vector yyB0,N_Vector ypB0)737 int IDAInitBS(void *ida_mem, int which, IDAResFnBS resS,
738 realtype tB0, N_Vector yyB0, N_Vector ypB0)
739 {
740 IDAadjMem IDAADJ_mem;
741 IDAMem IDA_mem;
742 IDABMem IDAB_mem;
743 void * ida_memB;
744 int flag;
745
746 if (ida_mem == NULL) {
747 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAInitBS", MSGAM_NULL_IDAMEM);
748 return(IDA_MEM_NULL);
749 }
750 IDA_mem = (IDAMem) ida_mem;
751
752 /* Is ASA initialized ? */
753 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
754 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAInitBS", MSGAM_NO_ADJ);
755 return(IDA_NO_ADJ);
756 }
757 IDAADJ_mem = IDA_mem->ida_adj_mem;
758
759 /* Check the initial time for this backward problem against the adjoint data. */
760 if ( (tB0 < IDAADJ_mem->ia_tinitial) || (tB0 > IDAADJ_mem->ia_tfinal) ) {
761 IDAProcessError(IDA_mem, IDA_BAD_TB0, "IDAA", "IDAInitBS", MSGAM_BAD_TB0);
762 return(IDA_BAD_TB0);
763 }
764
765 /* Were sensitivities active during the forward integration? */
766 if (!IDAADJ_mem->ia_storeSensi) {
767 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAInitBS", MSGAM_BAD_SENSI);
768 return(IDA_ILL_INPUT);
769 }
770
771 /* Check the value of which */
772 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
773 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAInitBS", MSGAM_BAD_WHICH);
774 return(IDA_ILL_INPUT);
775 }
776
777 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
778 IDAB_mem = IDAADJ_mem->IDAB_mem;
779 while (IDAB_mem != NULL) {
780 if( which == IDAB_mem->ida_index ) break;
781 /* advance */
782 IDAB_mem = IDAB_mem->ida_next;
783 }
784
785 /* Get the IDAMem corresponding to this backward problem. */
786 ida_memB = (void*) IDAB_mem->IDA_mem;
787
788 /* Allocate and set the IDAS object */
789 flag = IDAInit(ida_memB, IDAAres, tB0, yyB0, ypB0);
790
791 if (flag != IDA_SUCCESS) return(flag);
792
793 /* Copy residual function pointer in IDAB_mem. */
794 IDAB_mem->ida_res_withSensi = SUNTRUE;
795 IDAB_mem->ida_resS = resS;
796
797 /* Allocate space and initialize the yy and yp vectors. */
798 IDAB_mem->ida_t0 = tB0;
799 IDAB_mem->ida_yy = N_VClone(yyB0);
800 IDAB_mem->ida_yp = N_VClone(ypB0);
801 N_VScale(ONE, yyB0, IDAB_mem->ida_yy);
802 N_VScale(ONE, ypB0, IDAB_mem->ida_yp);
803
804 return(IDA_SUCCESS);
805 }
806
807
IDAReInitB(void * ida_mem,int which,realtype tB0,N_Vector yyB0,N_Vector ypB0)808 int IDAReInitB(void *ida_mem, int which,
809 realtype tB0, N_Vector yyB0, N_Vector ypB0)
810 {
811
812 IDAadjMem IDAADJ_mem;
813 IDAMem IDA_mem;
814 IDABMem IDAB_mem;
815 void * ida_memB;
816 int flag;
817
818 if (ida_mem == NULL) {
819 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAReInitB", MSGAM_NULL_IDAMEM);
820 return(IDA_MEM_NULL);
821 }
822 IDA_mem = (IDAMem) ida_mem;
823
824 /* Is ASA initialized ? */
825 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
826 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAReInitB", MSGAM_NO_ADJ);
827 return(IDA_NO_ADJ);
828 }
829 IDAADJ_mem = IDA_mem->ida_adj_mem;
830
831 /* Check the initial time for this backward problem against the adjoint data. */
832 if ( (tB0 < IDAADJ_mem->ia_tinitial) || (tB0 > IDAADJ_mem->ia_tfinal) ) {
833 IDAProcessError(IDA_mem, IDA_BAD_TB0, "IDAA", "IDAReInitB", MSGAM_BAD_TB0);
834 return(IDA_BAD_TB0);
835 }
836
837 /* Check the value of which */
838 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
839 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAReInitB", MSGAM_BAD_WHICH);
840 return(IDA_ILL_INPUT);
841 }
842
843 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
844 IDAB_mem = IDAADJ_mem->IDAB_mem;
845 while (IDAB_mem != NULL) {
846 if( which == IDAB_mem->ida_index ) break;
847 /* advance */
848 IDAB_mem = IDAB_mem->ida_next;
849 }
850
851 /* Get the IDAMem corresponding to this backward problem. */
852 ida_memB = (void*) IDAB_mem->IDA_mem;
853
854
855 /* Call the IDAReInit for this backward problem. */
856 flag = IDAReInit(ida_memB, tB0, yyB0, ypB0);
857 return(flag);
858 }
859
IDASStolerancesB(void * ida_mem,int which,realtype relTolB,realtype absTolB)860 int IDASStolerancesB(void *ida_mem, int which,
861 realtype relTolB, realtype absTolB)
862 {
863 IDAMem IDA_mem;
864 IDAadjMem IDAADJ_mem;
865 IDABMem IDAB_mem;
866 void *ida_memB;
867
868 /* Is ida_mem valid? */
869 if (ida_mem == NULL) {
870 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDASStolerancesB", MSGAM_NULL_IDAMEM);
871 return IDA_MEM_NULL;
872 }
873 IDA_mem = (IDAMem) ida_mem;
874
875 /* Is ASA initialized? */
876 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
877 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDASStolerancesB", MSGAM_NO_ADJ);
878 return(IDA_NO_ADJ);
879 }
880 IDAADJ_mem = IDA_mem->ida_adj_mem;
881
882 /* Check the value of which */
883 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
884 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDASStolerancesB", MSGAM_BAD_WHICH);
885 return(IDA_ILL_INPUT);
886 }
887
888 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
889 IDAB_mem = IDAADJ_mem->IDAB_mem;
890 while (IDAB_mem != NULL) {
891 if( which == IDAB_mem->ida_index ) break;
892 /* advance */
893 IDAB_mem = IDAB_mem->ida_next;
894 }
895
896 /* Get the IDAMem corresponding to this backward problem. */
897 ida_memB = (void*) IDAB_mem->IDA_mem;
898
899 /* Set tolerances and return. */
900 return IDASStolerances(ida_memB, relTolB, absTolB);
901
902 }
IDASVtolerancesB(void * ida_mem,int which,realtype relTolB,N_Vector absTolB)903 int IDASVtolerancesB(void *ida_mem, int which,
904 realtype relTolB, N_Vector absTolB)
905 {
906 IDAMem IDA_mem;
907 IDAadjMem IDAADJ_mem;
908 IDABMem IDAB_mem;
909 void *ida_memB;
910
911 /* Is ida_mem valid? */
912 if (ida_mem == NULL) {
913 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDASVtolerancesB", MSGAM_NULL_IDAMEM);
914 return IDA_MEM_NULL;
915 }
916 IDA_mem = (IDAMem) ida_mem;
917
918 /* Is ASA initialized? */
919 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
920 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDASVtolerancesB", MSGAM_NO_ADJ);
921 return(IDA_NO_ADJ);
922 }
923 IDAADJ_mem = IDA_mem->ida_adj_mem;
924
925 /* Check the value of which */
926 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
927 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDASVtolerancesB", MSGAM_BAD_WHICH);
928 return(IDA_ILL_INPUT);
929 }
930
931 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
932 IDAB_mem = IDAADJ_mem->IDAB_mem;
933 while (IDAB_mem != NULL) {
934 if( which == IDAB_mem->ida_index ) break;
935 /* advance */
936 IDAB_mem = IDAB_mem->ida_next;
937 }
938
939 /* Get the IDAMem corresponding to this backward problem. */
940 ida_memB = (void*) IDAB_mem->IDA_mem;
941
942 /* Set tolerances and return. */
943 return IDASVtolerances(ida_memB, relTolB, absTolB);
944 }
945
IDAQuadSStolerancesB(void * ida_mem,int which,realtype reltolQB,realtype abstolQB)946 int IDAQuadSStolerancesB(void *ida_mem, int which,
947 realtype reltolQB, realtype abstolQB)
948 {
949 IDAMem IDA_mem;
950 IDAadjMem IDAADJ_mem;
951 IDABMem IDAB_mem;
952 void *ida_memB;
953
954 /* Is ida_mem valid? */
955 if (ida_mem == NULL) {
956 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAQuadSStolerancesB", MSGAM_NULL_IDAMEM);
957 return IDA_MEM_NULL;
958 }
959 IDA_mem = (IDAMem) ida_mem;
960
961 /* Is ASA initialized? */
962 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
963 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAQuadSStolerancesB", MSGAM_NO_ADJ);
964 return(IDA_NO_ADJ);
965 }
966 IDAADJ_mem = IDA_mem->ida_adj_mem;
967
968 /* Check the value of which */
969 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
970 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAQuadSStolerancesB", MSGAM_BAD_WHICH);
971 return(IDA_ILL_INPUT);
972 }
973
974 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
975 IDAB_mem = IDAADJ_mem->IDAB_mem;
976 while (IDAB_mem != NULL) {
977 if( which == IDAB_mem->ida_index ) break;
978 /* advance */
979 IDAB_mem = IDAB_mem->ida_next;
980 }
981 ida_memB = (void *) IDAB_mem->IDA_mem;
982
983 return IDAQuadSStolerances(ida_memB, reltolQB, abstolQB);
984 }
985
986
IDAQuadSVtolerancesB(void * ida_mem,int which,realtype reltolQB,N_Vector abstolQB)987 int IDAQuadSVtolerancesB(void *ida_mem, int which,
988 realtype reltolQB, N_Vector abstolQB)
989 {
990 IDAMem IDA_mem;
991 IDAadjMem IDAADJ_mem;
992 IDABMem IDAB_mem;
993 void *ida_memB;
994
995 /* Is ida_mem valid? */
996 if (ida_mem == NULL) {
997 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAQuadSVtolerancesB", MSGAM_NULL_IDAMEM);
998 return IDA_MEM_NULL;
999 }
1000 IDA_mem = (IDAMem) ida_mem;
1001
1002 /* Is ASA initialized? */
1003 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
1004 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAQuadSVtolerancesB", MSGAM_NO_ADJ);
1005 return(IDA_NO_ADJ);
1006 }
1007 IDAADJ_mem = IDA_mem->ida_adj_mem;
1008
1009 /* Check the value of which */
1010 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
1011 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAQuadSVtolerancesB", MSGAM_BAD_WHICH);
1012 return(IDA_ILL_INPUT);
1013 }
1014
1015 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
1016 IDAB_mem = IDAADJ_mem->IDAB_mem;
1017 while (IDAB_mem != NULL) {
1018 if( which == IDAB_mem->ida_index ) break;
1019 /* advance */
1020 IDAB_mem = IDAB_mem->ida_next;
1021 }
1022 ida_memB = (void *) IDAB_mem->IDA_mem;
1023
1024 return IDAQuadSVtolerances(ida_memB, reltolQB, abstolQB);
1025 }
1026
1027
IDAQuadInitB(void * ida_mem,int which,IDAQuadRhsFnB rhsQB,N_Vector yQB0)1028 int IDAQuadInitB(void *ida_mem, int which, IDAQuadRhsFnB rhsQB, N_Vector yQB0)
1029 {
1030 IDAMem IDA_mem;
1031 IDAadjMem IDAADJ_mem;
1032 IDABMem IDAB_mem;
1033 void *ida_memB;
1034 int flag;
1035
1036 /* Is ida_mem valid? */
1037 if (ida_mem == NULL) {
1038 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAQuadInitB", MSGAM_NULL_IDAMEM);
1039 return IDA_MEM_NULL;
1040 }
1041 IDA_mem = (IDAMem) ida_mem;
1042
1043 /* Is ASA initialized? */
1044 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
1045 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAQuadInitB", MSGAM_NO_ADJ);
1046 return(IDA_NO_ADJ);
1047 }
1048 IDAADJ_mem = IDA_mem->ida_adj_mem;
1049
1050 /* Check the value of which */
1051 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
1052 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAQuadInitB", MSGAM_BAD_WHICH);
1053 return(IDA_ILL_INPUT);
1054 }
1055
1056 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
1057 IDAB_mem = IDAADJ_mem->IDAB_mem;
1058 while (IDAB_mem != NULL) {
1059 if( which == IDAB_mem->ida_index ) break;
1060 /* advance */
1061 IDAB_mem = IDAB_mem->ida_next;
1062 }
1063 ida_memB = (void *) IDAB_mem->IDA_mem;
1064
1065 flag = IDAQuadInit(ida_memB, IDAArhsQ, yQB0);
1066 if (IDA_SUCCESS != flag) return flag;
1067
1068 IDAB_mem->ida_rhsQ_withSensi = SUNFALSE;
1069 IDAB_mem->ida_rhsQ = rhsQB;
1070
1071 return(flag);
1072 }
1073
1074
IDAQuadInitBS(void * ida_mem,int which,IDAQuadRhsFnBS rhsQS,N_Vector yQB0)1075 int IDAQuadInitBS(void *ida_mem, int which,
1076 IDAQuadRhsFnBS rhsQS, N_Vector yQB0)
1077 {
1078 IDAadjMem IDAADJ_mem;
1079 IDAMem IDA_mem;
1080 IDABMem IDAB_mem;
1081 void * ida_memB;
1082 int flag;
1083
1084 if (ida_mem == NULL) {
1085 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAQuadInitBS", MSGAM_NULL_IDAMEM);
1086 return(IDA_MEM_NULL);
1087 }
1088 IDA_mem = (IDAMem) ida_mem;
1089
1090 /* Is ASA initialized ? */
1091 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
1092 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAQuadInitBS", MSGAM_NO_ADJ);
1093 return(IDA_NO_ADJ);
1094 }
1095 IDAADJ_mem = IDA_mem->ida_adj_mem;
1096
1097 /* Check the value of which */
1098 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
1099 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAQuadInitBS", MSGAM_BAD_WHICH);
1100 return(IDA_ILL_INPUT);
1101 }
1102
1103 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
1104 IDAB_mem = IDAADJ_mem->IDAB_mem;
1105 while (IDAB_mem != NULL) {
1106 if( which == IDAB_mem->ida_index ) break;
1107 /* advance */
1108 IDAB_mem = IDAB_mem->ida_next;
1109 }
1110
1111 /* Get the IDAMem corresponding to this backward problem. */
1112 ida_memB = (void*) IDAB_mem->IDA_mem;
1113
1114 /* Allocate and set the IDAS object */
1115 flag = IDAQuadInit(ida_memB, IDAArhsQ, yQB0);
1116
1117 if (flag != IDA_SUCCESS) return(flag);
1118
1119 /* Copy RHS function pointer in IDAB_mem and enable quad sensitivities. */
1120 IDAB_mem->ida_rhsQ_withSensi = SUNTRUE;
1121 IDAB_mem->ida_rhsQS = rhsQS;
1122
1123 return(IDA_SUCCESS);
1124 }
1125
1126
IDAQuadReInitB(void * ida_mem,int which,N_Vector yQB0)1127 int IDAQuadReInitB(void *ida_mem, int which, N_Vector yQB0)
1128 {
1129 IDAMem IDA_mem;
1130 IDAadjMem IDAADJ_mem;
1131 IDABMem IDAB_mem;
1132 void *ida_memB;
1133
1134 /* Is ida_mem valid? */
1135 if (ida_mem == NULL) {
1136 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAQuadInitB", MSGAM_NULL_IDAMEM);
1137 return IDA_MEM_NULL;
1138 }
1139 IDA_mem = (IDAMem) ida_mem;
1140
1141 /* Is ASA initialized? */
1142 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
1143 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAQuadInitB", MSGAM_NO_ADJ);
1144 return(IDA_NO_ADJ);
1145 }
1146 IDAADJ_mem = IDA_mem->ida_adj_mem;
1147
1148 /* Check the value of which */
1149 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
1150 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAQuadInitB", MSGAM_BAD_WHICH);
1151 return(IDA_ILL_INPUT);
1152 }
1153
1154 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
1155 IDAB_mem = IDAADJ_mem->IDAB_mem;
1156 while (IDAB_mem != NULL) {
1157 if( which == IDAB_mem->ida_index ) break;
1158 /* advance */
1159 IDAB_mem = IDAB_mem->ida_next;
1160 }
1161 ida_memB = (void *) IDAB_mem->IDA_mem;
1162
1163 return IDAQuadReInit(ida_memB, yQB0);
1164 }
1165
1166
1167 /*
1168 * ----------------------------------------------------------------
1169 * Function : IDACalcICB
1170 * ----------------------------------------------------------------
1171 * IDACalcIC calculates corrected initial conditions for a DAE
1172 * backward system (index-one in semi-implicit form).
1173 * It uses Newton iteration combined with a Linesearch algorithm.
1174 * Calling IDACalcICB is optional. It is only necessary when the
1175 * initial conditions do not solve the given system. I.e., if
1176 * yB0 and ypB0 are known to satisfy the backward problem, then
1177 * a call to IDACalcIC is NOT necessary (for index-one problems).
1178 */
1179
IDACalcICB(void * ida_mem,int which,realtype tout1,N_Vector yy0,N_Vector yp0)1180 int IDACalcICB(void *ida_mem, int which, realtype tout1,
1181 N_Vector yy0, N_Vector yp0)
1182 {
1183 IDAMem IDA_mem;
1184 IDAadjMem IDAADJ_mem;
1185 IDABMem IDAB_mem;
1186 void *ida_memB;
1187 int flag;
1188
1189 /* Is ida_mem valid? */
1190 if (ida_mem == NULL) {
1191 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDACalcICB", MSGAM_NULL_IDAMEM);
1192 return IDA_MEM_NULL;
1193 }
1194 IDA_mem = (IDAMem) ida_mem;
1195
1196 /* Is ASA initialized? */
1197 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
1198 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDACalcICB", MSGAM_NO_ADJ);
1199 return(IDA_NO_ADJ);
1200 }
1201 IDAADJ_mem = IDA_mem->ida_adj_mem;
1202
1203 /* Check the value of which */
1204 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
1205 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDACalcICB", MSGAM_BAD_WHICH);
1206 return(IDA_ILL_INPUT);
1207 }
1208
1209 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
1210 IDAB_mem = IDAADJ_mem->IDAB_mem;
1211 while (IDAB_mem != NULL) {
1212 if( which == IDAB_mem->ida_index ) break;
1213 /* advance */
1214 IDAB_mem = IDAB_mem->ida_next;
1215 }
1216 ida_memB = (void *) IDAB_mem->IDA_mem;
1217
1218 /* The wrapper for user supplied res function requires ia_bckpbCrt from
1219 IDAAdjMem to be set to curent problem. */
1220 IDAADJ_mem->ia_bckpbCrt = IDAB_mem;
1221
1222 /* Save (y, y') in yyTmp and ypTmp for use in the res wrapper.*/
1223 /* yyTmp and ypTmp workspaces are safe to use if IDAADataStore is not called.*/
1224 N_VScale(ONE, yy0, IDAADJ_mem->ia_yyTmp);
1225 N_VScale(ONE, yp0, IDAADJ_mem->ia_ypTmp);
1226
1227 /* Set noInterp flag to SUNTRUE, so IDAARes will use user provided values for
1228 y and y' and will not call the interpolation routine(s). */
1229 IDAADJ_mem->ia_noInterp = SUNTRUE;
1230
1231 flag = IDACalcIC(ida_memB, IDA_YA_YDP_INIT, tout1);
1232
1233 /* Set interpolation on in IDAARes. */
1234 IDAADJ_mem->ia_noInterp = SUNFALSE;
1235
1236 return(flag);
1237 }
1238
1239 /*
1240 * ----------------------------------------------------------------
1241 * Function : IDACalcICBS
1242 * ----------------------------------------------------------------
1243 * IDACalcIC calculates corrected initial conditions for a DAE
1244 * backward system (index-one in semi-implicit form) that also
1245 * dependes on the sensivities.
1246 *
1247 * It calls IDACalcIC for the 'which' backward problem.
1248 */
1249
IDACalcICBS(void * ida_mem,int which,realtype tout1,N_Vector yy0,N_Vector yp0,N_Vector * yyS0,N_Vector * ypS0)1250 int IDACalcICBS(void *ida_mem, int which, realtype tout1,
1251 N_Vector yy0, N_Vector yp0,
1252 N_Vector *yyS0, N_Vector *ypS0)
1253 {
1254 IDAMem IDA_mem;
1255 IDAadjMem IDAADJ_mem;
1256 IDABMem IDAB_mem;
1257 void *ida_memB;
1258 int flag, is, retval;
1259
1260 /* Is ida_mem valid? */
1261 if (ida_mem == NULL) {
1262 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDACalcICBS", MSGAM_NULL_IDAMEM);
1263 return IDA_MEM_NULL;
1264 }
1265 IDA_mem = (IDAMem) ida_mem;
1266
1267 /* Is ASA initialized? */
1268 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
1269 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDACalcICBS", MSGAM_NO_ADJ);
1270 return(IDA_NO_ADJ);
1271 }
1272 IDAADJ_mem = IDA_mem->ida_adj_mem;
1273
1274 /* Were sensitivities active during the forward integration? */
1275 if (!IDAADJ_mem->ia_storeSensi) {
1276 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDACalcICBS", MSGAM_BAD_SENSI);
1277 return(IDA_ILL_INPUT);
1278 }
1279
1280 /* Check the value of which */
1281 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
1282 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDACalcICBS", MSGAM_BAD_WHICH);
1283 return(IDA_ILL_INPUT);
1284 }
1285
1286 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
1287 IDAB_mem = IDAADJ_mem->IDAB_mem;
1288 while (IDAB_mem != NULL) {
1289 if( which == IDAB_mem->ida_index ) break;
1290 /* advance */
1291 IDAB_mem = IDAB_mem->ida_next;
1292 }
1293 ida_memB = (void *) IDAB_mem->IDA_mem;
1294
1295 /* Was InitBS called for this problem? */
1296 if (!IDAB_mem->ida_res_withSensi) {
1297 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDACalcICBS", MSGAM_NO_INITBS);
1298 return(IDA_ILL_INPUT);
1299 }
1300
1301 /* The wrapper for user supplied res function requires ia_bckpbCrt from
1302 IDAAdjMem to be set to curent problem. */
1303 IDAADJ_mem->ia_bckpbCrt = IDAB_mem;
1304
1305 /* Save (y, y') and (y_p, y'_p) in yyTmp, ypTmp and yySTmp, ypSTmp.The wrapper
1306 for residual will use these values instead of calling interpolation routine.*/
1307
1308 /* The four workspaces variables are safe to use if IDAADataStore is not called.*/
1309 N_VScale(ONE, yy0, IDAADJ_mem->ia_yyTmp);
1310 N_VScale(ONE, yp0, IDAADJ_mem->ia_ypTmp);
1311
1312 for (is=0; is<IDA_mem->ida_Ns; is++)
1313 IDA_mem->ida_cvals[is] = ONE;
1314
1315 retval = N_VScaleVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_cvals,
1316 yyS0, IDAADJ_mem->ia_yySTmp);
1317 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
1318
1319 retval = N_VScaleVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_cvals,
1320 ypS0, IDAADJ_mem->ia_ypSTmp);
1321 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
1322
1323 /* Set noInterp flag to SUNTRUE, so IDAARes will use user provided values for
1324 y and y' and will not call the interpolation routine(s). */
1325 IDAADJ_mem->ia_noInterp = SUNTRUE;
1326
1327 flag = IDACalcIC(ida_memB, IDA_YA_YDP_INIT, tout1);
1328
1329 /* Set interpolation on in IDAARes. */
1330 IDAADJ_mem->ia_noInterp = SUNFALSE;
1331
1332 return(flag);
1333 }
1334
1335
1336 /*
1337 * IDASolveB
1338 *
1339 * This routine performs the backward integration from tB0
1340 * to tinitial through a sequence of forward-backward runs in
1341 * between consecutive check points. It returns the values of
1342 * the adjoint variables and any existing quadrature variables
1343 * at tinitial.
1344 *
1345 * On a successful return, IDASolveB returns IDA_SUCCESS.
1346 *
1347 * NOTE that IDASolveB DOES NOT return the solution for the
1348 * backward problem(s). Use IDAGetB to extract the solution
1349 * for any given backward problem.
1350 *
1351 * If there are multiple backward problems and multiple check points,
1352 * IDASolveB may not succeed in getting all problems to take one step
1353 * when called in ONE_STEP mode.
1354 */
1355
IDASolveB(void * ida_mem,realtype tBout,int itaskB)1356 int IDASolveB(void *ida_mem, realtype tBout, int itaskB)
1357 {
1358 IDAMem IDA_mem;
1359 IDAadjMem IDAADJ_mem;
1360 CkpntMem ck_mem;
1361 IDABMem IDAB_mem, tmp_IDAB_mem;
1362 int flag=0, sign;
1363 realtype tfuzz, tBret, tBn;
1364 booleantype gotCkpnt, reachedTBout, isActive;
1365
1366 /* Is the mem OK? */
1367 if (ida_mem == NULL) {
1368 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDASolveB", MSGAM_NULL_IDAMEM);
1369 return(IDA_MEM_NULL);
1370 }
1371 IDA_mem = (IDAMem) ida_mem;
1372
1373 /* Is ASA initialized ? */
1374 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
1375 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDASolveB", MSGAM_NO_ADJ);
1376 return(IDA_NO_ADJ);
1377 }
1378 IDAADJ_mem = IDA_mem->ida_adj_mem;
1379
1380 if ( IDAADJ_mem->ia_nbckpbs == 0 ) {
1381 IDAProcessError(IDA_mem, IDA_NO_BCK, "IDAA", "IDASolveB", MSGAM_NO_BCK);
1382 return(IDA_NO_BCK);
1383 }
1384 IDAB_mem = IDAADJ_mem->IDAB_mem;
1385
1386 /* Check whether IDASolveF has been called */
1387 if ( IDAADJ_mem->ia_firstIDAFcall ) {
1388 IDAProcessError(IDA_mem, IDA_NO_FWD, "IDAA", "IDASolveB", MSGAM_NO_FWD);
1389 return(IDA_NO_FWD);
1390 }
1391 sign = (IDAADJ_mem->ia_tfinal - IDAADJ_mem->ia_tinitial > ZERO) ? 1 : -1;
1392
1393 /* If this is the first call, loop over all backward problems and
1394 * - check that tB0 is valid
1395 * - check that tBout is ahead of tB0 in the backward direction
1396 * - check whether we need to interpolate forward sensitivities
1397 */
1398 if (IDAADJ_mem->ia_firstIDABcall) {
1399
1400 /* First IDABMem struct. */
1401 tmp_IDAB_mem = IDAB_mem;
1402
1403 while (tmp_IDAB_mem != NULL) {
1404
1405 tBn = tmp_IDAB_mem->IDA_mem->ida_tn;
1406
1407 if ( (sign*(tBn-IDAADJ_mem->ia_tinitial) < ZERO) || (sign*(IDAADJ_mem->ia_tfinal-tBn) < ZERO) ) {
1408 IDAProcessError(IDA_mem, IDA_BAD_TB0, "IDAA", "IDASolveB",
1409 MSGAM_BAD_TB0, tmp_IDAB_mem->ida_index);
1410 return(IDA_BAD_TB0);
1411 }
1412
1413 if (sign*(tBn-tBout) <= ZERO) {
1414 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDASolveB", MSGAM_BAD_TBOUT,
1415 tmp_IDAB_mem->ida_index);
1416 return(IDA_ILL_INPUT);
1417 }
1418
1419 if ( tmp_IDAB_mem->ida_res_withSensi ||
1420 tmp_IDAB_mem->ida_rhsQ_withSensi )
1421 IDAADJ_mem->ia_interpSensi = SUNTRUE;
1422
1423 /* Advance in list. */
1424 tmp_IDAB_mem = tmp_IDAB_mem->ida_next;
1425 }
1426
1427 if ( IDAADJ_mem->ia_interpSensi && !IDAADJ_mem->ia_storeSensi) {
1428 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDASolveB", MSGAM_BAD_SENSI);
1429 return(IDA_ILL_INPUT);
1430 }
1431
1432 IDAADJ_mem->ia_firstIDABcall = SUNFALSE;
1433 }
1434
1435 /* Check for valid itask */
1436 if ( (itaskB != IDA_NORMAL) && (itaskB != IDA_ONE_STEP) ) {
1437 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDASolveB", MSG_BAD_ITASK);
1438 return(IDA_ILL_INPUT);
1439 }
1440
1441 /* Check if tBout is legal */
1442 if ( (sign*(tBout-IDAADJ_mem->ia_tinitial) < ZERO) || (sign*(IDAADJ_mem->ia_tfinal-tBout) < ZERO) ) {
1443 tfuzz = HUNDRED * IDA_mem->ida_uround *
1444 (SUNRabs(IDAADJ_mem->ia_tinitial) + SUNRabs(IDAADJ_mem->ia_tfinal));
1445 if ( (sign*(tBout-IDAADJ_mem->ia_tinitial) < ZERO) && (SUNRabs(tBout-IDAADJ_mem->ia_tinitial) < tfuzz) ) {
1446 tBout = IDAADJ_mem->ia_tinitial;
1447 } else {
1448 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDASolveB", MSGAM_BAD_TBOUT);
1449 return(IDA_ILL_INPUT);
1450 }
1451 }
1452
1453 /* Loop through the check points and stop as soon as a backward
1454 * problem has its tn value behind the current check point's t0_
1455 * value (in the backward direction) */
1456
1457 ck_mem = IDAADJ_mem->ck_mem;
1458
1459 gotCkpnt = SUNFALSE;
1460
1461 for(;;) {
1462 tmp_IDAB_mem = IDAB_mem;
1463 while(tmp_IDAB_mem != NULL) {
1464 tBn = tmp_IDAB_mem->IDA_mem->ida_tn;
1465
1466 if ( sign*(tBn-ck_mem->ck_t0) > ZERO ) {
1467 gotCkpnt = SUNTRUE;
1468 break;
1469 }
1470
1471 if ( (itaskB == IDA_NORMAL) && (tBn == ck_mem->ck_t0) && (sign*(tBout-ck_mem->ck_t0) >= ZERO) ) {
1472 gotCkpnt = SUNTRUE;
1473 break;
1474 }
1475
1476 tmp_IDAB_mem = tmp_IDAB_mem->ida_next;
1477 }
1478
1479 if (gotCkpnt) break;
1480
1481 if (ck_mem->ck_next == NULL) break;
1482
1483 ck_mem = ck_mem->ck_next;
1484 }
1485
1486 /* Loop while propagating backward problems */
1487 for(;;) {
1488
1489 /* Store interpolation data if not available.
1490 This is the 2nd forward integration pass */
1491 if (ck_mem != IDAADJ_mem->ia_ckpntData) {
1492
1493 flag = IDAAdataStore(IDA_mem, ck_mem);
1494 if (flag != IDA_SUCCESS) break;
1495 }
1496
1497 /* Starting with the current check point from above, loop over check points
1498 while propagating backward problems */
1499
1500 tmp_IDAB_mem = IDAB_mem;
1501 while (tmp_IDAB_mem != NULL) {
1502
1503 /* Decide if current backward problem is "active" in this check point */
1504 isActive = SUNTRUE;
1505
1506 tBn = tmp_IDAB_mem->IDA_mem->ida_tn;
1507
1508 if ( (tBn == ck_mem->ck_t0) && (sign*(tBout-ck_mem->ck_t0) < ZERO ) ) isActive = SUNFALSE;
1509 if ( (tBn == ck_mem->ck_t0) && (itaskB == IDA_ONE_STEP) ) isActive = SUNFALSE;
1510 if ( sign*(tBn - ck_mem->ck_t0) < ZERO ) isActive = SUNFALSE;
1511
1512 if ( isActive ) {
1513 /* Store the address of current backward problem memory
1514 * in IDAADJ_mem to be used in the wrapper functions */
1515 IDAADJ_mem->ia_bckpbCrt = tmp_IDAB_mem;
1516
1517 /* Integrate current backward problem */
1518 IDASetStopTime(tmp_IDAB_mem->IDA_mem, ck_mem->ck_t0);
1519 flag = IDASolve(tmp_IDAB_mem->IDA_mem, tBout, &tBret,
1520 tmp_IDAB_mem->ida_yy, tmp_IDAB_mem->ida_yp,
1521 itaskB);
1522
1523 /* Set the time at which we will report solution and/or quadratures */
1524 tmp_IDAB_mem->ida_tout = tBret;
1525
1526 /* If an error occurred, exit while loop */
1527 if (flag < 0) break;
1528
1529 } else {
1530
1531 flag = IDA_SUCCESS;
1532 tmp_IDAB_mem->ida_tout = tBn;
1533 }
1534
1535 /* Move to next backward problem */
1536 tmp_IDAB_mem = tmp_IDAB_mem->ida_next;
1537 } /* End of while: iteration through backward problems. */
1538
1539 /* If an error occurred, return now */
1540 if (flag <0) {
1541 IDAProcessError(IDA_mem, flag, "IDAA", "IDASolveB",
1542 MSGAM_BACK_ERROR, tmp_IDAB_mem->ida_index);
1543 return(flag);
1544 }
1545
1546 /* If in IDA_ONE_STEP mode, return now (flag = IDA_SUCCESS) */
1547 if (itaskB == IDA_ONE_STEP) break;
1548
1549 /* If all backward problems have succesfully reached tBout, return now */
1550 reachedTBout = SUNTRUE;
1551
1552 tmp_IDAB_mem = IDAB_mem;
1553 while(tmp_IDAB_mem != NULL) {
1554 if ( sign*(tmp_IDAB_mem->ida_tout - tBout) > ZERO ) {
1555 reachedTBout = SUNFALSE;
1556 break;
1557 }
1558 tmp_IDAB_mem = tmp_IDAB_mem->ida_next;
1559 }
1560
1561 if ( reachedTBout ) break;
1562
1563 /* Move check point in linked list to next one */
1564 ck_mem = ck_mem->ck_next;
1565
1566 } /* End of loop. */
1567
1568 return(flag);
1569 }
1570
1571
1572 /*
1573 * IDAGetB
1574 *
1575 * IDAGetB returns the state variables at the same time (also returned
1576 * in tret) as that at which IDASolveBreturned the solution.
1577 */
1578
IDAGetB(void * ida_mem,int which,realtype * tret,N_Vector yy,N_Vector yp)1579 int IDAGetB(void* ida_mem, int which, realtype *tret,
1580 N_Vector yy, N_Vector yp)
1581 {
1582 IDAMem IDA_mem;
1583 IDAadjMem IDAADJ_mem;
1584 IDABMem IDAB_mem;
1585
1586 /* Is ida_mem valid? */
1587 if (ida_mem == NULL) {
1588 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAGetB", MSGAM_NULL_IDAMEM);
1589 return IDA_MEM_NULL;
1590 }
1591 IDA_mem = (IDAMem) ida_mem;
1592
1593 /* Is ASA initialized? */
1594 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
1595 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAGetB", MSGAM_NO_ADJ);
1596 return(IDA_NO_ADJ);
1597 }
1598 IDAADJ_mem = IDA_mem->ida_adj_mem;
1599
1600 /* Check the value of which */
1601 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
1602 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAGetB", MSGAM_BAD_WHICH);
1603 return(IDA_ILL_INPUT);
1604 }
1605
1606 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
1607 IDAB_mem = IDAADJ_mem->IDAB_mem;
1608 while (IDAB_mem != NULL) {
1609 if( which == IDAB_mem->ida_index ) break;
1610 /* advance */
1611 IDAB_mem = IDAB_mem->ida_next;
1612 }
1613
1614 N_VScale(ONE, IDAB_mem->ida_yy, yy);
1615 N_VScale(ONE, IDAB_mem->ida_yp, yp);
1616 *tret = IDAB_mem->ida_tout;
1617
1618 return(IDA_SUCCESS);
1619 }
1620
1621
1622
1623 /*
1624 * IDAGetQuadB
1625 *
1626 * IDAGetQuadB returns the quadrature variables at the same
1627 * time (also returned in tret) as that at which IDASolveB
1628 * returned the solution.
1629 */
1630
IDAGetQuadB(void * ida_mem,int which,realtype * tret,N_Vector qB)1631 int IDAGetQuadB(void *ida_mem, int which, realtype *tret, N_Vector qB)
1632 {
1633 IDAMem IDA_mem;
1634 IDAadjMem IDAADJ_mem;
1635 IDABMem IDAB_mem;
1636 void *ida_memB;
1637 int flag;
1638 long int nstB;
1639
1640 /* Is ida_mem valid? */
1641 if (ida_mem == NULL) {
1642 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAGetQuadB", MSGAM_NULL_IDAMEM);
1643 return IDA_MEM_NULL;
1644 }
1645 IDA_mem = (IDAMem) ida_mem;
1646
1647 /* Is ASA initialized? */
1648 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
1649 IDAProcessError(IDA_mem, IDA_NO_ADJ, "IDAA", "IDAGetQuadB", MSGAM_NO_ADJ);
1650 return(IDA_NO_ADJ);
1651 }
1652 IDAADJ_mem = IDA_mem->ida_adj_mem;
1653
1654 /* Check the value of which */
1655 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
1656 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDAA", "IDAGetQuadB", MSGAM_BAD_WHICH);
1657 return(IDA_ILL_INPUT);
1658 }
1659
1660 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
1661 IDAB_mem = IDAADJ_mem->IDAB_mem;
1662 while (IDAB_mem != NULL) {
1663 if( which == IDAB_mem->ida_index ) break;
1664 /* advance */
1665 IDAB_mem = IDAB_mem->ida_next;
1666 }
1667 ida_memB = (void *) IDAB_mem->IDA_mem;
1668
1669 /* If the integration for this backward problem has not started yet,
1670 * simply return the current value of qB (i.e. the final conditions) */
1671
1672 flag = IDAGetNumSteps(ida_memB, &nstB);
1673 if (IDA_SUCCESS != flag) return(flag);
1674
1675 if (nstB == 0) {
1676 N_VScale(ONE, IDAB_mem->IDA_mem->ida_phiQ[0], qB);
1677 *tret = IDAB_mem->ida_tout;
1678 } else {
1679 flag = IDAGetQuad(ida_memB, tret, qB);
1680 }
1681 return(flag);
1682 }
1683
1684 /*=================================================================*/
1685 /* Private Functions Implementation */
1686 /*=================================================================*/
1687
1688 /*
1689 * IDAAckpntInit
1690 *
1691 * This routine initializes the check point linked list with
1692 * information from the initial time.
1693 */
1694
IDAAckpntInit(IDAMem IDA_mem)1695 static CkpntMem IDAAckpntInit(IDAMem IDA_mem)
1696 {
1697 CkpntMem ck_mem;
1698
1699 /* Allocate space for ckdata */
1700 ck_mem = (CkpntMem) malloc(sizeof(struct CkpntMemRec));
1701 if (NULL==ck_mem) return(NULL);
1702
1703 ck_mem->ck_t0 = IDA_mem->ida_tn;
1704 ck_mem->ck_nst = 0;
1705 ck_mem->ck_kk = 1;
1706 ck_mem->ck_hh = ZERO;
1707
1708 /* Test if we need to carry quadratures */
1709 ck_mem->ck_quadr = IDA_mem->ida_quadr && IDA_mem->ida_errconQ;
1710
1711 /* Test if we need to carry sensitivities */
1712 ck_mem->ck_sensi = IDA_mem->ida_sensi;
1713 if(ck_mem->ck_sensi) ck_mem->ck_Ns = IDA_mem->ida_Ns;
1714
1715 /* Test if we need to carry quadrature sensitivities */
1716 ck_mem->ck_quadr_sensi = IDA_mem->ida_quadr_sensi && IDA_mem->ida_errconQS;
1717
1718 /* Alloc 3: current order, i.e. 1, + 2. */
1719 ck_mem->ck_phi_alloc = 3;
1720
1721 if (!IDAAckpntAllocVectors(IDA_mem, ck_mem)) {
1722 free(ck_mem); ck_mem = NULL;
1723 return(NULL);
1724 }
1725 /* Save phi* vectors from IDA_mem to ck_mem. */
1726 IDAAckpntCopyVectors(IDA_mem, ck_mem);
1727
1728 /* Next in list */
1729 ck_mem->ck_next = NULL;
1730
1731 return(ck_mem);
1732 }
1733
1734 /*
1735 * IDAAckpntNew
1736 *
1737 * This routine allocates space for a new check point and sets
1738 * its data from current values in IDA_mem.
1739 */
1740
IDAAckpntNew(IDAMem IDA_mem)1741 static CkpntMem IDAAckpntNew(IDAMem IDA_mem)
1742 {
1743 CkpntMem ck_mem;
1744 int j;
1745
1746 /* Allocate space for ckdata */
1747 ck_mem = (CkpntMem) malloc(sizeof(struct CkpntMemRec));
1748 if (ck_mem == NULL) return(NULL);
1749
1750 ck_mem->ck_nst = IDA_mem->ida_nst;
1751 ck_mem->ck_tretlast = IDA_mem->ida_tretlast;
1752 ck_mem->ck_kk = IDA_mem->ida_kk;
1753 ck_mem->ck_kused = IDA_mem->ida_kused;
1754 ck_mem->ck_knew = IDA_mem->ida_knew;
1755 ck_mem->ck_phase = IDA_mem->ida_phase;
1756 ck_mem->ck_ns = IDA_mem->ida_ns;
1757 ck_mem->ck_hh = IDA_mem->ida_hh;
1758 ck_mem->ck_hused = IDA_mem->ida_hused;
1759 ck_mem->ck_rr = IDA_mem->ida_rr;
1760 ck_mem->ck_cj = IDA_mem->ida_cj;
1761 ck_mem->ck_cjlast = IDA_mem->ida_cjlast;
1762 ck_mem->ck_cjold = IDA_mem->ida_cjold;
1763 ck_mem->ck_cjratio = IDA_mem->ida_cjratio;
1764 ck_mem->ck_ss = IDA_mem->ida_ss;
1765 ck_mem->ck_ssS = IDA_mem->ida_ssS;
1766 ck_mem->ck_t0 = IDA_mem->ida_tn;
1767
1768 for (j=0; j<MXORDP1; j++) {
1769 ck_mem->ck_psi[j] = IDA_mem->ida_psi[j];
1770 ck_mem->ck_alpha[j] = IDA_mem->ida_alpha[j];
1771 ck_mem->ck_beta[j] = IDA_mem->ida_beta[j];
1772 ck_mem->ck_sigma[j] = IDA_mem->ida_sigma[j];
1773 ck_mem->ck_gamma[j] = IDA_mem->ida_gamma[j];
1774 }
1775
1776 /* Test if we need to carry quadratures */
1777 ck_mem->ck_quadr = IDA_mem->ida_quadr && IDA_mem->ida_errconQ;
1778
1779 /* Test if we need to carry sensitivities */
1780 ck_mem->ck_sensi = IDA_mem->ida_sensi;
1781 if(ck_mem->ck_sensi) ck_mem->ck_Ns = IDA_mem->ida_Ns;
1782
1783 /* Test if we need to carry quadrature sensitivities */
1784 ck_mem->ck_quadr_sensi = IDA_mem->ida_quadr_sensi && IDA_mem->ida_errconQS;
1785
1786 ck_mem->ck_phi_alloc = (IDA_mem->ida_kk+2 < MXORDP1) ?
1787 IDA_mem->ida_kk+2 : MXORDP1;
1788
1789 if (!IDAAckpntAllocVectors(IDA_mem, ck_mem)) {
1790 free(ck_mem); ck_mem = NULL;
1791 return(NULL);
1792 }
1793
1794 /* Save phi* vectors from IDA_mem to ck_mem. */
1795 IDAAckpntCopyVectors(IDA_mem, ck_mem);
1796
1797 return(ck_mem);
1798 }
1799
1800 /* IDAAckpntDelete
1801 *
1802 * This routine deletes the first check point in list.
1803 */
1804
IDAAckpntDelete(CkpntMem * ck_memPtr)1805 static void IDAAckpntDelete(CkpntMem *ck_memPtr)
1806 {
1807 CkpntMem tmp;
1808 int j;
1809
1810 if (*ck_memPtr != NULL) {
1811 /* store head of list */
1812 tmp = *ck_memPtr;
1813 /* move head of list */
1814 *ck_memPtr = (*ck_memPtr)->ck_next;
1815
1816 /* free N_Vectors in tmp */
1817 for (j=0; j<tmp->ck_phi_alloc; j++)
1818 N_VDestroy(tmp->ck_phi[j]);
1819
1820 /* free N_Vectors for quadratures in tmp */
1821 if (tmp->ck_quadr) {
1822 for (j=0; j<tmp->ck_phi_alloc; j++)
1823 N_VDestroy(tmp->ck_phiQ[j]);
1824 }
1825
1826 /* Free sensitivity related data. */
1827 if (tmp->ck_sensi) {
1828 for (j=0; j<tmp->ck_phi_alloc; j++)
1829 N_VDestroyVectorArray(tmp->ck_phiS[j], tmp->ck_Ns);
1830 }
1831
1832 if (tmp->ck_quadr_sensi) {
1833 for (j=0; j<tmp->ck_phi_alloc; j++)
1834 N_VDestroyVectorArray(tmp->ck_phiQS[j], tmp->ck_Ns);
1835 }
1836
1837 free(tmp); tmp=NULL;
1838 }
1839 }
1840
1841 /*
1842 * IDAAckpntAllocVectors
1843 *
1844 * Allocate checkpoint's phi, phiQ, phiS, phiQS vectors needed to save
1845 * current state of IDAMem.
1846 *
1847 */
IDAAckpntAllocVectors(IDAMem IDA_mem,CkpntMem ck_mem)1848 static booleantype IDAAckpntAllocVectors(IDAMem IDA_mem, CkpntMem ck_mem)
1849 {
1850 int j, jj;
1851
1852 for (j=0; j<ck_mem->ck_phi_alloc; j++) {
1853 ck_mem->ck_phi[j] = N_VClone(IDA_mem->ida_tempv1);
1854 if(ck_mem->ck_phi[j] == NULL) {
1855 for(jj=0; jj<j; jj++) N_VDestroy(ck_mem->ck_phi[jj]);
1856 return(SUNFALSE);
1857 }
1858 }
1859
1860 /* Do we need to carry quadratures? */
1861 if(ck_mem->ck_quadr) {
1862 for (j=0; j<ck_mem->ck_phi_alloc; j++) {
1863 ck_mem->ck_phiQ[j] = N_VClone(IDA_mem->ida_eeQ);
1864 if(ck_mem->ck_phiQ[j] == NULL) {
1865 for (jj=0; jj<j; jj++) N_VDestroy(ck_mem->ck_phiQ[jj]);
1866
1867 for(jj=0; jj<ck_mem->ck_phi_alloc; jj++)
1868 N_VDestroy(ck_mem->ck_phi[jj]);
1869
1870 return(SUNFALSE);
1871 }
1872 }
1873 }
1874
1875 /* Do we need to carry sensitivities? */
1876 if(ck_mem->ck_sensi) {
1877
1878 for (j=0; j<ck_mem->ck_phi_alloc; j++) {
1879 ck_mem->ck_phiS[j] = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_tempv1);
1880 if (ck_mem->ck_phiS[j] == NULL) {
1881 for (jj=0; jj<j; jj++)
1882 N_VDestroyVectorArray(ck_mem->ck_phiS[jj], IDA_mem->ida_Ns);
1883
1884 if (ck_mem->ck_quadr)
1885 for (jj=0; jj<ck_mem->ck_phi_alloc; jj++)
1886 N_VDestroy(ck_mem->ck_phiQ[jj]);
1887
1888 for (jj=0; jj<ck_mem->ck_phi_alloc; jj++)
1889 N_VDestroy(ck_mem->ck_phi[jj]);
1890
1891 return(SUNFALSE);
1892 }
1893 }
1894 }
1895
1896 /* Do we need to carry quadrature sensitivities? */
1897 if (ck_mem->ck_quadr_sensi) {
1898
1899 for (j=0; j<ck_mem->ck_phi_alloc; j++) {
1900 ck_mem->ck_phiQS[j] = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_eeQ);
1901 if (ck_mem->ck_phiQS[j] == NULL) {
1902
1903 for (jj=0; jj<j; jj++)
1904 N_VDestroyVectorArray(ck_mem->ck_phiQS[jj], IDA_mem->ida_Ns);
1905
1906 for (jj=0; jj<ck_mem->ck_phi_alloc; jj++)
1907 N_VDestroyVectorArray(ck_mem->ck_phiS[jj], IDA_mem->ida_Ns);
1908
1909 if (ck_mem->ck_quadr)
1910 for (jj=0; jj<ck_mem->ck_phi_alloc; jj++)
1911 N_VDestroy(ck_mem->ck_phiQ[jj]);
1912
1913 for (jj=0; jj<ck_mem->ck_phi_alloc; jj++)
1914 N_VDestroy(ck_mem->ck_phi[jj]);
1915
1916 return(SUNFALSE);
1917 }
1918 }
1919 }
1920 return(SUNTRUE);
1921 }
1922
1923 /*
1924 * IDAAckpntCopyVectors
1925 *
1926 * Copy phi* vectors from IDAMem in the corresponding vectors from checkpoint
1927 *
1928 */
IDAAckpntCopyVectors(IDAMem IDA_mem,CkpntMem ck_mem)1929 static void IDAAckpntCopyVectors(IDAMem IDA_mem, CkpntMem ck_mem)
1930 {
1931 int j, is;
1932
1933 /* Save phi* arrays from IDA_mem */
1934
1935 for (j=0; j<ck_mem->ck_phi_alloc; j++)
1936 IDA_mem->ida_cvals[j] = ONE;
1937
1938 (void) N_VScaleVectorArray(ck_mem->ck_phi_alloc, IDA_mem->ida_cvals,
1939 IDA_mem->ida_phi, ck_mem->ck_phi);
1940
1941 if (ck_mem->ck_quadr)
1942 (void) N_VScaleVectorArray(ck_mem->ck_phi_alloc, IDA_mem->ida_cvals,
1943 IDA_mem->ida_phiQ, ck_mem->ck_phiQ);
1944
1945 if (ck_mem->ck_sensi || ck_mem->ck_quadr_sensi) {
1946 for (j=0; j<ck_mem->ck_phi_alloc; j++) {
1947 for (is=0; is<IDA_mem->ida_Ns; is++) {
1948 IDA_mem->ida_cvals[j*IDA_mem->ida_Ns + is] = ONE;
1949 }
1950 }
1951 }
1952
1953 if (ck_mem->ck_sensi) {
1954 for (j=0; j<ck_mem->ck_phi_alloc; j++) {
1955 for (is=0; is<IDA_mem->ida_Ns; is++) {
1956 IDA_mem->ida_Xvecs[j*IDA_mem->ida_Ns + is] = IDA_mem->ida_phiS[j][is];
1957 IDA_mem->ida_Zvecs[j*IDA_mem->ida_Ns + is] = ck_mem->ck_phiS[j][is];
1958 }
1959 }
1960
1961 (void) N_VScaleVectorArray(ck_mem->ck_phi_alloc * IDA_mem->ida_Ns,
1962 IDA_mem->ida_cvals,
1963 IDA_mem->ida_Xvecs, IDA_mem->ida_Zvecs);
1964 }
1965
1966 if(ck_mem->ck_quadr_sensi) {
1967 for (j=0; j<ck_mem->ck_phi_alloc; j++) {
1968 for (is=0; is<IDA_mem->ida_Ns; is++) {
1969 IDA_mem->ida_Xvecs[j*IDA_mem->ida_Ns + is] = IDA_mem->ida_phiQS[j][is];
1970 IDA_mem->ida_Zvecs[j*IDA_mem->ida_Ns + is] = ck_mem->ck_phiQS[j][is];
1971 }
1972 }
1973
1974 (void) N_VScaleVectorArray(ck_mem->ck_phi_alloc * IDA_mem->ida_Ns,
1975 IDA_mem->ida_cvals,
1976 IDA_mem->ida_Xvecs, IDA_mem->ida_Zvecs);
1977 }
1978
1979 }
1980
1981 /*
1982 * IDAAdataMalloc
1983 *
1984 * This routine allocates memory for storing information at all
1985 * intermediate points between two consecutive check points.
1986 * This data is then used to interpolate the forward solution
1987 * at any other time.
1988 */
1989
IDAAdataMalloc(IDAMem IDA_mem)1990 static booleantype IDAAdataMalloc(IDAMem IDA_mem)
1991 {
1992 IDAadjMem IDAADJ_mem;
1993 DtpntMem *dt_mem;
1994 long int i, j;
1995
1996 IDAADJ_mem = IDA_mem->ida_adj_mem;
1997 IDAADJ_mem->dt_mem = NULL;
1998
1999 dt_mem = (DtpntMem *)malloc((IDAADJ_mem->ia_nsteps+1)*sizeof(struct DtpntMemRec *));
2000 if (dt_mem==NULL) return(SUNFALSE);
2001
2002 for (i=0; i<=IDAADJ_mem->ia_nsteps; i++) {
2003
2004 dt_mem[i] = (DtpntMem)malloc(sizeof(struct DtpntMemRec));
2005
2006 /* On failure, free any allocated memory and return NULL. */
2007 if (dt_mem[i] == NULL) {
2008
2009 for(j=0; j<i; j++)
2010 free(dt_mem[j]);
2011
2012 free(dt_mem);
2013 return(SUNFALSE);
2014 }
2015 dt_mem[i]->content = NULL;
2016 }
2017 /* Attach the allocated dt_mem to IDAADJ_mem. */
2018 IDAADJ_mem->dt_mem = dt_mem;
2019 return(SUNTRUE);
2020 }
2021
2022 /*
2023 * IDAAdataFree
2024 *
2025 * This routine frees the memory allocated for data storage.
2026 */
2027
IDAAdataFree(IDAMem IDA_mem)2028 static void IDAAdataFree(IDAMem IDA_mem)
2029 {
2030 IDAadjMem IDAADJ_mem;
2031 long int i;
2032
2033 IDAADJ_mem = IDA_mem->ida_adj_mem;
2034
2035 if (IDAADJ_mem == NULL) return;
2036
2037 /* Destroy data points by calling the interpolation's 'free' routine. */
2038 IDAADJ_mem->ia_free(IDA_mem);
2039
2040 for (i=0; i<=IDAADJ_mem->ia_nsteps; i++) {
2041 free(IDAADJ_mem->dt_mem[i]);
2042 IDAADJ_mem->dt_mem[i] = NULL;
2043 }
2044
2045 free(IDAADJ_mem->dt_mem);
2046 IDAADJ_mem->dt_mem = NULL;
2047 }
2048
2049
2050 /*
2051 * IDAAdataStore
2052 *
2053 * This routine integrates the forward model starting at the check
2054 * point ck_mem and stores y and yprime at all intermediate
2055 * steps.
2056 *
2057 * Return values:
2058 * - the flag that IDASolve may return on error
2059 * - IDA_REIFWD_FAIL if no check point is available for this hot start
2060 * - IDA_SUCCESS
2061 */
2062
IDAAdataStore(IDAMem IDA_mem,CkpntMem ck_mem)2063 static int IDAAdataStore(IDAMem IDA_mem, CkpntMem ck_mem)
2064 {
2065 IDAadjMem IDAADJ_mem;
2066 DtpntMem *dt_mem;
2067 realtype t;
2068 long int i;
2069 int flag, sign;
2070
2071 IDAADJ_mem = IDA_mem->ida_adj_mem;
2072 dt_mem = IDAADJ_mem->dt_mem;
2073
2074 /* Initialize IDA_mem with data from ck_mem. */
2075 flag = IDAAckpntGet(IDA_mem, ck_mem);
2076 if (flag != IDA_SUCCESS)
2077 return(IDA_REIFWD_FAIL);
2078
2079 /* Set first structure in dt_mem[0] */
2080 dt_mem[0]->t = ck_mem->ck_t0;
2081 IDAADJ_mem->ia_storePnt(IDA_mem, dt_mem[0]);
2082
2083 /* Decide whether TSTOP must be activated */
2084 if (IDAADJ_mem->ia_tstopIDAFcall) {
2085 IDASetStopTime(IDA_mem, IDAADJ_mem->ia_tstopIDAF);
2086 }
2087
2088 sign = (IDAADJ_mem->ia_tfinal - IDAADJ_mem->ia_tinitial > ZERO) ? 1 : -1;
2089
2090 /* Run IDASolve in IDA_ONE_STEP mode to set following structures in dt_mem[i]. */
2091 i = 1;
2092 do {
2093
2094 flag = IDASolve(IDA_mem, ck_mem->ck_t1, &t, IDAADJ_mem->ia_yyTmp,
2095 IDAADJ_mem->ia_ypTmp, IDA_ONE_STEP);
2096 if (flag < 0) return(IDA_FWD_FAIL);
2097
2098 dt_mem[i]->t = t;
2099 IDAADJ_mem->ia_storePnt(IDA_mem, dt_mem[i]);
2100
2101 i++;
2102 } while ( sign*(ck_mem->ck_t1 - t) > ZERO );
2103
2104 /* New data is now available. */
2105 IDAADJ_mem->ia_ckpntData = ck_mem;
2106 IDAADJ_mem->ia_newData = SUNTRUE;
2107 IDAADJ_mem->ia_np = i;
2108
2109 return(IDA_SUCCESS);
2110 }
2111
2112 /*
2113 * CVAckpntGet
2114 *
2115 * This routine prepares IDAS for a hot restart from
2116 * the check point ck_mem
2117 */
2118
IDAAckpntGet(IDAMem IDA_mem,CkpntMem ck_mem)2119 static int IDAAckpntGet(IDAMem IDA_mem, CkpntMem ck_mem)
2120 {
2121 int flag, j, is;
2122
2123 if (ck_mem->ck_next == NULL) {
2124
2125 /* In this case, we just call the reinitialization routine,
2126 * but make sure we use the same initial stepsize as on
2127 * the first run. */
2128
2129 IDASetInitStep(IDA_mem, IDA_mem->ida_h0u);
2130
2131 flag = IDAReInit(IDA_mem, ck_mem->ck_t0, ck_mem->ck_phi[0], ck_mem->ck_phi[1]);
2132 if (flag != IDA_SUCCESS) return(flag);
2133
2134 if (ck_mem->ck_quadr) {
2135 flag = IDAQuadReInit(IDA_mem, ck_mem->ck_phiQ[0]);
2136 if (flag != IDA_SUCCESS) return(flag);
2137 }
2138
2139 if (ck_mem->ck_sensi) {
2140 flag = IDASensReInit(IDA_mem, IDA_mem->ida_ism, ck_mem->ck_phiS[0], ck_mem->ck_phiS[1]);
2141 if (flag != IDA_SUCCESS) return(flag);
2142 }
2143
2144 if (ck_mem->ck_quadr_sensi) {
2145 flag = IDAQuadSensReInit(IDA_mem, ck_mem->ck_phiQS[0]);
2146 if (flag != IDA_SUCCESS) return(flag);
2147 }
2148
2149 } else {
2150
2151 /* Copy parameters from check point data structure */
2152 IDA_mem->ida_nst = ck_mem->ck_nst;
2153 IDA_mem->ida_tretlast = ck_mem->ck_tretlast;
2154 IDA_mem->ida_kk = ck_mem->ck_kk;
2155 IDA_mem->ida_kused = ck_mem->ck_kused;
2156 IDA_mem->ida_knew = ck_mem->ck_knew;
2157 IDA_mem->ida_phase = ck_mem->ck_phase;
2158 IDA_mem->ida_ns = ck_mem->ck_ns;
2159 IDA_mem->ida_hh = ck_mem->ck_hh;
2160 IDA_mem->ida_hused = ck_mem->ck_hused;
2161 IDA_mem->ida_rr = ck_mem->ck_rr;
2162 IDA_mem->ida_cj = ck_mem->ck_cj;
2163 IDA_mem->ida_cjlast = ck_mem->ck_cjlast;
2164 IDA_mem->ida_cjold = ck_mem->ck_cjold;
2165 IDA_mem->ida_cjratio = ck_mem->ck_cjratio;
2166 IDA_mem->ida_tn = ck_mem->ck_t0;
2167 IDA_mem->ida_ss = ck_mem->ck_ss;
2168 IDA_mem->ida_ssS = ck_mem->ck_ssS;
2169
2170
2171 /* Copy the arrays from check point data structure */
2172 for (j=0; j<ck_mem->ck_phi_alloc; j++)
2173 N_VScale(ONE, ck_mem->ck_phi[j], IDA_mem->ida_phi[j]);
2174
2175 if(ck_mem->ck_quadr) {
2176 for (j=0; j<ck_mem->ck_phi_alloc; j++)
2177 N_VScale(ONE, ck_mem->ck_phiQ[j], IDA_mem->ida_phiQ[j]);
2178 }
2179
2180 if (ck_mem->ck_sensi) {
2181 for (is=0; is<IDA_mem->ida_Ns; is++) {
2182 for (j=0; j<ck_mem->ck_phi_alloc; j++)
2183 N_VScale(ONE, ck_mem->ck_phiS[j][is], IDA_mem->ida_phiS[j][is]);
2184 }
2185 }
2186
2187 if (ck_mem->ck_quadr_sensi) {
2188 for (is=0; is<IDA_mem->ida_Ns; is++) {
2189 for (j=0; j<ck_mem->ck_phi_alloc; j++)
2190 N_VScale(ONE, ck_mem->ck_phiQS[j][is], IDA_mem->ida_phiQS[j][is]);
2191 }
2192 }
2193
2194 for (j=0; j<MXORDP1; j++) {
2195 IDA_mem->ida_psi[j] = ck_mem->ck_psi[j];
2196 IDA_mem->ida_alpha[j] = ck_mem->ck_alpha[j];
2197 IDA_mem->ida_beta[j] = ck_mem->ck_beta[j];
2198 IDA_mem->ida_sigma[j] = ck_mem->ck_sigma[j];
2199 IDA_mem->ida_gamma[j] = ck_mem->ck_gamma[j];
2200 }
2201
2202 /* Force a call to setup */
2203 IDA_mem->ida_forceSetup = SUNTRUE;
2204 }
2205
2206 return(IDA_SUCCESS);
2207 }
2208
2209
2210 /*
2211 * -----------------------------------------------------------------
2212 * Functions specific to cubic Hermite interpolation
2213 * -----------------------------------------------------------------
2214 */
2215
2216 /*
2217 * IDAAhermiteMalloc
2218 *
2219 * This routine allocates memory for storing information at all
2220 * intermediate points between two consecutive check points.
2221 * This data is then used to interpolate the forward solution
2222 * at any other time.
2223 */
2224
IDAAhermiteMalloc(IDAMem IDA_mem)2225 static booleantype IDAAhermiteMalloc(IDAMem IDA_mem)
2226 {
2227 IDAadjMem IDAADJ_mem;
2228 DtpntMem *dt_mem;
2229 HermiteDataMem content;
2230 long int i, ii=0;
2231 booleantype allocOK;
2232
2233 allocOK = SUNTRUE;
2234
2235 IDAADJ_mem = IDA_mem->ida_adj_mem;
2236
2237 /* Allocate space for the vectors yyTmp and ypTmp. */
2238 IDAADJ_mem->ia_yyTmp = N_VClone(IDA_mem->ida_tempv1);
2239 if (IDAADJ_mem->ia_yyTmp == NULL) {
2240 return(SUNFALSE);
2241 }
2242 IDAADJ_mem->ia_ypTmp = N_VClone(IDA_mem->ida_tempv1);
2243 if (IDAADJ_mem->ia_ypTmp == NULL) {
2244 return(SUNFALSE);
2245 }
2246
2247 /* Allocate space for sensitivities temporary vectors. */
2248 if (IDAADJ_mem->ia_storeSensi) {
2249
2250 IDAADJ_mem->ia_yySTmp = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_tempv1);
2251 if (IDAADJ_mem->ia_yySTmp == NULL) {
2252 N_VDestroy(IDAADJ_mem->ia_yyTmp);
2253 N_VDestroy(IDAADJ_mem->ia_ypTmp);
2254 return(SUNFALSE);
2255 }
2256
2257 IDAADJ_mem->ia_ypSTmp = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_tempv1);
2258 if (IDAADJ_mem->ia_ypSTmp == NULL) {
2259 N_VDestroy(IDAADJ_mem->ia_yyTmp);
2260 N_VDestroy(IDAADJ_mem->ia_ypTmp);
2261 N_VDestroyVectorArray(IDAADJ_mem->ia_yySTmp, IDA_mem->ida_Ns);
2262 return(SUNFALSE);
2263
2264 }
2265 }
2266
2267 /* Allocate space for the content field of the dt structures */
2268
2269 dt_mem = IDAADJ_mem->dt_mem;
2270
2271 for (i=0; i<=IDAADJ_mem->ia_nsteps; i++) {
2272
2273 content = NULL;
2274 content = (HermiteDataMem) malloc(sizeof(struct HermiteDataMemRec));
2275 if (content == NULL) {
2276 ii = i;
2277 allocOK = SUNFALSE;
2278 break;
2279 }
2280
2281 content->y = N_VClone(IDA_mem->ida_tempv1);
2282 if (content->y == NULL) {
2283 free(content); content = NULL;
2284 ii = i;
2285 allocOK = SUNFALSE;
2286 break;
2287 }
2288
2289 content->yd = N_VClone(IDA_mem->ida_tempv1);
2290 if (content->yd == NULL) {
2291 N_VDestroy(content->y);
2292 free(content); content = NULL;
2293 ii = i;
2294 allocOK = SUNFALSE;
2295 break;
2296 }
2297
2298 if (IDAADJ_mem->ia_storeSensi) {
2299
2300 content->yS = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_tempv1);
2301 if (content->yS == NULL) {
2302 N_VDestroy(content->y);
2303 N_VDestroy(content->yd);
2304 free(content); content = NULL;
2305 ii = i;
2306 allocOK = SUNFALSE;
2307 break;
2308 }
2309
2310 content->ySd = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_tempv1);
2311 if (content->ySd == NULL) {
2312 N_VDestroy(content->y);
2313 N_VDestroy(content->yd);
2314 N_VDestroyVectorArray(content->yS, IDA_mem->ida_Ns);
2315 free(content); content = NULL;
2316 ii = i;
2317 allocOK = SUNFALSE;
2318 break;
2319 }
2320 }
2321
2322 dt_mem[i]->content = content;
2323
2324 }
2325
2326 /* If an error occurred, deallocate and return */
2327
2328 if (!allocOK) {
2329
2330 N_VDestroy(IDAADJ_mem->ia_yyTmp);
2331 N_VDestroy(IDAADJ_mem->ia_ypTmp);
2332
2333 if (IDAADJ_mem->ia_storeSensi) {
2334 N_VDestroyVectorArray(IDAADJ_mem->ia_yySTmp, IDA_mem->ida_Ns);
2335 N_VDestroyVectorArray(IDAADJ_mem->ia_ypSTmp, IDA_mem->ida_Ns);
2336 }
2337
2338 for (i=0; i<ii; i++) {
2339 content = (HermiteDataMem) (dt_mem[i]->content);
2340 N_VDestroy(content->y);
2341 N_VDestroy(content->yd);
2342
2343 if (IDAADJ_mem->ia_storeSensi) {
2344 N_VDestroyVectorArray(content->yS, IDA_mem->ida_Ns);
2345 N_VDestroyVectorArray(content->ySd, IDA_mem->ida_Ns);
2346 }
2347
2348 free(dt_mem[i]->content); dt_mem[i]->content = NULL;
2349 }
2350
2351 }
2352
2353 return(allocOK);
2354 }
2355
2356 /*
2357 * IDAAhermiteFree
2358 *
2359 * This routine frees the memory allocated for data storage.
2360 */
2361
IDAAhermiteFree(IDAMem IDA_mem)2362 static void IDAAhermiteFree(IDAMem IDA_mem)
2363 {
2364 IDAadjMem IDAADJ_mem;
2365 DtpntMem *dt_mem;
2366 HermiteDataMem content;
2367 long int i;
2368
2369 IDAADJ_mem = IDA_mem->ida_adj_mem;
2370
2371 N_VDestroy(IDAADJ_mem->ia_yyTmp);
2372 N_VDestroy(IDAADJ_mem->ia_ypTmp);
2373
2374 if (IDAADJ_mem->ia_storeSensi) {
2375 N_VDestroyVectorArray(IDAADJ_mem->ia_yySTmp, IDA_mem->ida_Ns);
2376 N_VDestroyVectorArray(IDAADJ_mem->ia_ypSTmp, IDA_mem->ida_Ns);
2377 }
2378
2379 dt_mem = IDAADJ_mem->dt_mem;
2380
2381 for (i=0; i<=IDAADJ_mem->ia_nsteps; i++) {
2382
2383 content = (HermiteDataMem) (dt_mem[i]->content);
2384 /* content might be NULL, if IDAAdjInit was called but IDASolveF was not. */
2385 if(content) {
2386
2387 N_VDestroy(content->y);
2388 N_VDestroy(content->yd);
2389
2390 if (IDAADJ_mem->ia_storeSensi) {
2391 N_VDestroyVectorArray(content->yS, IDA_mem->ida_Ns);
2392 N_VDestroyVectorArray(content->ySd, IDA_mem->ida_Ns);
2393 }
2394 free(dt_mem[i]->content);
2395 dt_mem[i]->content = NULL;
2396 }
2397 }
2398 }
2399
2400 /*
2401 * IDAAhermiteStorePnt
2402 *
2403 * This routine stores a new point (y,yd) in the structure d for use
2404 * in the cubic Hermite interpolation.
2405 * Note that the time is already stored.
2406 */
2407
IDAAhermiteStorePnt(IDAMem IDA_mem,DtpntMem d)2408 static int IDAAhermiteStorePnt(IDAMem IDA_mem, DtpntMem d)
2409 {
2410 IDAadjMem IDAADJ_mem;
2411 HermiteDataMem content;
2412 int is, retval;
2413
2414 IDAADJ_mem = IDA_mem->ida_adj_mem;
2415
2416 content = (HermiteDataMem) d->content;
2417
2418 /* Load solution(s) */
2419 N_VScale(ONE, IDA_mem->ida_phi[0], content->y);
2420
2421 if (IDAADJ_mem->ia_storeSensi) {
2422 for (is=0; is<IDA_mem->ida_Ns; is++)
2423 IDA_mem->ida_cvals[is] = ONE;
2424
2425 retval = N_VScaleVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_cvals,
2426 IDA_mem->ida_phiS[0], content->yS);
2427 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2428 }
2429
2430 /* Load derivative(s). */
2431 IDAAGettnSolutionYp(IDA_mem, content->yd);
2432
2433 if (IDAADJ_mem->ia_storeSensi) {
2434 IDAAGettnSolutionYpS(IDA_mem, content->ySd);
2435 }
2436
2437 return(0);
2438 }
2439
2440
2441 /*
2442 * IDAAhermiteGetY
2443 *
2444 * This routine uses cubic piece-wise Hermite interpolation for
2445 * the forward solution vector.
2446 * It is typically called by the wrapper routines before calling
2447 * user provided routines (fB, djacB, bjacB, jtimesB, psolB) but
2448 * can be directly called by the user through IDAGetAdjY
2449 */
2450
IDAAhermiteGetY(IDAMem IDA_mem,realtype t,N_Vector yy,N_Vector yp,N_Vector * yyS,N_Vector * ypS)2451 static int IDAAhermiteGetY(IDAMem IDA_mem, realtype t,
2452 N_Vector yy, N_Vector yp,
2453 N_Vector *yyS, N_Vector *ypS)
2454 {
2455 IDAadjMem IDAADJ_mem;
2456 DtpntMem *dt_mem;
2457 HermiteDataMem content0, content1;
2458
2459 realtype t0, t1, delta;
2460 realtype factor1, factor2, factor3;
2461
2462 N_Vector y0, yd0, y1, yd1;
2463 N_Vector *yS0=NULL, *ySd0=NULL, *yS1, *ySd1;
2464
2465 int flag, is, NS;
2466 long int indx;
2467 booleantype newpoint;
2468
2469 /* local variables for fused vector oerations */
2470 int retval;
2471 realtype cvals[4];
2472 N_Vector Xvecs[4];
2473 N_Vector* XXvecs[4];
2474
2475 IDAADJ_mem = IDA_mem->ida_adj_mem;
2476 dt_mem = IDAADJ_mem->dt_mem;
2477
2478 /* Local value of Ns */
2479 NS = (IDAADJ_mem->ia_interpSensi && (yyS != NULL)) ? IDA_mem->ida_Ns : 0;
2480
2481 /* Get the index in dt_mem */
2482 flag = IDAAfindIndex(IDA_mem, t, &indx, &newpoint);
2483 if (flag != IDA_SUCCESS) return(flag);
2484
2485 /* If we are beyond the left limit but close enough,
2486 then return y at the left limit. */
2487
2488 if (indx == 0) {
2489 content0 = (HermiteDataMem) (dt_mem[0]->content);
2490 N_VScale(ONE, content0->y, yy);
2491 N_VScale(ONE, content0->yd, yp);
2492
2493 if (NS > 0) {
2494 for (is=0; is<NS; is++)
2495 IDA_mem->ida_cvals[is] = ONE;
2496
2497 retval = N_VScaleVectorArray(NS, IDA_mem->ida_cvals, content0->yS, yyS);
2498 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2499
2500 retval = N_VScaleVectorArray(NS, IDA_mem->ida_cvals, content0->ySd, ypS);
2501 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2502 }
2503
2504 return(IDA_SUCCESS);
2505 }
2506
2507 /* Extract stuff from the appropriate data points */
2508 t0 = dt_mem[indx-1]->t;
2509 t1 = dt_mem[indx]->t;
2510 delta = t1 - t0;
2511
2512 content0 = (HermiteDataMem) (dt_mem[indx-1]->content);
2513 y0 = content0->y;
2514 yd0 = content0->yd;
2515 if (IDAADJ_mem->ia_interpSensi) {
2516 yS0 = content0->yS;
2517 ySd0 = content0->ySd;
2518 }
2519
2520 if (newpoint) {
2521
2522 /* Recompute Y0 and Y1 */
2523 content1 = (HermiteDataMem) (dt_mem[indx]->content);
2524
2525 y1 = content1->y;
2526 yd1 = content1->yd;
2527
2528 /* Y1 = delta (yd1 + yd0) - 2 (y1 - y0) */
2529 cvals[0] = -TWO; Xvecs[0] = y1;
2530 cvals[1] = TWO; Xvecs[1] = y0;
2531 cvals[2] = delta; Xvecs[2] = yd1;
2532 cvals[3] = delta; Xvecs[3] = yd0;
2533
2534 retval = N_VLinearCombination(4, cvals, Xvecs, IDAADJ_mem->ia_Y[1]);
2535 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2536
2537 /* Y0 = y1 - y0 - delta * yd0 */
2538 cvals[0] = ONE; Xvecs[0] = y1;
2539 cvals[1] = -ONE; Xvecs[1] = y0;
2540 cvals[2] = -delta; Xvecs[2] = yd0;
2541
2542 retval = N_VLinearCombination(3, cvals, Xvecs, IDAADJ_mem->ia_Y[0]);
2543 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2544
2545 /* Recompute YS0 and YS1, if needed */
2546
2547 if (NS > 0) {
2548
2549 yS1 = content1->yS;
2550 ySd1 = content1->ySd;
2551
2552 /* YS1 = delta (ySd1 + ySd0) - 2 (yS1 - yS0) */
2553 cvals[0] = -TWO; XXvecs[0] = yS1;
2554 cvals[1] = TWO; XXvecs[1] = yS0;
2555 cvals[2] = delta; XXvecs[2] = ySd1;
2556 cvals[3] = delta; XXvecs[3] = ySd0;
2557
2558 retval = N_VLinearCombinationVectorArray(NS, 4, cvals, XXvecs, IDAADJ_mem->ia_YS[1]);
2559 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2560
2561 /* YS0 = yS1 - yS0 - delta * ySd0 */
2562 cvals[0] = ONE; XXvecs[0] = yS1;
2563 cvals[1] = -ONE; XXvecs[1] = yS0;
2564 cvals[2] = -delta; XXvecs[2] = ySd0;
2565
2566 retval = N_VLinearCombinationVectorArray(NS, 3, cvals, XXvecs, IDAADJ_mem->ia_YS[0]);
2567 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2568
2569 }
2570
2571 }
2572
2573 /* Perform the actual interpolation. */
2574
2575 /* For y. */
2576 factor1 = t - t0;
2577
2578 factor2 = factor1/delta;
2579 factor2 = factor2*factor2;
2580
2581 factor3 = factor2*(t-t1)/delta;
2582
2583 cvals[0] = ONE;
2584 cvals[1] = factor1;
2585 cvals[2] = factor2;
2586 cvals[3] = factor3;
2587
2588 /* y = y0 + factor1 yd0 + factor2 * Y[0] + factor3 Y[1] */
2589 Xvecs[0] = y0;
2590 Xvecs[1] = yd0;
2591 Xvecs[2] = IDAADJ_mem->ia_Y[0];
2592 Xvecs[3] = IDAADJ_mem->ia_Y[1];
2593
2594 retval = N_VLinearCombination(4, cvals, Xvecs, yy);
2595 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2596
2597 /* Sensi Interpolation. */
2598
2599 /* yS = yS0 + factor1 ySd0 + factor2 * YS[0] + factor3 YS[1], if needed */
2600 if (NS > 0) {
2601
2602 XXvecs[0] = yS0;
2603 XXvecs[1] = ySd0;
2604 XXvecs[2] = IDAADJ_mem->ia_YS[0];
2605 XXvecs[3] = IDAADJ_mem->ia_YS[1];
2606
2607 retval = N_VLinearCombinationVectorArray(NS, 4, cvals, XXvecs, yyS);
2608 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2609
2610 }
2611
2612 /* For y'. */
2613 factor1 = factor1/delta/delta; /* factor1 = 2(t-t0)/(t1-t0)^2 */
2614 factor2 = factor1*((3*t-2*t1-t0)/delta); /* factor2 = (t-t0)(3*t-2*t1-t0)/(t1-t0)^3 */
2615 factor1 *= 2;
2616
2617 cvals[0] = ONE;
2618 cvals[1] = factor1;
2619 cvals[2] = factor2;
2620
2621 /* yp = yd0 + factor1 Y[0] + factor 2 Y[1] */
2622 Xvecs[0] = yd0;
2623 Xvecs[1] = IDAADJ_mem->ia_Y[0];
2624 Xvecs[2] = IDAADJ_mem->ia_Y[1];
2625
2626 retval = N_VLinearCombination(3, cvals, Xvecs, yp);
2627 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2628
2629 /* Sensi interpolation for 1st derivative. */
2630
2631 /* ypS = ySd0 + factor1 YS[0] + factor 2 YS[1], if needed */
2632 if (NS > 0) {
2633
2634 XXvecs[0] = ySd0;
2635 XXvecs[1] = IDAADJ_mem->ia_YS[0];
2636 XXvecs[2] = IDAADJ_mem->ia_YS[1];
2637
2638 retval = N_VLinearCombinationVectorArray(NS, 3, cvals, XXvecs, ypS);
2639 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2640
2641 }
2642
2643 return(IDA_SUCCESS);
2644 }
2645
2646 /*
2647 * -----------------------------------------------------------------
2648 * Functions specific to Polynomial interpolation
2649 * -----------------------------------------------------------------
2650 */
2651
2652 /*
2653 * IDAApolynomialMalloc
2654 *
2655 * This routine allocates memory for storing information at all
2656 * intermediate points between two consecutive check points.
2657 * This data is then used to interpolate the forward solution
2658 * at any other time.
2659 *
2660 * Information about the first derivative is stored only for the first
2661 * data point.
2662 */
2663
IDAApolynomialMalloc(IDAMem IDA_mem)2664 static booleantype IDAApolynomialMalloc(IDAMem IDA_mem)
2665 {
2666 IDAadjMem IDAADJ_mem;
2667 DtpntMem *dt_mem;
2668 PolynomialDataMem content;
2669 long int i, ii=0;
2670 booleantype allocOK;
2671
2672 allocOK = SUNTRUE;
2673
2674 IDAADJ_mem = IDA_mem->ida_adj_mem;
2675
2676 /* Allocate space for the vectors yyTmp and ypTmp */
2677 IDAADJ_mem->ia_yyTmp = N_VClone(IDA_mem->ida_tempv1);
2678 if (IDAADJ_mem->ia_yyTmp == NULL) {
2679 return(SUNFALSE);
2680 }
2681 IDAADJ_mem->ia_ypTmp = N_VClone(IDA_mem->ida_tempv1);
2682 if (IDAADJ_mem->ia_ypTmp == NULL) {
2683 return(SUNFALSE);
2684 }
2685
2686 if (IDAADJ_mem->ia_storeSensi) {
2687
2688 IDAADJ_mem->ia_yySTmp = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_tempv1);
2689 if (IDAADJ_mem->ia_yySTmp == NULL) {
2690 N_VDestroy(IDAADJ_mem->ia_yyTmp);
2691 N_VDestroy(IDAADJ_mem->ia_ypTmp);
2692 return(SUNFALSE);
2693 }
2694
2695 IDAADJ_mem->ia_ypSTmp = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_tempv1);
2696 if (IDAADJ_mem->ia_ypSTmp == NULL) {
2697 N_VDestroy(IDAADJ_mem->ia_yyTmp);
2698 N_VDestroy(IDAADJ_mem->ia_ypTmp);
2699 N_VDestroyVectorArray(IDAADJ_mem->ia_yySTmp, IDA_mem->ida_Ns);
2700 return(SUNFALSE);
2701
2702 }
2703 }
2704
2705 /* Allocate space for the content field of the dt structures */
2706 dt_mem = IDAADJ_mem->dt_mem;
2707
2708 for (i=0; i<=IDAADJ_mem->ia_nsteps; i++) {
2709
2710 content = NULL;
2711 content = (PolynomialDataMem) malloc(sizeof(struct PolynomialDataMemRec));
2712 if (content == NULL) {
2713 ii = i;
2714 allocOK = SUNFALSE;
2715 break;
2716 }
2717
2718 content->y = N_VClone(IDA_mem->ida_tempv1);
2719 if (content->y == NULL) {
2720 free(content); content = NULL;
2721 ii = i;
2722 allocOK = SUNFALSE;
2723 break;
2724 }
2725
2726 /* Allocate space for yp also. Needed for the most left point interpolation. */
2727 if (i == 0) {
2728 content->yd = N_VClone(IDA_mem->ida_tempv1);
2729
2730 /* Memory allocation failure ? */
2731 if (content->yd == NULL) {
2732 N_VDestroy(content->y);
2733 free(content); content = NULL;
2734 ii = i;
2735 allocOK = SUNFALSE;
2736 }
2737 } else {
2738 /* Not the first data point. */
2739 content->yd = NULL;
2740 }
2741
2742 if (IDAADJ_mem->ia_storeSensi) {
2743
2744 content->yS = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_tempv1);
2745 if (content->yS == NULL) {
2746 N_VDestroy(content->y);
2747 if (content->yd) N_VDestroy(content->yd);
2748 free(content); content = NULL;
2749 ii = i;
2750 allocOK = SUNFALSE;
2751 break;
2752 }
2753
2754 if (i==0) {
2755 content->ySd = N_VCloneVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_tempv1);
2756 if (content->ySd == NULL) {
2757 N_VDestroy(content->y);
2758 if (content->yd) N_VDestroy(content->yd);
2759 N_VDestroyVectorArray(content->yS, IDA_mem->ida_Ns);
2760 free(content); content = NULL;
2761 ii = i;
2762 allocOK = SUNFALSE;
2763 }
2764 } else {
2765 content->ySd = NULL;
2766 }
2767 }
2768
2769 dt_mem[i]->content = content;
2770 }
2771
2772 /* If an error occurred, deallocate and return */
2773 if (!allocOK) {
2774
2775 N_VDestroy(IDAADJ_mem->ia_yyTmp);
2776 N_VDestroy(IDAADJ_mem->ia_ypTmp);
2777 if (IDAADJ_mem->ia_storeSensi) {
2778
2779 N_VDestroyVectorArray(IDAADJ_mem->ia_yySTmp, IDA_mem->ida_Ns);
2780 N_VDestroyVectorArray(IDAADJ_mem->ia_ypSTmp, IDA_mem->ida_Ns);
2781 }
2782
2783 for (i=0; i<ii; i++) {
2784 content = (PolynomialDataMem) (dt_mem[i]->content);
2785 N_VDestroy(content->y);
2786
2787 if (content->yd) N_VDestroy(content->yd);
2788
2789 if (IDAADJ_mem->ia_storeSensi) {
2790
2791 N_VDestroyVectorArray(content->yS, IDA_mem->ida_Ns);
2792
2793 if (content->ySd)
2794 N_VDestroyVectorArray(content->ySd, IDA_mem->ida_Ns);
2795 }
2796 free(dt_mem[i]->content); dt_mem[i]->content = NULL;
2797 }
2798
2799 }
2800 return(allocOK);
2801 }
2802
2803 /*
2804 * IDAApolynomialFree
2805 *
2806 * This routine frees the memory allocated for data storage.
2807 */
2808
IDAApolynomialFree(IDAMem IDA_mem)2809 static void IDAApolynomialFree(IDAMem IDA_mem)
2810 {
2811 IDAadjMem IDAADJ_mem;
2812 DtpntMem *dt_mem;
2813 PolynomialDataMem content;
2814 long int i;
2815
2816 IDAADJ_mem = IDA_mem->ida_adj_mem;
2817
2818 N_VDestroy(IDAADJ_mem->ia_yyTmp);
2819 N_VDestroy(IDAADJ_mem->ia_ypTmp);
2820
2821 if (IDAADJ_mem->ia_storeSensi) {
2822 N_VDestroyVectorArray(IDAADJ_mem->ia_yySTmp, IDA_mem->ida_Ns);
2823 N_VDestroyVectorArray(IDAADJ_mem->ia_ypSTmp, IDA_mem->ida_Ns);
2824 }
2825
2826 dt_mem = IDAADJ_mem->dt_mem;
2827
2828 for (i=0; i<=IDAADJ_mem->ia_nsteps; i++) {
2829
2830 content = (PolynomialDataMem) (dt_mem[i]->content);
2831
2832 /* content might be NULL, if IDAAdjInit was called but IDASolveF was not. */
2833 if(content) {
2834 N_VDestroy(content->y);
2835
2836 if (content->yd) N_VDestroy(content->yd);
2837
2838 if (IDAADJ_mem->ia_storeSensi) {
2839
2840 N_VDestroyVectorArray(content->yS, IDA_mem->ida_Ns);
2841
2842 if (content->ySd)
2843 N_VDestroyVectorArray(content->ySd, IDA_mem->ida_Ns);
2844 }
2845 free(dt_mem[i]->content); dt_mem[i]->content = NULL;
2846 }
2847 }
2848 }
2849
2850 /*
2851 * IDAApolynomialStorePnt
2852 *
2853 * This routine stores a new point y in the structure d for use
2854 * in the Polynomial interpolation.
2855 *
2856 * Note that the time is already stored. Information about the
2857 * first derivative is available only for the first data point,
2858 * in which case content->yp is non-null.
2859 */
2860
IDAApolynomialStorePnt(IDAMem IDA_mem,DtpntMem d)2861 static int IDAApolynomialStorePnt(IDAMem IDA_mem, DtpntMem d)
2862 {
2863 IDAadjMem IDAADJ_mem;
2864 PolynomialDataMem content;
2865 int is, retval;
2866
2867 IDAADJ_mem = IDA_mem->ida_adj_mem;
2868 content = (PolynomialDataMem) d->content;
2869
2870 N_VScale(ONE, IDA_mem->ida_phi[0], content->y);
2871
2872 /* copy also the derivative for the first data point (in this case
2873 content->yp is non-null). */
2874 if (content->yd)
2875 IDAAGettnSolutionYp(IDA_mem, content->yd);
2876
2877 if (IDAADJ_mem->ia_storeSensi) {
2878
2879 for (is=0; is<IDA_mem->ida_Ns; is++)
2880 IDA_mem->ida_cvals[is] = ONE;
2881
2882 retval = N_VScaleVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_cvals,
2883 IDA_mem->ida_phiS[0], content->yS);
2884 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2885
2886 /* store the derivative if it is the first data point. */
2887 if(content->ySd)
2888 IDAAGettnSolutionYpS(IDA_mem, content->ySd);
2889 }
2890
2891 content->order = IDA_mem->ida_kused;
2892
2893 return(0);
2894 }
2895
2896 /*
2897 * IDAApolynomialGetY
2898 *
2899 * This routine uses polynomial interpolation for the forward solution vector.
2900 * It is typically called by the wrapper routines before calling
2901 * user provided routines (fB, djacB, bjacB, jtimesB, psolB)) but
2902 * can be directly called by the user through CVodeGetAdjY.
2903 */
2904
IDAApolynomialGetY(IDAMem IDA_mem,realtype t,N_Vector yy,N_Vector yp,N_Vector * yyS,N_Vector * ypS)2905 static int IDAApolynomialGetY(IDAMem IDA_mem, realtype t,
2906 N_Vector yy, N_Vector yp,
2907 N_Vector *yyS, N_Vector *ypS)
2908 {
2909 IDAadjMem IDAADJ_mem;
2910 DtpntMem *dt_mem;
2911 PolynomialDataMem content;
2912
2913 int flag, dir, order, i, j, is, NS, retval;
2914 long int indx, base;
2915 booleantype newpoint;
2916 realtype delt, factor, Psi, Psiprime;
2917
2918 IDAADJ_mem = IDA_mem->ida_adj_mem;
2919 dt_mem = IDAADJ_mem->dt_mem;
2920
2921 /* Local value of Ns */
2922 NS = (IDAADJ_mem->ia_interpSensi && (yyS != NULL)) ? IDA_mem->ida_Ns : 0;
2923
2924 /* Get the index in dt_mem */
2925 flag = IDAAfindIndex(IDA_mem, t, &indx, &newpoint);
2926 if (flag != IDA_SUCCESS) return(flag);
2927
2928 /* If we are beyond the left limit but close enough,
2929 then return y at the left limit. */
2930
2931 if (indx == 0) {
2932 content = (PolynomialDataMem) (dt_mem[0]->content);
2933 N_VScale(ONE, content->y, yy);
2934 N_VScale(ONE, content->yd, yp);
2935
2936 if (NS > 0) {
2937 for (is=0; is<NS; is++)
2938 IDA_mem->ida_cvals[is] = ONE;
2939
2940 retval = N_VScaleVectorArray(NS, IDA_mem->ida_cvals, content->yS, yyS);
2941 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2942
2943 retval = N_VScaleVectorArray(NS, IDA_mem->ida_cvals, content->ySd, ypS);
2944 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2945 }
2946
2947 return(IDA_SUCCESS);
2948 }
2949
2950 /* Scaling factor */
2951 delt = SUNRabs(dt_mem[indx]->t - dt_mem[indx-1]->t);
2952
2953 /* Find the direction of the forward integration */
2954 dir = (IDAADJ_mem->ia_tfinal - IDAADJ_mem->ia_tinitial > ZERO) ? 1 : -1;
2955
2956 /* Establish the base point depending on the integration direction.
2957 Modify the base if there are not enough points for the current order */
2958
2959 if (dir == 1) {
2960 base = indx;
2961 content = (PolynomialDataMem) (dt_mem[base]->content);
2962 order = content->order;
2963 if(indx < order) base += order-indx;
2964 } else {
2965 base = indx-1;
2966 content = (PolynomialDataMem) (dt_mem[base]->content);
2967 order = content->order;
2968 if (IDAADJ_mem->ia_np-indx > order) base -= indx+order-IDAADJ_mem->ia_np;
2969 }
2970
2971 /* Recompute Y (divided differences for Newton polynomial) if needed */
2972
2973 if (newpoint) {
2974
2975 /* Store 0-th order DD */
2976 if (dir == 1) {
2977 for(j=0;j<=order;j++) {
2978 IDAADJ_mem->ia_T[j] = dt_mem[base-j]->t;
2979 content = (PolynomialDataMem) (dt_mem[base-j]->content);
2980 N_VScale(ONE, content->y, IDAADJ_mem->ia_Y[j]);
2981
2982 if (NS > 0) {
2983 for (is=0; is<NS; is++)
2984 IDA_mem->ida_cvals[is] = ONE;
2985 retval = N_VScaleVectorArray(NS, IDA_mem->ida_cvals,
2986 content->yS, IDAADJ_mem->ia_YS[j]);
2987 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
2988 }
2989 }
2990 } else {
2991 for(j=0;j<=order;j++) {
2992 IDAADJ_mem->ia_T[j] = dt_mem[base-1+j]->t;
2993 content = (PolynomialDataMem) (dt_mem[base-1+j]->content);
2994 N_VScale(ONE, content->y, IDAADJ_mem->ia_Y[j]);
2995
2996 if (NS > 0) {
2997 for (is=0; is<NS; is++)
2998 IDA_mem->ida_cvals[is] = ONE;
2999 retval = N_VScaleVectorArray(NS, IDA_mem->ida_cvals,
3000 content->yS, IDAADJ_mem->ia_YS[j]);
3001 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
3002 }
3003 }
3004 }
3005
3006 /* Compute higher-order DD */
3007 for(i=1;i<=order;i++) {
3008 for(j=order;j>=i;j--) {
3009 factor = delt/(IDAADJ_mem->ia_T[j]-IDAADJ_mem->ia_T[j-i]);
3010 N_VLinearSum(factor, IDAADJ_mem->ia_Y[j], -factor, IDAADJ_mem->ia_Y[j-1], IDAADJ_mem->ia_Y[j]);
3011
3012 for (is=0; is<NS; is++)
3013 N_VLinearSum(factor, IDAADJ_mem->ia_YS[j][is], -factor, IDAADJ_mem->ia_YS[j-1][is], IDAADJ_mem->ia_YS[j][is]);
3014
3015 }
3016 }
3017 }
3018
3019 /* Perform the actual interpolation for yy using nested multiplications */
3020
3021 IDA_mem->ida_cvals[0] = ONE;
3022 for (i=0; i<order; i++)
3023 IDA_mem->ida_cvals[i+1] = IDA_mem->ida_cvals[i] * (t-IDAADJ_mem->ia_T[i]) / delt;
3024
3025 retval = N_VLinearCombination(order+1, IDA_mem->ida_cvals, IDAADJ_mem->ia_Y, yy);
3026 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
3027
3028 if (NS > 0) {
3029 retval = N_VLinearCombinationVectorArray(NS, order+1, IDA_mem->ida_cvals, IDAADJ_mem->ia_YS, yyS);
3030 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
3031 }
3032
3033 /* Perform the actual interpolation for yp.
3034
3035 Writing p(t) = y0 + (t-t0)*f[t0,t1] + ... + (t-t0)(t-t1)...(t-tn)*f[t0,t1,...tn],
3036 denote psi_k(t) = (t-t0)(t-t1)...(t-tk).
3037
3038 The formula used for p'(t) is:
3039 - p'(t) = f[t0,t1] + psi_1'(t)*f[t0,t1,t2] + ... + psi_n'(t)*f[t0,t1,...,tn]
3040
3041 We reccursively compute psi_k'(t) from:
3042 - psi_k'(t) = (t-tk)*psi_{k-1}'(t) + psi_{k-1}
3043
3044 psi_k is rescaled with 1/delt each time is computed, because the Newton DDs from Y were
3045 scaled with delt.
3046 */
3047
3048 Psi = ONE;
3049 Psiprime = ZERO;
3050
3051 for(i=1; i<=order; i++) {
3052 factor = (t-IDAADJ_mem->ia_T[i-1])/delt;
3053
3054 Psiprime = Psi/delt + factor * Psiprime;
3055 Psi = Psi * factor;
3056
3057 IDA_mem->ida_cvals[i-1] = Psiprime;
3058 }
3059
3060 retval = N_VLinearCombination(order, IDA_mem->ida_cvals, IDAADJ_mem->ia_Y+1, yp);
3061 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
3062
3063 if (NS > 0) {
3064 retval = N_VLinearCombinationVectorArray(NS, order, IDA_mem->ida_cvals, IDAADJ_mem->ia_YS+1, ypS);
3065 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
3066 }
3067
3068 return(IDA_SUCCESS);
3069 }
3070
3071 /*
3072 * IDAAGettnSolutionYp
3073 *
3074 * Evaluates the first derivative of the solution at the last time returned by
3075 * IDASolve (tretlast).
3076 *
3077 * The function implements the same algorithm as in IDAGetSolution but in the
3078 * particular case when t=tn (i.e. delta=0).
3079 *
3080 * This function was implemented to avoid calls to IDAGetSolution which computes
3081 * y by doing a loop that is not necessary for this particular situation.
3082 */
3083
IDAAGettnSolutionYp(IDAMem IDA_mem,N_Vector yp)3084 static int IDAAGettnSolutionYp(IDAMem IDA_mem, N_Vector yp)
3085 {
3086 int j, kord, retval;
3087 realtype C, D, gam;
3088
3089 if (IDA_mem->ida_nst==0) {
3090
3091 /* If no integration was done, return the yp supplied by user.*/
3092 N_VScale(ONE, IDA_mem->ida_phi[1], yp);
3093
3094 return(0);
3095 }
3096
3097 /* Compute yp as in IDAGetSolution for this particular case when t=tn. */
3098
3099 kord = IDA_mem->ida_kused;
3100 if(IDA_mem->ida_kused==0) kord=1;
3101
3102 C = ONE; D = ZERO;
3103 gam = ZERO;
3104 for (j=1; j <= kord; j++) {
3105 D = D*gam + C/IDA_mem->ida_psi[j-1];
3106 C = C*gam;
3107 gam = IDA_mem->ida_psi[j-1] / IDA_mem->ida_psi[j];
3108
3109 IDA_mem->ida_dvals[j-1] = D;
3110 }
3111
3112 retval = N_VLinearCombination(kord, IDA_mem->ida_dvals,
3113 IDA_mem->ida_phi+1, yp);
3114 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
3115
3116 return(0);
3117 }
3118
3119
3120 /*
3121 * IDAAGettnSolutionYpS
3122 *
3123 * Same as IDAAGettnSolutionYp, but for first derivative of the sensitivities.
3124 *
3125 */
3126
IDAAGettnSolutionYpS(IDAMem IDA_mem,N_Vector * ypS)3127 static int IDAAGettnSolutionYpS(IDAMem IDA_mem, N_Vector *ypS)
3128 {
3129 int j, kord, is, retval;
3130 realtype C, D, gam;
3131
3132 if (IDA_mem->ida_nst==0) {
3133
3134 /* If no integration was done, return the ypS supplied by user.*/
3135 for (is=0; is<IDA_mem->ida_Ns; is++)
3136 IDA_mem->ida_cvals[is] = ONE;
3137
3138 retval = N_VScaleVectorArray(IDA_mem->ida_Ns, IDA_mem->ida_cvals,
3139 IDA_mem->ida_phiS[1], ypS);
3140 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
3141
3142 return(0);
3143 }
3144
3145 kord = IDA_mem->ida_kused;
3146 if(IDA_mem->ida_kused==0) kord=1;
3147
3148 C = ONE; D = ZERO;
3149 gam = ZERO;
3150 for (j=1; j <= kord; j++) {
3151 D = D*gam + C/IDA_mem->ida_psi[j-1];
3152 C = C*gam;
3153 gam = IDA_mem->ida_psi[j-1] / IDA_mem->ida_psi[j];
3154
3155 IDA_mem->ida_dvals[j-1] = D;
3156 }
3157
3158 retval = N_VLinearCombinationVectorArray(IDA_mem->ida_Ns, kord,
3159 IDA_mem->ida_dvals,
3160 IDA_mem->ida_phiS+1, ypS);
3161 if (retval != IDA_SUCCESS) return (IDA_VECTOROP_ERR);
3162
3163 return(0);
3164 }
3165
3166
3167
3168 /*
3169 * IDAAfindIndex
3170 *
3171 * Finds the index in the array of data point strctures such that
3172 * dt_mem[indx-1].t <= t < dt_mem[indx].t
3173 * If indx is changed from the previous invocation, then newpoint = SUNTRUE
3174 *
3175 * If t is beyond the leftmost limit, but close enough, indx=0.
3176 *
3177 * Returns IDA_SUCCESS if successful and IDA_GETY_BADT if unable to
3178 * find indx (t is too far beyond limits).
3179 */
3180
IDAAfindIndex(IDAMem ida_mem,realtype t,long int * indx,booleantype * newpoint)3181 static int IDAAfindIndex(IDAMem ida_mem, realtype t,
3182 long int *indx, booleantype *newpoint)
3183 {
3184 IDAadjMem IDAADJ_mem;
3185 IDAMem IDA_mem;
3186 DtpntMem *dt_mem;
3187 int sign;
3188 booleantype to_left, to_right;
3189
3190 IDA_mem = (IDAMem) ida_mem;
3191 IDAADJ_mem = IDA_mem->ida_adj_mem;
3192 dt_mem = IDAADJ_mem->dt_mem;
3193
3194 *newpoint = SUNFALSE;
3195
3196 /* Find the direction of integration */
3197 sign = (IDAADJ_mem->ia_tfinal - IDAADJ_mem->ia_tinitial > ZERO) ? 1 : -1;
3198
3199 /* If this is the first time we use new data */
3200 if (IDAADJ_mem->ia_newData) {
3201 IDAADJ_mem->ia_ilast = IDAADJ_mem->ia_np-1;
3202 *newpoint = SUNTRUE;
3203 IDAADJ_mem->ia_newData = SUNFALSE;
3204 }
3205
3206 /* Search for indx starting from ilast */
3207 to_left = ( sign*(t - dt_mem[IDAADJ_mem->ia_ilast-1]->t) < ZERO);
3208 to_right = ( sign*(t - dt_mem[IDAADJ_mem->ia_ilast]->t) > ZERO);
3209
3210 if ( to_left ) {
3211 /* look for a new indx to the left */
3212
3213 *newpoint = SUNTRUE;
3214
3215 *indx = IDAADJ_mem->ia_ilast;
3216 for(;;) {
3217 if ( *indx == 0 ) break;
3218 if ( sign*(t - dt_mem[*indx-1]->t) <= ZERO ) (*indx)--;
3219 else break;
3220 }
3221
3222 if ( *indx == 0 )
3223 IDAADJ_mem->ia_ilast = 1;
3224 else
3225 IDAADJ_mem->ia_ilast = *indx;
3226
3227 if ( *indx == 0 ) {
3228 /* t is beyond leftmost limit. Is it too far? */
3229 if ( SUNRabs(t - dt_mem[0]->t) > FUZZ_FACTOR * IDA_mem->ida_uround ) {
3230 return(IDA_GETY_BADT);
3231 }
3232 }
3233
3234 } else if ( to_right ) {
3235 /* look for a new indx to the right */
3236
3237 *newpoint = SUNTRUE;
3238
3239 *indx = IDAADJ_mem->ia_ilast;
3240 for(;;) {
3241 if ( sign*(t - dt_mem[*indx]->t) > ZERO) (*indx)++;
3242 else break;
3243 }
3244
3245 IDAADJ_mem->ia_ilast = *indx;
3246
3247 } else {
3248 /* ilast is still OK */
3249
3250 *indx = IDAADJ_mem->ia_ilast;
3251
3252 }
3253 return(IDA_SUCCESS);
3254 }
3255
3256
3257 /*
3258 * IDAGetAdjY
3259 *
3260 * This routine returns the interpolated forward solution at time t.
3261 * The user must allocate space for y.
3262 */
3263
IDAGetAdjY(void * ida_mem,realtype t,N_Vector yy,N_Vector yp)3264 int IDAGetAdjY(void *ida_mem, realtype t, N_Vector yy, N_Vector yp)
3265 {
3266 IDAMem IDA_mem;
3267 IDAadjMem IDAADJ_mem;
3268 int flag;
3269
3270 if (ida_mem == NULL) {
3271 IDAProcessError(NULL, IDA_MEM_NULL, "IDAA", "IDAGetAdjY", MSG_NO_MEM);
3272 return(IDA_MEM_NULL);
3273 }
3274 IDA_mem = (IDAMem) ida_mem;
3275 IDAADJ_mem = IDA_mem->ida_adj_mem;
3276
3277 flag = IDAADJ_mem->ia_getY(IDA_mem, t, yy, yp, NULL, NULL);
3278
3279 return(flag);
3280 }
3281
3282 /*=================================================================*/
3283 /* Wrappers for adjoint system */
3284 /*=================================================================*/
3285
3286 /*
3287 * IDAAres
3288 *
3289 * This routine interfaces to the RhsFnB routine provided by
3290 * the user.
3291 */
3292
IDAAres(realtype tt,N_Vector yyB,N_Vector ypB,N_Vector rrB,void * ida_mem)3293 static int IDAAres(realtype tt,
3294 N_Vector yyB, N_Vector ypB, N_Vector rrB,
3295 void *ida_mem)
3296 {
3297 IDAadjMem IDAADJ_mem;
3298 IDABMem IDAB_mem;
3299 IDAMem IDA_mem;
3300 int flag, retval;
3301
3302 IDA_mem = (IDAMem) ida_mem;
3303
3304 IDAADJ_mem = IDA_mem->ida_adj_mem;
3305
3306 /* Get the current backward problem. */
3307 IDAB_mem = IDAADJ_mem->ia_bckpbCrt;
3308
3309 /* Get forward solution from interpolation. */
3310 if( IDAADJ_mem->ia_noInterp == SUNFALSE) {
3311 if (IDAADJ_mem->ia_interpSensi)
3312 flag = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp, IDAADJ_mem->ia_ypTmp, IDAADJ_mem->ia_yySTmp, IDAADJ_mem->ia_ypSTmp);
3313 else
3314 flag = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp, IDAADJ_mem->ia_ypTmp, NULL, NULL);
3315
3316 if (flag != IDA_SUCCESS) {
3317 IDAProcessError(IDA_mem, -1, "IDAA", "IDAAres", MSGAM_BAD_TINTERP, tt);
3318 return(-1);
3319 }
3320 }
3321
3322 /* Call the user supplied residual. */
3323 if(IDAB_mem->ida_res_withSensi) {
3324 retval = IDAB_mem->ida_resS(tt, IDAADJ_mem->ia_yyTmp, IDAADJ_mem->ia_ypTmp,
3325 IDAADJ_mem->ia_yySTmp, IDAADJ_mem->ia_ypSTmp,
3326 yyB, ypB,
3327 rrB, IDAB_mem->ida_user_data);
3328 }else {
3329 retval = IDAB_mem->ida_res(tt, IDAADJ_mem->ia_yyTmp, IDAADJ_mem->ia_ypTmp, yyB, ypB, rrB, IDAB_mem->ida_user_data);
3330 }
3331 return(retval);
3332 }
3333
3334 /*
3335 *IDAArhsQ
3336 *
3337 * This routine interfaces to the IDAQuadRhsFnB routine provided by
3338 * the user.
3339 *
3340 * It is passed to IDAQuadInit calls for backward problem, so it must
3341 * be of IDAQuadRhsFn type.
3342 */
3343
IDAArhsQ(realtype tt,N_Vector yyB,N_Vector ypB,N_Vector resvalQB,void * ida_mem)3344 static int IDAArhsQ(realtype tt,
3345 N_Vector yyB, N_Vector ypB,
3346 N_Vector resvalQB, void *ida_mem)
3347 {
3348 IDAMem IDA_mem;
3349 IDAadjMem IDAADJ_mem;
3350 IDABMem IDAB_mem;
3351 int retval, flag;
3352
3353 IDA_mem = (IDAMem) ida_mem;
3354 IDAADJ_mem = IDA_mem->ida_adj_mem;
3355
3356 /* Get current backward problem. */
3357 IDAB_mem = IDAADJ_mem->ia_bckpbCrt;
3358
3359 retval = IDA_SUCCESS;
3360
3361 /* Get forward solution from interpolation. */
3362 if (IDAADJ_mem->ia_noInterp == SUNFALSE) {
3363 if (IDAADJ_mem->ia_interpSensi) {
3364 flag = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp, IDAADJ_mem->ia_ypTmp, IDAADJ_mem->ia_yySTmp, IDAADJ_mem->ia_ypSTmp);
3365 } else {
3366 flag = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp, IDAADJ_mem->ia_ypTmp, NULL, NULL);
3367 }
3368
3369 if (flag != IDA_SUCCESS) {
3370 IDAProcessError(IDA_mem, -1, "IDAA", "IDAArhsQ", MSGAM_BAD_TINTERP, tt);
3371 return(-1);
3372 }
3373 }
3374
3375 /* Call user's adjoint quadrature RHS routine */
3376 if (IDAB_mem->ida_rhsQ_withSensi) {
3377 retval = IDAB_mem->ida_rhsQS(tt, IDAADJ_mem->ia_yyTmp, IDAADJ_mem->ia_ypTmp, IDAADJ_mem->ia_yySTmp, IDAADJ_mem->ia_ypSTmp,
3378 yyB, ypB,
3379 resvalQB, IDAB_mem->ida_user_data);
3380 } else {
3381 retval = IDAB_mem->ida_rhsQ(tt,
3382 IDAADJ_mem->ia_yyTmp, IDAADJ_mem->ia_ypTmp,
3383 yyB, ypB,
3384 resvalQB, IDAB_mem->ida_user_data);
3385 }
3386 return(retval);
3387 }
3388