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