1 /* -----------------------------------------------------------------
2 * Programmer(s): Alan Hindmarsh, Radu Serban and Aaron Collier @ LLNL
3 * -----------------------------------------------------------------
4 * SUNDIALS Copyright Start
5 * Copyright (c) 2002-2020, Lawrence Livermore National Security
6 * and Southern Methodist University.
7 * All rights reserved.
8 *
9 * See the top-level LICENSE and NOTICE files for details.
10 *
11 * SPDX-License-Identifier: BSD-3-Clause
12 * SUNDIALS Copyright End
13 * -----------------------------------------------------------------
14 * This is the implementation file for the main IDA solver.
15 * -----------------------------------------------------------------
16 *
17 * EXPORTED FUNCTIONS
18 * ------------------
19 * Creation, allocation and re-initialization functions
20 * IDACreate
21 * IDAInit
22 * IDAReInit
23 * IDARootInit
24 * Main solver function
25 * IDASolve
26 * Interpolated output and extraction functions
27 * IDAGetDky
28 * Deallocation functions
29 * IDAFree
30 *
31 * PRIVATE FUNCTIONS
32 * -----------------
33 * IDACheckNvector
34 * Memory allocation/deallocation
35 * IDAAllocVectors
36 * IDAFreeVectors
37 * Initial setup
38 * IDAInitialSetup
39 * IDAEwtSet
40 * IDAEwtSetSS
41 * IDAEwtSetSV
42 * Stopping tests
43 * IDAStopTest1
44 * IDAStopTest2
45 * Error handler
46 * IDAHandleFailure
47 * Main IDAStep function
48 * IDAStep
49 * IDASetCoeffs
50 * Nonlinear solver functions
51 * IDANls
52 * IDAPredict
53 * Error test
54 * IDATestError
55 * IDARestore
56 * Handler for convergence and/or error test failures
57 * IDAHandleNFlag
58 * IDAReset
59 * Function called after a successful step
60 * IDACompleteStep
61 * Get solution
62 * IDAGetSolution
63 * Norm functions
64 * IDAWrmsNorm
65 * Functions for rootfinding
66 * IDARcheck1
67 * IDARcheck2
68 * IDARcheck3
69 * IDARootfind
70 * IDA Error message handling functions
71 * IDAProcessError
72 * IDAErrHandler
73 * -----------------------------------------------------------------
74 */
75
76 /*
77 * =================================================================
78 * IMPORTED HEADER FILES
79 * =================================================================
80 */
81
82 #include <stdio.h>
83 #include <stdlib.h>
84 #include <stdarg.h>
85 #include <string.h>
86
87 #include "ida_impl.h"
88 #include <sundials/sundials_math.h>
89 #include <sunnonlinsol/sunnonlinsol_newton.h>
90
91 /*
92 * =================================================================
93 * IDAS PRIVATE CONSTANTS
94 * =================================================================
95 */
96
97 #define ZERO RCONST(0.0) /* real 0.0 */
98 #define HALF RCONST(0.5) /* real 0.5 */
99 #define QUARTER RCONST(0.25) /* real 0.25 */
100 #define TWOTHIRDS RCONST(0.667) /* real 2/3 */
101 #define ONE RCONST(1.0) /* real 1.0 */
102 #define ONEPT5 RCONST(1.5) /* real 1.5 */
103 #define TWO RCONST(2.0) /* real 2.0 */
104 #define FOUR RCONST(4.0) /* real 4.0 */
105 #define FIVE RCONST(5.0) /* real 5.0 */
106 #define TEN RCONST(10.0) /* real 10.0 */
107 #define TWELVE RCONST(12.0) /* real 12.0 */
108 #define TWENTY RCONST(20.0) /* real 20.0 */
109 #define HUNDRED RCONST(100.0) /* real 100.0 */
110 #define PT9 RCONST(0.9) /* real 0.9 */
111 #define PT99 RCONST(0.99) /* real 0.99 */
112 #define PT1 RCONST(0.1) /* real 0.1 */
113 #define PT01 RCONST(0.01) /* real 0.01 */
114 #define PT001 RCONST(0.001) /* real 0.001 */
115 #define PT0001 RCONST(0.0001) /* real 0.0001 */
116
117 /*
118 * =================================================================
119 * IDAS ROUTINE-SPECIFIC CONSTANTS
120 * =================================================================
121 */
122
123 /*
124 * Control constants for lower-level functions used by IDASolve
125 * ------------------------------------------------------------
126 */
127
128 /* IDAStep control constants */
129
130 #define PREDICT_AGAIN 20
131
132 /* Return values for lower level routines used by IDASolve */
133
134 #define CONTINUE_STEPS +99
135
136 /* IDACompleteStep constants */
137
138 #define UNSET -1
139 #define LOWER +1
140 #define RAISE +2
141 #define MAINTAIN +3
142
143 /* IDATestError constants */
144
145 #define ERROR_TEST_FAIL +7
146
147 /*
148 * Control constants for lower-level rootfinding functions
149 * -------------------------------------------------------
150 */
151
152 #define RTFOUND +1
153 #define CLOSERT +3
154
155 /*
156 * Control constants for tolerances
157 * --------------------------------
158 */
159
160 #define IDA_NN 0
161 #define IDA_SS 1
162 #define IDA_SV 2
163 #define IDA_WF 3
164
165 /*
166 * Algorithmic constants
167 * ---------------------
168 */
169
170 #define MXNCF 10 /* max number of convergence failures allowed */
171 #define MXNEF 10 /* max number of error test failures allowed */
172 #define MAXNH 5 /* max. number of h tries in IC calc. */
173 #define MAXNJ 4 /* max. number of J tries in IC calc. */
174 #define MAXNI 10 /* max. Newton iterations in IC calc. */
175 #define EPCON RCONST(0.33) /* Newton convergence test constant */
176 #define MAXBACKS 100 /* max backtracks per Newton step in IDACalcIC */
177 #define XRATE RCONST(0.25) /* constant for updating Jacobian/preconditioner */
178
179 /*
180 * =================================================================
181 * PRIVATE FUNCTION PROTOTYPES
182 * =================================================================
183 */
184
185 static booleantype IDACheckNvector(N_Vector tmpl);
186
187 /* Memory allocation/deallocation */
188
189 static booleantype IDAAllocVectors(IDAMem IDA_mem, N_Vector tmpl);
190 static void IDAFreeVectors(IDAMem IDA_mem);
191
192 /* Initial setup */
193
194 int IDAInitialSetup(IDAMem IDA_mem);
195
196 static int IDAEwtSetSS(IDAMem IDA_mem, N_Vector ycur, N_Vector weight);
197 static int IDAEwtSetSV(IDAMem IDA_mem, N_Vector ycur, N_Vector weight);
198
199 /* Main IDAStep function */
200
201 static int IDAStep(IDAMem IDA_mem);
202
203 /* Function called at beginning of step */
204
205 static void IDASetCoeffs(IDAMem IDA_mem, realtype *ck);
206
207 /* Nonlinear solver functions */
208
209 static void IDAPredict(IDAMem IDA_mem);
210 static int IDANls(IDAMem IDA_mem);
211
212 /* Error test */
213
214 static int IDATestError(IDAMem IDA_mem, realtype ck,
215 realtype *err_k, realtype *err_km1);
216
217 /* Handling of convergence and/or error test failures */
218
219 static void IDARestore(IDAMem IDA_mem, realtype saved_t);
220 static int IDAHandleNFlag(IDAMem IDA_mem, int nflag, realtype err_k, realtype err_km1,
221 long int *ncfnPtr, int *ncfPtr, long int *netfPtr, int *nefPtr);
222 static void IDAReset(IDAMem IDA_mem);
223
224 /* Function called after a successful step */
225
226 static void IDACompleteStep(IDAMem IDA_mem, realtype err_k, realtype err_km1);
227
228 /* Function called to evaluate the solutions y(t) and y'(t) at t */
229
230 int IDAGetSolution(void *ida_mem, realtype t, N_Vector yret, N_Vector ypret);
231
232 /* Stopping tests and failure handling */
233
234 static int IDAStopTest1(IDAMem IDA_mem, realtype tout,realtype *tret,
235 N_Vector yret, N_Vector ypret, int itask);
236 static int IDAStopTest2(IDAMem IDA_mem, realtype tout, realtype *tret,
237 N_Vector yret, N_Vector ypret, int itask);
238 static int IDAHandleFailure(IDAMem IDA_mem, int sflag);
239
240 /* Functions for rootfinding */
241
242 static int IDARcheck1(IDAMem IDA_mem);
243 static int IDARcheck2(IDAMem IDA_mem);
244 static int IDARcheck3(IDAMem IDA_mem);
245 static int IDARootfind(IDAMem IDA_mem);
246
247 /*
248 * =================================================================
249 * EXPORTED FUNCTIONS IMPLEMENTATION
250 * =================================================================
251 */
252
253 /*
254 * -----------------------------------------------------------------
255 * Creation, allocation and re-initialization functions
256 * -----------------------------------------------------------------
257 */
258
259 /*
260 * IDACreate
261 *
262 * IDACreate creates an internal memory block for a problem to
263 * be solved by IDA.
264 * If successful, IDACreate returns a pointer to the problem memory.
265 * This pointer should be passed to IDAInit.
266 * If an initialization error occurs, IDACreate prints an error
267 * message to standard err and returns NULL.
268 */
269
IDACreate(void)270 void *IDACreate(void)
271 {
272 IDAMem IDA_mem;
273
274 IDA_mem = NULL;
275 IDA_mem = (IDAMem) malloc(sizeof(struct IDAMemRec));
276 if (IDA_mem == NULL) {
277 IDAProcessError(NULL, 0, "IDA", "IDACreate", MSG_MEM_FAIL);
278 return (NULL);
279 }
280
281 /* Zero out ida_mem */
282 memset(IDA_mem, 0, sizeof(struct IDAMemRec));
283
284 /* Set unit roundoff in IDA_mem */
285 IDA_mem->ida_uround = UNIT_ROUNDOFF;
286
287 /* Set default values for integrator optional inputs */
288 IDA_mem->ida_res = NULL;
289 IDA_mem->ida_user_data = NULL;
290 IDA_mem->ida_itol = IDA_NN;
291 IDA_mem->ida_atolmin0 = SUNTRUE;
292 IDA_mem->ida_user_efun = SUNFALSE;
293 IDA_mem->ida_efun = NULL;
294 IDA_mem->ida_edata = NULL;
295 IDA_mem->ida_ehfun = IDAErrHandler;
296 IDA_mem->ida_eh_data = IDA_mem;
297 IDA_mem->ida_errfp = stderr;
298 IDA_mem->ida_maxord = MAXORD_DEFAULT;
299 IDA_mem->ida_mxstep = MXSTEP_DEFAULT;
300 IDA_mem->ida_hmax_inv = HMAX_INV_DEFAULT;
301 IDA_mem->ida_hin = ZERO;
302 IDA_mem->ida_epcon = EPCON;
303 IDA_mem->ida_maxnef = MXNEF;
304 IDA_mem->ida_maxncf = MXNCF;
305 IDA_mem->ida_suppressalg = SUNFALSE;
306 IDA_mem->ida_id = NULL;
307 IDA_mem->ida_constraints = NULL;
308 IDA_mem->ida_constraintsSet = SUNFALSE;
309 IDA_mem->ida_tstopset = SUNFALSE;
310
311 /* set the saved value maxord_alloc */
312 IDA_mem->ida_maxord_alloc = MAXORD_DEFAULT;
313
314 /* Set default values for IC optional inputs */
315 IDA_mem->ida_epiccon = PT01 * EPCON;
316 IDA_mem->ida_maxnh = MAXNH;
317 IDA_mem->ida_maxnj = MAXNJ;
318 IDA_mem->ida_maxnit = MAXNI;
319 IDA_mem->ida_maxbacks = MAXBACKS;
320 IDA_mem->ida_lsoff = SUNFALSE;
321 IDA_mem->ida_steptol = SUNRpowerR(IDA_mem->ida_uround, TWOTHIRDS);
322
323 /* Initialize lrw and liw */
324 IDA_mem->ida_lrw = 25 + 5*MXORDP1;
325 IDA_mem->ida_liw = 38;
326
327 /* No mallocs have been done yet */
328 IDA_mem->ida_VatolMallocDone = SUNFALSE;
329 IDA_mem->ida_constraintsMallocDone = SUNFALSE;
330 IDA_mem->ida_idMallocDone = SUNFALSE;
331 IDA_mem->ida_MallocDone = SUNFALSE;
332
333 /* Initialize nonlinear solver variables */
334 IDA_mem->NLS = NULL;
335 IDA_mem->ownNLS = SUNFALSE;
336
337 /* Return pointer to IDA memory block */
338 return((void *)IDA_mem);
339 }
340
341 /*-----------------------------------------------------------------*/
342
343 /*
344 * IDAInit
345 *
346 * IDAInit allocates and initializes memory for a problem. All
347 * problem specification inputs are checked for errors. If any
348 * error occurs during initialization, it is reported to the
349 * error handler function.
350 */
351
IDAInit(void * ida_mem,IDAResFn res,realtype t0,N_Vector yy0,N_Vector yp0)352 int IDAInit(void *ida_mem, IDAResFn res,
353 realtype t0, N_Vector yy0, N_Vector yp0)
354 {
355 int retval;
356 IDAMem IDA_mem;
357 booleantype nvectorOK, allocOK;
358 sunindextype lrw1, liw1;
359 SUNNonlinearSolver NLS;
360
361 /* Check ida_mem */
362
363 if (ida_mem == NULL) {
364 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDAInit", MSG_NO_MEM);
365 return(IDA_MEM_NULL);
366 }
367 IDA_mem = (IDAMem) ida_mem;
368
369 /* Check for legal input parameters */
370
371 if (yy0 == NULL) {
372 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAInit", MSG_Y0_NULL);
373 return(IDA_ILL_INPUT);
374 }
375
376 if (yp0 == NULL) {
377 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAInit", MSG_YP0_NULL);
378 return(IDA_ILL_INPUT);
379 }
380
381 if (res == NULL) {
382 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAInit", MSG_RES_NULL);
383 return(IDA_ILL_INPUT);
384 }
385
386 /* Test if all required vector operations are implemented */
387
388 nvectorOK = IDACheckNvector(yy0);
389 if (!nvectorOK) {
390 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAInit", MSG_BAD_NVECTOR);
391 return(IDA_ILL_INPUT);
392 }
393
394 /* Set space requirements for one N_Vector */
395
396 if (yy0->ops->nvspace != NULL) {
397 N_VSpace(yy0, &lrw1, &liw1);
398 } else {
399 lrw1 = 0;
400 liw1 = 0;
401 }
402 IDA_mem->ida_lrw1 = lrw1;
403 IDA_mem->ida_liw1 = liw1;
404
405 /* Allocate the vectors (using yy0 as a template) */
406
407 allocOK = IDAAllocVectors(IDA_mem, yy0);
408 if (!allocOK) {
409 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDAInit", MSG_MEM_FAIL);
410 return(IDA_MEM_FAIL);
411 }
412
413 /* create a Newton nonlinear solver object by default */
414 NLS = SUNNonlinSol_Newton(yy0);
415
416 /* check that nonlinear solver is non-NULL */
417 if (NLS == NULL) {
418 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDAInit", MSG_MEM_FAIL);
419 IDAFreeVectors(IDA_mem);
420 return(IDA_MEM_FAIL);
421 }
422
423 /* attach the nonlinear solver to the IDA memory */
424 retval = IDASetNonlinearSolver(IDA_mem, NLS);
425
426 /* check that the nonlinear solver was successfully attached */
427 if (retval != IDA_SUCCESS) {
428 IDAProcessError(IDA_mem, retval, "IDA", "IDAInit",
429 "Setting the nonlinear solver failed");
430 IDAFreeVectors(IDA_mem);
431 SUNNonlinSolFree(NLS);
432 return(IDA_MEM_FAIL);
433 }
434
435 /* set ownership flag */
436 IDA_mem->ownNLS = SUNTRUE;
437
438 /* All error checking is complete at this point */
439
440 /* Copy the input parameters into IDA memory block */
441
442 IDA_mem->ida_res = res;
443 IDA_mem->ida_tn = t0;
444
445 /* Set the linear solver addresses to NULL */
446
447 IDA_mem->ida_linit = NULL;
448 IDA_mem->ida_lsetup = NULL;
449 IDA_mem->ida_lsolve = NULL;
450 IDA_mem->ida_lperf = NULL;
451 IDA_mem->ida_lfree = NULL;
452 IDA_mem->ida_lmem = NULL;
453
454 /* Initialize the phi array */
455
456 N_VScale(ONE, yy0, IDA_mem->ida_phi[0]);
457 N_VScale(ONE, yp0, IDA_mem->ida_phi[1]);
458
459 /* Initialize all the counters and other optional output values */
460
461 IDA_mem->ida_nst = 0;
462 IDA_mem->ida_nre = 0;
463 IDA_mem->ida_ncfn = 0;
464 IDA_mem->ida_netf = 0;
465 IDA_mem->ida_nni = 0;
466 IDA_mem->ida_nsetups = 0;
467
468 IDA_mem->ida_kused = 0;
469 IDA_mem->ida_hused = ZERO;
470 IDA_mem->ida_tolsf = ONE;
471
472 IDA_mem->ida_nge = 0;
473
474 IDA_mem->ida_irfnd = 0;
475
476 /* Initialize root-finding variables */
477
478 IDA_mem->ida_glo = NULL;
479 IDA_mem->ida_ghi = NULL;
480 IDA_mem->ida_grout = NULL;
481 IDA_mem->ida_iroots = NULL;
482 IDA_mem->ida_rootdir = NULL;
483 IDA_mem->ida_gfun = NULL;
484 IDA_mem->ida_nrtfn = 0;
485 IDA_mem->ida_gactive = NULL;
486 IDA_mem->ida_mxgnull = 1;
487
488 /* Initial setup not done yet */
489
490 IDA_mem->ida_SetupDone = SUNFALSE;
491
492 /* Problem memory has been successfully allocated */
493
494 IDA_mem->ida_MallocDone = SUNTRUE;
495
496 return(IDA_SUCCESS);
497 }
498
499 /*-----------------------------------------------------------------*/
500
501 /*
502 * IDAReInit
503 *
504 * IDAReInit re-initializes IDA's memory for a problem, assuming
505 * it has already beeen allocated in a prior IDAInit call.
506 * All problem specification inputs are checked for errors.
507 * The problem size Neq is assumed to be unchanged since the call
508 * to IDAInit, and the maximum order maxord must not be larger.
509 * If any error occurs during reinitialization, it is reported to
510 * the error handler function.
511 * The return value is IDA_SUCCESS = 0 if no errors occurred, or
512 * a negative value otherwise.
513 */
514
IDAReInit(void * ida_mem,realtype t0,N_Vector yy0,N_Vector yp0)515 int IDAReInit(void *ida_mem,
516 realtype t0, N_Vector yy0, N_Vector yp0)
517 {
518 IDAMem IDA_mem;
519
520 /* Check for legal input parameters */
521
522 if (ida_mem == NULL) {
523 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDAReInit", MSG_NO_MEM);
524 return(IDA_MEM_NULL);
525 }
526 IDA_mem = (IDAMem) ida_mem;
527
528 /* Check if problem was malloc'ed */
529
530 if (IDA_mem->ida_MallocDone == SUNFALSE) {
531 IDAProcessError(IDA_mem, IDA_NO_MALLOC, "IDA", "IDAReInit", MSG_NO_MALLOC);
532 return(IDA_NO_MALLOC);
533 }
534
535 /* Check for legal input parameters */
536
537 if (yy0 == NULL) {
538 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAReInit", MSG_Y0_NULL);
539 return(IDA_ILL_INPUT);
540 }
541
542 if (yp0 == NULL) {
543 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAReInit", MSG_YP0_NULL);
544 return(IDA_ILL_INPUT);
545 }
546
547 /* Copy the input parameters into IDA memory block */
548
549 IDA_mem->ida_tn = t0;
550
551 /* Initialize the phi array */
552
553 N_VScale(ONE, yy0, IDA_mem->ida_phi[0]);
554 N_VScale(ONE, yp0, IDA_mem->ida_phi[1]);
555
556 /* Initialize all the counters and other optional output values */
557
558 IDA_mem->ida_nst = 0;
559 IDA_mem->ida_nre = 0;
560 IDA_mem->ida_ncfn = 0;
561 IDA_mem->ida_netf = 0;
562 IDA_mem->ida_nni = 0;
563 IDA_mem->ida_nsetups = 0;
564
565 IDA_mem->ida_kused = 0;
566 IDA_mem->ida_hused = ZERO;
567 IDA_mem->ida_tolsf = ONE;
568
569 IDA_mem->ida_nge = 0;
570
571 IDA_mem->ida_irfnd = 0;
572
573 /* Initial setup not done yet */
574
575 IDA_mem->ida_SetupDone = SUNFALSE;
576
577 /* Problem has been successfully re-initialized */
578
579 return(IDA_SUCCESS);
580 }
581
582 /*-----------------------------------------------------------------*/
583
584 /*
585 * IDASStolerances
586 * IDASVtolerances
587 * IDAWFtolerances
588 *
589 * These functions specify the integration tolerances. One of them
590 * MUST be called before the first call to IDA.
591 *
592 * IDASStolerances specifies scalar relative and absolute tolerances.
593 * IDASVtolerances specifies scalar relative tolerance and a vector
594 * absolute tolerance (a potentially different absolute tolerance
595 * for each vector component).
596 * IDAWFtolerances specifies a user-provides function (of type IDAEwtFn)
597 * which will be called to set the error weight vector.
598 */
599
IDASStolerances(void * ida_mem,realtype reltol,realtype abstol)600 int IDASStolerances(void *ida_mem, realtype reltol, realtype abstol)
601 {
602 IDAMem IDA_mem;
603
604 if (ida_mem==NULL) {
605 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDASStolerances", MSG_NO_MEM);
606 return(IDA_MEM_NULL);
607 }
608 IDA_mem = (IDAMem) ida_mem;
609
610 if (IDA_mem->ida_MallocDone == SUNFALSE) {
611 IDAProcessError(IDA_mem, IDA_NO_MALLOC, "IDA", "IDASStolerances", MSG_NO_MALLOC);
612 return(IDA_NO_MALLOC);
613 }
614
615 /* Check inputs */
616
617 if (reltol < ZERO) {
618 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASStolerances", MSG_BAD_RTOL);
619 return(IDA_ILL_INPUT);
620 }
621
622 if (abstol < ZERO) {
623 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASStolerances", MSG_BAD_ATOL);
624 return(IDA_ILL_INPUT);
625 }
626
627 /* Copy tolerances into memory */
628
629 IDA_mem->ida_rtol = reltol;
630 IDA_mem->ida_Satol = abstol;
631 IDA_mem->ida_atolmin0 = (abstol == ZERO);
632
633 IDA_mem->ida_itol = IDA_SS;
634
635 IDA_mem->ida_user_efun = SUNFALSE;
636 IDA_mem->ida_efun = IDAEwtSet;
637 IDA_mem->ida_edata = NULL; /* will be set to ida_mem in InitialSetup */
638
639 return(IDA_SUCCESS);
640 }
641
642
IDASVtolerances(void * ida_mem,realtype reltol,N_Vector abstol)643 int IDASVtolerances(void *ida_mem, realtype reltol, N_Vector abstol)
644 {
645 IDAMem IDA_mem;
646 realtype atolmin;
647
648 if (ida_mem==NULL) {
649 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDASVtolerances", MSG_NO_MEM);
650 return(IDA_MEM_NULL);
651 }
652 IDA_mem = (IDAMem) ida_mem;
653
654 if (IDA_mem->ida_MallocDone == SUNFALSE) {
655 IDAProcessError(IDA_mem, IDA_NO_MALLOC, "IDA", "IDASVtolerances", MSG_NO_MALLOC);
656 return(IDA_NO_MALLOC);
657 }
658
659 /* Check inputs */
660
661 if (reltol < ZERO) {
662 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASVtolerances", MSG_BAD_RTOL);
663 return(IDA_ILL_INPUT);
664 }
665
666 atolmin = N_VMin(abstol);
667 if (atolmin < ZERO) {
668 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASVtolerances", MSG_BAD_ATOL);
669 return(IDA_ILL_INPUT);
670 }
671
672 /* Copy tolerances into memory */
673
674 if ( !(IDA_mem->ida_VatolMallocDone) ) {
675 IDA_mem->ida_Vatol = N_VClone(IDA_mem->ida_ewt);
676 IDA_mem->ida_lrw += IDA_mem->ida_lrw1;
677 IDA_mem->ida_liw += IDA_mem->ida_liw1;
678 IDA_mem->ida_VatolMallocDone = SUNTRUE;
679 }
680
681 IDA_mem->ida_rtol = reltol;
682 N_VScale(ONE, abstol, IDA_mem->ida_Vatol);
683 IDA_mem->ida_atolmin0 = (atolmin == ZERO);
684
685 IDA_mem->ida_itol = IDA_SV;
686
687 IDA_mem->ida_user_efun = SUNFALSE;
688 IDA_mem->ida_efun = IDAEwtSet;
689 IDA_mem->ida_edata = NULL; /* will be set to ida_mem in InitialSetup */
690
691 return(IDA_SUCCESS);
692 }
693
694
IDAWFtolerances(void * ida_mem,IDAEwtFn efun)695 int IDAWFtolerances(void *ida_mem, IDAEwtFn efun)
696 {
697 IDAMem IDA_mem;
698
699 if (ida_mem==NULL) {
700 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDAWFtolerances", MSG_NO_MEM);
701 return(IDA_MEM_NULL);
702 }
703 IDA_mem = (IDAMem) ida_mem;
704
705 if (IDA_mem->ida_MallocDone == SUNFALSE) {
706 IDAProcessError(IDA_mem, IDA_NO_MALLOC, "IDA", "IDAWFtolerances", MSG_NO_MALLOC);
707 return(IDA_NO_MALLOC);
708 }
709
710 IDA_mem->ida_itol = IDA_WF;
711
712 IDA_mem->ida_user_efun = SUNTRUE;
713 IDA_mem->ida_efun = efun;
714 IDA_mem->ida_edata = NULL; /* will be set to user_data in InitialSetup */
715
716 return(IDA_SUCCESS);
717 }
718
719 /*-----------------------------------------------------------------*/
720
721 /*
722 * IDARootInit
723 *
724 * IDARootInit initializes a rootfinding problem to be solved
725 * during the integration of the DAE system. It loads the root
726 * function pointer and the number of root functions, and allocates
727 * workspace memory. The return value is IDA_SUCCESS = 0 if no
728 * errors occurred, or a negative value otherwise.
729 */
730
IDARootInit(void * ida_mem,int nrtfn,IDARootFn g)731 int IDARootInit(void *ida_mem, int nrtfn, IDARootFn g)
732 {
733 IDAMem IDA_mem;
734 int i, nrt;
735
736 /* Check ida_mem pointer */
737 if (ida_mem == NULL) {
738 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDARootInit", MSG_NO_MEM);
739 return(IDA_MEM_NULL);
740 }
741 IDA_mem = (IDAMem) ida_mem;
742
743 nrt = (nrtfn < 0) ? 0 : nrtfn;
744
745 /* If rerunning IDARootInit() with a different number of root
746 functions (changing number of gfun components), then free
747 currently held memory resources */
748 if ((nrt != IDA_mem->ida_nrtfn) && (IDA_mem->ida_nrtfn > 0)) {
749
750 free(IDA_mem->ida_glo); IDA_mem->ida_glo = NULL;
751 free(IDA_mem->ida_ghi); IDA_mem->ida_ghi = NULL;
752 free(IDA_mem->ida_grout); IDA_mem->ida_grout = NULL;
753 free(IDA_mem->ida_iroots); IDA_mem->ida_iroots = NULL;
754 free(IDA_mem->ida_rootdir); IDA_mem->ida_rootdir = NULL;
755 free(IDA_mem->ida_gactive); IDA_mem->ida_gactive = NULL;
756
757 IDA_mem->ida_lrw -= 3 * (IDA_mem->ida_nrtfn);
758 IDA_mem->ida_liw -= 3 * (IDA_mem->ida_nrtfn);
759
760 }
761
762 /* If IDARootInit() was called with nrtfn == 0, then set ida_nrtfn to
763 zero and ida_gfun to NULL before returning */
764 if (nrt == 0) {
765 IDA_mem->ida_nrtfn = nrt;
766 IDA_mem->ida_gfun = NULL;
767 return(IDA_SUCCESS);
768 }
769
770 /* If rerunning IDARootInit() with the same number of root functions
771 (not changing number of gfun components), then check if the root
772 function argument has changed */
773 /* If g != NULL then return as currently reserved memory resources
774 will suffice */
775 if (nrt == IDA_mem->ida_nrtfn) {
776 if (g != IDA_mem->ida_gfun) {
777 if (g == NULL) {
778 free(IDA_mem->ida_glo); IDA_mem->ida_glo = NULL;
779 free(IDA_mem->ida_ghi); IDA_mem->ida_ghi = NULL;
780 free(IDA_mem->ida_grout); IDA_mem->ida_grout = NULL;
781 free(IDA_mem->ida_iroots); IDA_mem->ida_iroots = NULL;
782 free(IDA_mem->ida_rootdir); IDA_mem->ida_rootdir = NULL;
783 free(IDA_mem->ida_gactive); IDA_mem->ida_gactive = NULL;
784
785 IDA_mem->ida_lrw -= 3*nrt;
786 IDA_mem->ida_liw -= 3*nrt;
787
788 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDARootInit", MSG_ROOT_FUNC_NULL);
789 return(IDA_ILL_INPUT);
790 }
791 else {
792 IDA_mem->ida_gfun = g;
793 return(IDA_SUCCESS);
794 }
795 }
796 else return(IDA_SUCCESS);
797 }
798
799 /* Set variable values in IDA memory block */
800 IDA_mem->ida_nrtfn = nrt;
801 if (g == NULL) {
802 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDARootInit", MSG_ROOT_FUNC_NULL);
803 return(IDA_ILL_INPUT);
804 }
805 else IDA_mem->ida_gfun = g;
806
807 /* Allocate necessary memory and return */
808 IDA_mem->ida_glo = NULL;
809 IDA_mem->ida_glo = (realtype *) malloc(nrt*sizeof(realtype));
810 if (IDA_mem->ida_glo == NULL) {
811 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDARootInit", MSG_MEM_FAIL);
812 return(IDA_MEM_FAIL);
813 }
814
815 IDA_mem->ida_ghi = NULL;
816 IDA_mem->ida_ghi = (realtype *) malloc(nrt*sizeof(realtype));
817 if (IDA_mem->ida_ghi == NULL) {
818 free(IDA_mem->ida_glo); IDA_mem->ida_glo = NULL;
819 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDARootInit", MSG_MEM_FAIL);
820 return(IDA_MEM_FAIL);
821 }
822
823 IDA_mem->ida_grout = NULL;
824 IDA_mem->ida_grout = (realtype *) malloc(nrt*sizeof(realtype));
825 if (IDA_mem->ida_grout == NULL) {
826 free(IDA_mem->ida_glo); IDA_mem->ida_glo = NULL;
827 free(IDA_mem->ida_ghi); IDA_mem->ida_ghi = NULL;
828 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDARootInit", MSG_MEM_FAIL);
829 return(IDA_MEM_FAIL);
830 }
831
832 IDA_mem->ida_iroots = NULL;
833 IDA_mem->ida_iroots = (int *) malloc(nrt*sizeof(int));
834 if (IDA_mem->ida_iroots == NULL) {
835 free(IDA_mem->ida_glo); IDA_mem->ida_glo = NULL;
836 free(IDA_mem->ida_ghi); IDA_mem->ida_ghi = NULL;
837 free(IDA_mem->ida_grout); IDA_mem->ida_grout = NULL;
838 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDARootInit", MSG_MEM_FAIL);
839 return(IDA_MEM_FAIL);
840 }
841
842 IDA_mem->ida_rootdir = NULL;
843 IDA_mem->ida_rootdir = (int *) malloc(nrt*sizeof(int));
844 if (IDA_mem->ida_rootdir == NULL) {
845 free(IDA_mem->ida_glo); IDA_mem->ida_glo = NULL;
846 free(IDA_mem->ida_ghi); IDA_mem->ida_ghi = NULL;
847 free(IDA_mem->ida_grout); IDA_mem->ida_grout = NULL;
848 free(IDA_mem->ida_iroots); IDA_mem->ida_iroots = NULL;
849 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDARootInit", MSG_MEM_FAIL);
850 return(IDA_MEM_FAIL);
851 }
852
853 IDA_mem->ida_gactive = NULL;
854 IDA_mem->ida_gactive = (booleantype *) malloc(nrt*sizeof(booleantype));
855 if (IDA_mem->ida_gactive == NULL) {
856 free(IDA_mem->ida_glo); IDA_mem->ida_glo = NULL;
857 free(IDA_mem->ida_ghi); IDA_mem->ida_ghi = NULL;
858 free(IDA_mem->ida_grout); IDA_mem->ida_grout = NULL;
859 free(IDA_mem->ida_iroots); IDA_mem->ida_iroots = NULL;
860 free(IDA_mem->ida_rootdir); IDA_mem->ida_rootdir = NULL;
861 IDAProcessError(IDA_mem, IDA_MEM_FAIL, "IDA", "IDARootInit", MSG_MEM_FAIL);
862 return(IDA_MEM_FAIL);
863 }
864
865 /* Set default values for rootdir (both directions) */
866 for(i=0; i<nrt; i++) IDA_mem->ida_rootdir[i] = 0;
867
868 /* Set default values for gactive (all active) */
869 for(i=0; i<nrt; i++) IDA_mem->ida_gactive[i] = SUNTRUE;
870
871 IDA_mem->ida_lrw += 3*nrt;
872 IDA_mem->ida_liw += 3*nrt;
873
874 return(IDA_SUCCESS);
875 }
876
877
878 /*
879 * -----------------------------------------------------------------
880 * Main solver function
881 * -----------------------------------------------------------------
882 */
883
884 /*
885 * IDASolve
886 *
887 * This routine is the main driver of the IDA package.
888 *
889 * It integrates over an independent variable interval defined by the user,
890 * by calling IDAStep to take internal independent variable steps.
891 *
892 * The first time that IDASolve is called for a successfully initialized
893 * problem, it computes a tentative initial step size.
894 *
895 * IDASolve supports two modes, specified by itask:
896 * In the IDA_NORMAL mode, the solver steps until it passes tout and then
897 * interpolates to obtain y(tout) and yp(tout).
898 * In the IDA_ONE_STEP mode, it takes one internal step and returns.
899 *
900 * IDASolve returns integer values corresponding to success and failure as below:
901 *
902 * successful returns:
903 *
904 * IDA_SUCCESS
905 * IDA_TSTOP_RETURN
906 *
907 * failed returns:
908 *
909 * IDA_ILL_INPUT
910 * IDA_TOO_MUCH_WORK
911 * IDA_MEM_NULL
912 * IDA_TOO_MUCH_ACC
913 * IDA_CONV_FAIL
914 * IDA_LSETUP_FAIL
915 * IDA_LSOLVE_FAIL
916 * IDA_CONSTR_FAIL
917 * IDA_ERR_FAIL
918 * IDA_REP_RES_ERR
919 * IDA_RES_FAIL
920 */
921
IDASolve(void * ida_mem,realtype tout,realtype * tret,N_Vector yret,N_Vector ypret,int itask)922 int IDASolve(void *ida_mem, realtype tout, realtype *tret,
923 N_Vector yret, N_Vector ypret, int itask)
924 {
925 long int nstloc;
926 int sflag, istate, ier, irfndp, ir;
927 realtype tdist, troundoff, ypnorm, rh, nrm;
928 IDAMem IDA_mem;
929 booleantype inactive_roots;
930
931 /* Check for legal inputs in all cases. */
932
933 if (ida_mem == NULL) {
934 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDASolve", MSG_NO_MEM);
935 return(IDA_MEM_NULL);
936 }
937 IDA_mem = (IDAMem) ida_mem;
938
939 /* Check if problem was malloc'ed */
940
941 if (IDA_mem->ida_MallocDone == SUNFALSE) {
942 IDAProcessError(IDA_mem, IDA_NO_MALLOC, "IDA", "IDASolve", MSG_NO_MALLOC);
943 return(IDA_NO_MALLOC);
944 }
945
946 /* Check for legal arguments */
947
948 if (yret == NULL) {
949 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve", MSG_YRET_NULL);
950 return(IDA_ILL_INPUT);
951 }
952 IDA_mem->ida_yy = yret;
953
954 if (ypret == NULL) {
955 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve", MSG_YPRET_NULL);
956 return(IDA_ILL_INPUT);
957 }
958 IDA_mem->ida_yp = ypret;
959
960 if (tret == NULL) {
961 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve", MSG_TRET_NULL);
962 return(IDA_ILL_INPUT);
963 }
964
965 if ((itask != IDA_NORMAL) && (itask != IDA_ONE_STEP)) {
966 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve", MSG_BAD_ITASK);
967 return(IDA_ILL_INPUT);
968 }
969
970 if (itask == IDA_NORMAL) IDA_mem->ida_toutc = tout;
971 IDA_mem->ida_taskc = itask;
972
973 if (IDA_mem->ida_nst == 0) { /* This is the first call */
974
975 /* Check inputs to IDA for correctness and consistency */
976
977 if (IDA_mem->ida_SetupDone == SUNFALSE) {
978 ier = IDAInitialSetup(IDA_mem);
979 if (ier != IDA_SUCCESS) return(ier);
980 IDA_mem->ida_SetupDone = SUNTRUE;
981 }
982
983 /* On first call, check for tout - tn too small, set initial hh,
984 check for approach to tstop, and scale phi[1] by hh.
985 Also check for zeros of root function g at and near t0. */
986
987 tdist = SUNRabs(tout - IDA_mem->ida_tn);
988 if (tdist == ZERO) {
989 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve", MSG_TOO_CLOSE);
990 return(IDA_ILL_INPUT);
991 }
992 troundoff = TWO * IDA_mem->ida_uround * (SUNRabs(IDA_mem->ida_tn) + SUNRabs(tout));
993 if (tdist < troundoff) {
994 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve", MSG_TOO_CLOSE);
995 return(IDA_ILL_INPUT);
996 }
997
998 IDA_mem->ida_hh = IDA_mem->ida_hin;
999 if ( (IDA_mem->ida_hh != ZERO) && ((tout-IDA_mem->ida_tn)*IDA_mem->ida_hh < ZERO) ) {
1000 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve", MSG_BAD_HINIT);
1001 return(IDA_ILL_INPUT);
1002 }
1003
1004 if (IDA_mem->ida_hh == ZERO) {
1005 IDA_mem->ida_hh = PT001*tdist;
1006 ypnorm = IDAWrmsNorm(IDA_mem, IDA_mem->ida_phi[1],
1007 IDA_mem->ida_ewt, IDA_mem->ida_suppressalg);
1008 if (ypnorm > HALF / IDA_mem->ida_hh) IDA_mem->ida_hh = HALF/ypnorm;
1009 if (tout < IDA_mem->ida_tn) IDA_mem->ida_hh = -IDA_mem->ida_hh;
1010 }
1011
1012 rh = SUNRabs(IDA_mem->ida_hh) * IDA_mem->ida_hmax_inv;
1013 if (rh > ONE) IDA_mem->ida_hh /= rh;
1014
1015 if (IDA_mem->ida_tstopset) {
1016 if ( (IDA_mem->ida_tstop - IDA_mem->ida_tn)*IDA_mem->ida_hh <= ZERO) {
1017 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve",
1018 MSG_BAD_TSTOP, IDA_mem->ida_tstop, IDA_mem->ida_tn);
1019 return(IDA_ILL_INPUT);
1020 }
1021 if ( (IDA_mem->ida_tn + IDA_mem->ida_hh - IDA_mem->ida_tstop)*IDA_mem->ida_hh > ZERO)
1022 IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn)*(ONE - FOUR * IDA_mem->ida_uround);
1023 }
1024
1025 IDA_mem->ida_h0u = IDA_mem->ida_hh;
1026 IDA_mem->ida_kk = 0;
1027 IDA_mem->ida_kused = 0; /* set in case of an error return before a step */
1028
1029 /* Check for exact zeros of the root functions at or near t0. */
1030 if (IDA_mem->ida_nrtfn > 0) {
1031 ier = IDARcheck1(IDA_mem);
1032 if (ier == IDA_RTFUNC_FAIL) {
1033 IDAProcessError(IDA_mem, IDA_RTFUNC_FAIL, "IDA", "IDARcheck1",
1034 MSG_RTFUNC_FAILED, IDA_mem->ida_tn);
1035 return(IDA_RTFUNC_FAIL);
1036 }
1037 }
1038
1039 /* set phi[1] = hh*y' */
1040 N_VScale(IDA_mem->ida_hh, IDA_mem->ida_phi[1], IDA_mem->ida_phi[1]);
1041
1042 /* Set the convergence test constants epsNewt and toldel */
1043 IDA_mem->ida_epsNewt = IDA_mem->ida_epcon;
1044 IDA_mem->ida_toldel = PT0001 * IDA_mem->ida_epsNewt;
1045
1046 } /* end of first-call block. */
1047
1048 /* Call lperf function and set nstloc for later performance testing. */
1049
1050 if (IDA_mem->ida_lperf != NULL)
1051 IDA_mem->ida_lperf(IDA_mem, 0);
1052 nstloc = 0;
1053
1054 /* If not the first call, perform all stopping tests. */
1055
1056 if (IDA_mem->ida_nst > 0) {
1057
1058 /* First, check for a root in the last step taken, other than the
1059 last root found, if any. If itask = IDA_ONE_STEP and y(tn) was not
1060 returned because of an intervening root, return y(tn) now. */
1061
1062 if (IDA_mem->ida_nrtfn > 0) {
1063
1064 irfndp = IDA_mem->ida_irfnd;
1065
1066 ier = IDARcheck2(IDA_mem);
1067
1068 if (ier == CLOSERT) {
1069 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDARcheck2",
1070 MSG_CLOSE_ROOTS, IDA_mem->ida_tlo);
1071 return(IDA_ILL_INPUT);
1072 } else if (ier == IDA_RTFUNC_FAIL) {
1073 IDAProcessError(IDA_mem, IDA_RTFUNC_FAIL, "IDA", "IDARcheck2",
1074 MSG_RTFUNC_FAILED, IDA_mem->ida_tlo);
1075 return(IDA_RTFUNC_FAIL);
1076 } else if (ier == RTFOUND) {
1077 IDA_mem->ida_tretlast = *tret = IDA_mem->ida_tlo;
1078 return(IDA_ROOT_RETURN);
1079 }
1080
1081 /* If tn is distinct from tretlast (within roundoff),
1082 check remaining interval for roots */
1083 troundoff = HUNDRED * IDA_mem->ida_uround * (SUNRabs(IDA_mem->ida_tn) + SUNRabs(IDA_mem->ida_hh));
1084 if ( SUNRabs(IDA_mem->ida_tn - IDA_mem->ida_tretlast) > troundoff ) {
1085 ier = IDARcheck3(IDA_mem);
1086 if (ier == IDA_SUCCESS) { /* no root found */
1087 IDA_mem->ida_irfnd = 0;
1088 if ((irfndp == 1) && (itask == IDA_ONE_STEP)) {
1089 IDA_mem->ida_tretlast = *tret = IDA_mem->ida_tn;
1090 ier = IDAGetSolution(IDA_mem, IDA_mem->ida_tn, yret, ypret);
1091 return(IDA_SUCCESS);
1092 }
1093 } else if (ier == RTFOUND) { /* a new root was found */
1094 IDA_mem->ida_irfnd = 1;
1095 IDA_mem->ida_tretlast = *tret = IDA_mem->ida_tlo;
1096 return(IDA_ROOT_RETURN);
1097 } else if (ier == IDA_RTFUNC_FAIL) { /* g failed */
1098 IDAProcessError(IDA_mem, IDA_RTFUNC_FAIL, "IDA", "IDARcheck3",
1099 MSG_RTFUNC_FAILED, IDA_mem->ida_tlo);
1100 return(IDA_RTFUNC_FAIL);
1101 }
1102 }
1103
1104 } /* end of root stop check */
1105
1106
1107 /* Now test for all other stop conditions. */
1108
1109 istate = IDAStopTest1(IDA_mem, tout, tret, yret, ypret, itask);
1110 if (istate != CONTINUE_STEPS) return(istate);
1111 }
1112
1113 /* Looping point for internal steps. */
1114
1115 for(;;) {
1116
1117 /* Check for too many steps taken. */
1118
1119 if ( (IDA_mem->ida_mxstep>0) && (nstloc >= IDA_mem->ida_mxstep) ) {
1120 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve",
1121 MSG_MAX_STEPS, IDA_mem->ida_tn);
1122 istate = IDA_TOO_MUCH_WORK;
1123 *tret = IDA_mem->ida_tretlast = IDA_mem->ida_tn;
1124 break; /* Here yy=yret and yp=ypret already have the current solution. */
1125 }
1126
1127 /* Call lperf to generate warnings of poor performance. */
1128
1129 if (IDA_mem->ida_lperf != NULL)
1130 IDA_mem->ida_lperf(IDA_mem, 1);
1131
1132 /* Reset and check ewt (if not first call). */
1133
1134 if (IDA_mem->ida_nst > 0) {
1135
1136 ier = IDA_mem->ida_efun(IDA_mem->ida_phi[0], IDA_mem->ida_ewt,
1137 IDA_mem->ida_edata);
1138
1139 if (ier != 0) {
1140
1141 if (IDA_mem->ida_itol == IDA_WF)
1142 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve",
1143 MSG_EWT_NOW_FAIL, IDA_mem->ida_tn);
1144 else
1145 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve",
1146 MSG_EWT_NOW_BAD, IDA_mem->ida_tn);
1147
1148 istate = IDA_ILL_INPUT;
1149 ier = IDAGetSolution(IDA_mem, IDA_mem->ida_tn, yret, ypret);
1150 *tret = IDA_mem->ida_tretlast = IDA_mem->ida_tn;
1151 break;
1152
1153 }
1154
1155 }
1156
1157 /* Check for too much accuracy requested. */
1158
1159 nrm = IDAWrmsNorm(IDA_mem, IDA_mem->ida_phi[0], IDA_mem->ida_ewt,
1160 IDA_mem->ida_suppressalg);
1161 IDA_mem->ida_tolsf = IDA_mem->ida_uround * nrm;
1162 if (IDA_mem->ida_tolsf > ONE) {
1163 IDA_mem->ida_tolsf *= TEN;
1164 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve",
1165 MSG_TOO_MUCH_ACC, IDA_mem->ida_tn);
1166 istate = IDA_TOO_MUCH_ACC;
1167 *tret = IDA_mem->ida_tretlast = IDA_mem->ida_tn;
1168 if (IDA_mem->ida_nst > 0) ier = IDAGetSolution(IDA_mem, IDA_mem->ida_tn, yret, ypret);
1169 break;
1170 }
1171
1172 /* Call IDAStep to take a step. */
1173
1174 sflag = IDAStep(IDA_mem);
1175
1176 /* Process all failed-step cases, and exit loop. */
1177
1178 if (sflag != IDA_SUCCESS) {
1179 istate = IDAHandleFailure(IDA_mem, sflag);
1180 *tret = IDA_mem->ida_tretlast = IDA_mem->ida_tn;
1181 ier = IDAGetSolution(IDA_mem, IDA_mem->ida_tn, yret, ypret);
1182 break;
1183 }
1184
1185 nstloc++;
1186
1187 /* After successful step, check for stop conditions; continue or break. */
1188
1189 /* First check for root in the last step taken. */
1190
1191 if (IDA_mem->ida_nrtfn > 0) {
1192
1193 ier = IDARcheck3(IDA_mem);
1194
1195 if (ier == RTFOUND) { /* A new root was found */
1196 IDA_mem->ida_irfnd = 1;
1197 istate = IDA_ROOT_RETURN;
1198 IDA_mem->ida_tretlast = *tret = IDA_mem->ida_tlo;
1199 break;
1200 } else if (ier == IDA_RTFUNC_FAIL) { /* g failed */
1201 IDAProcessError(IDA_mem, IDA_RTFUNC_FAIL, "IDA", "IDARcheck3",
1202 MSG_RTFUNC_FAILED, IDA_mem->ida_tlo);
1203 istate = IDA_RTFUNC_FAIL;
1204 break;
1205 }
1206
1207 /* If we are at the end of the first step and we still have
1208 * some event functions that are inactive, issue a warning
1209 * as this may indicate a user error in the implementation
1210 * of the root function. */
1211
1212 if (IDA_mem->ida_nst==1) {
1213 inactive_roots = SUNFALSE;
1214 for (ir=0; ir<IDA_mem->ida_nrtfn; ir++) {
1215 if (!IDA_mem->ida_gactive[ir]) {
1216 inactive_roots = SUNTRUE;
1217 break;
1218 }
1219 }
1220 if ((IDA_mem->ida_mxgnull > 0) && inactive_roots) {
1221 IDAProcessError(IDA_mem, IDA_WARNING, "IDA", "IDASolve",
1222 MSG_INACTIVE_ROOTS);
1223 }
1224 }
1225
1226 }
1227
1228 /* Now check all other stop conditions. */
1229
1230 istate = IDAStopTest2(IDA_mem, tout, tret, yret, ypret, itask);
1231 if (istate != CONTINUE_STEPS) break;
1232
1233 } /* End of step loop */
1234
1235 return(istate);
1236 }
1237
1238 /*
1239 * -----------------------------------------------------------------
1240 * Interpolated output and extraction functions
1241 * -----------------------------------------------------------------
1242 */
1243
1244 /*
1245 * IDAGetDky
1246 *
1247 * This routine evaluates the k-th derivative of y(t) as the value of
1248 * the k-th derivative of the interpolating polynomial at the independent
1249 * variable t, and stores the results in the vector dky. It uses the current
1250 * independent variable value, tn, and the method order last used, kused.
1251 *
1252 * The return values are:
1253 * IDA_SUCCESS if t is legal
1254 * IDA_BAD_T if t is not within the interval of the last step taken
1255 * IDA_BAD_DKY if the dky vector is NULL
1256 * IDA_BAD_K if the requested k is not in the range [0,order used]
1257 * IDA_VECTOROP_ERR if the fused vector operation fails
1258 *
1259 */
1260
IDAGetDky(void * ida_mem,realtype t,int k,N_Vector dky)1261 int IDAGetDky(void *ida_mem, realtype t, int k, N_Vector dky)
1262 {
1263 IDAMem IDA_mem;
1264 realtype tfuzz, tp, delt, psij_1;
1265 int i, j, retval;
1266 realtype cjk [MXORDP1];
1267 realtype cjk_1[MXORDP1];
1268
1269 /* Check ida_mem */
1270 if (ida_mem == NULL) {
1271 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDAGetDky", MSG_NO_MEM);
1272 return (IDA_MEM_NULL);
1273 }
1274 IDA_mem = (IDAMem) ida_mem;
1275
1276 if (dky == NULL) {
1277 IDAProcessError(IDA_mem, IDA_BAD_DKY, "IDA", "IDAGetDky", MSG_NULL_DKY);
1278 return(IDA_BAD_DKY);
1279 }
1280
1281 if ((k < 0) || (k > IDA_mem->ida_kused)) {
1282 IDAProcessError(IDA_mem, IDA_BAD_K, "IDA", "IDAGetDky", MSG_BAD_K);
1283 return(IDA_BAD_K);
1284 }
1285
1286 /* Check t for legality. Here tn - hused is t_{n-1}. */
1287
1288 tfuzz = HUNDRED * IDA_mem->ida_uround * (SUNRabs(IDA_mem->ida_tn) + SUNRabs(IDA_mem->ida_hh));
1289 if (IDA_mem->ida_hh < ZERO) tfuzz = - tfuzz;
1290 tp = IDA_mem->ida_tn - IDA_mem->ida_hused - tfuzz;
1291 if ((t - tp)*IDA_mem->ida_hh < ZERO) {
1292 IDAProcessError(IDA_mem, IDA_BAD_T, "IDA", "IDAGetDky", MSG_BAD_T, t,
1293 IDA_mem->ida_tn-IDA_mem->ida_hused, IDA_mem->ida_tn);
1294 return(IDA_BAD_T);
1295 }
1296
1297 /* Initialize the c_j^(k) and c_k^(k-1) */
1298 for(i=0; i<MXORDP1; i++) {
1299 cjk [i] = 0;
1300 cjk_1[i] = 0;
1301 }
1302
1303 delt = t-IDA_mem->ida_tn;
1304
1305 for(i=0; i<=k; i++) {
1306
1307 /* The below reccurence is used to compute the k-th derivative of the solution:
1308 c_j^(k) = ( k * c_{j-1}^(k-1) + c_{j-1}^{k} (Delta+psi_{j-1}) ) / psi_j
1309
1310 Translated in indexes notation:
1311 cjk[j] = ( k*cjk_1[j-1] + cjk[j-1]*(delt+psi[j-2]) ) / psi[j-1]
1312
1313 For k=0, j=1: c_1 = c_0^(-1) + (delt+psi[-1]) / psi[0]
1314
1315 In order to be able to deal with k=0 in the same way as for k>0, the
1316 following conventions were adopted:
1317 - c_0(t) = 1 , c_0^(-1)(t)=0
1318 - psij_1 stands for psi[-1]=0 when j=1
1319 for psi[j-2] when j>1
1320 */
1321 if(i==0) {
1322
1323 cjk[i] = 1;
1324 psij_1 = 0;
1325 }else {
1326 /* i i-1 1
1327 c_i^(i) can be always updated since c_i^(i) = ----- -------- ... -----
1328 psi_j psi_{j-1} psi_1
1329 */
1330 cjk[i] = cjk[i-1]*i / IDA_mem->ida_psi[i-1];
1331 psij_1 = IDA_mem->ida_psi[i-1];
1332 }
1333
1334 /* update c_j^(i) */
1335
1336 /*j does not need to go till kused */
1337 for(j=i+1; j<=IDA_mem->ida_kused-k+i; j++) {
1338
1339 cjk[j] = ( i* cjk_1[j-1] + cjk[j-1] * (delt + psij_1) ) / IDA_mem->ida_psi[j-1];
1340 psij_1 = IDA_mem->ida_psi[j-1];
1341 }
1342
1343 /* save existing c_j^(i)'s */
1344 for(j=i+1; j<=IDA_mem->ida_kused-k+i; j++) cjk_1[j] = cjk[j];
1345 }
1346
1347 /* Compute sum (c_j(t) * phi(t)) */
1348
1349 /* Sum j=k to j<=IDA_mem->ida_kused */
1350 retval = N_VLinearCombination(IDA_mem->ida_kused-k+1, cjk+k,
1351 IDA_mem->ida_phi+k, dky);
1352 if (retval != IDA_SUCCESS) return(IDA_VECTOROP_ERR);
1353
1354 return(IDA_SUCCESS);
1355 }
1356
1357 /*
1358 * IDAComputeY
1359 *
1360 * Computes y based on the current prediction and given correction.
1361 */
IDAComputeY(void * ida_mem,N_Vector ycor,N_Vector y)1362 int IDAComputeY(void *ida_mem, N_Vector ycor, N_Vector y)
1363 {
1364 IDAMem IDA_mem;
1365
1366 if (ida_mem==NULL) {
1367 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDAComputeY", MSG_NO_MEM);
1368 return(IDA_MEM_NULL);
1369 }
1370
1371 IDA_mem = (IDAMem) ida_mem;
1372
1373 N_VLinearSum(ONE, IDA_mem->ida_yypredict, ONE, ycor, y);
1374
1375 return(IDA_SUCCESS);
1376 }
1377
1378 /*
1379 * IDAComputeYp
1380 *
1381 * Computes y' based on the current prediction and given correction.
1382 */
IDAComputeYp(void * ida_mem,N_Vector ycor,N_Vector yp)1383 int IDAComputeYp(void *ida_mem, N_Vector ycor, N_Vector yp)
1384 {
1385 IDAMem IDA_mem;
1386
1387 if (ida_mem==NULL) {
1388 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDAComputeYp", MSG_NO_MEM);
1389 return(IDA_MEM_NULL);
1390 }
1391
1392 IDA_mem = (IDAMem) ida_mem;
1393
1394 N_VLinearSum(ONE, IDA_mem->ida_yppredict, IDA_mem->ida_cj, ycor, yp);
1395
1396 return(IDA_SUCCESS);
1397 }
1398
1399 /*
1400 * -----------------------------------------------------------------
1401 * Deallocation function
1402 * -----------------------------------------------------------------
1403 */
1404
1405 /*
1406 * IDAFree
1407 *
1408 * This routine frees the problem memory allocated by IDAInit
1409 * Such memory includes all the vectors allocated by IDAAllocVectors,
1410 * and the memory lmem for the linear solver (deallocated by a call
1411 * to lfree).
1412 */
1413
IDAFree(void ** ida_mem)1414 void IDAFree(void **ida_mem)
1415 {
1416 IDAMem IDA_mem;
1417
1418 if (*ida_mem == NULL) return;
1419
1420 IDA_mem = (IDAMem) (*ida_mem);
1421
1422 IDAFreeVectors(IDA_mem);
1423
1424 /* if IDA created the NLS object then free it */
1425 if (IDA_mem->ownNLS) {
1426 SUNNonlinSolFree(IDA_mem->NLS);
1427 IDA_mem->ownNLS = SUNFALSE;
1428 IDA_mem->NLS = NULL;
1429 }
1430
1431 if (IDA_mem->ida_lfree != NULL)
1432 IDA_mem->ida_lfree(IDA_mem);
1433
1434 if (IDA_mem->ida_nrtfn > 0) {
1435 free(IDA_mem->ida_glo); IDA_mem->ida_glo = NULL;
1436 free(IDA_mem->ida_ghi); IDA_mem->ida_ghi = NULL;
1437 free(IDA_mem->ida_grout); IDA_mem->ida_grout = NULL;
1438 free(IDA_mem->ida_iroots); IDA_mem->ida_iroots = NULL;
1439 free(IDA_mem->ida_rootdir); IDA_mem->ida_rootdir = NULL;
1440 free(IDA_mem->ida_gactive); IDA_mem->ida_gactive = NULL;
1441 }
1442
1443 free(*ida_mem);
1444 *ida_mem = NULL;
1445 }
1446
1447 /*
1448 * =================================================================
1449 * PRIVATE FUNCTIONS
1450 * =================================================================
1451 */
1452
1453 /*
1454 * IDACheckNvector
1455 *
1456 * This routine checks if all required vector operations are present.
1457 * If any of them is missing it returns SUNFALSE.
1458 */
1459
IDACheckNvector(N_Vector tmpl)1460 static booleantype IDACheckNvector(N_Vector tmpl)
1461 {
1462 if ((tmpl->ops->nvclone == NULL) ||
1463 (tmpl->ops->nvdestroy == NULL) ||
1464 (tmpl->ops->nvlinearsum == NULL) ||
1465 (tmpl->ops->nvconst == NULL) ||
1466 (tmpl->ops->nvprod == NULL) ||
1467 (tmpl->ops->nvscale == NULL) ||
1468 (tmpl->ops->nvabs == NULL) ||
1469 (tmpl->ops->nvinv == NULL) ||
1470 (tmpl->ops->nvaddconst == NULL) ||
1471 (tmpl->ops->nvwrmsnorm == NULL) ||
1472 (tmpl->ops->nvmin == NULL))
1473 return(SUNFALSE);
1474 else
1475 return(SUNTRUE);
1476 }
1477
1478 /*
1479 * -----------------------------------------------------------------
1480 * Memory allocation/deallocation
1481 * -----------------------------------------------------------------
1482 */
1483
1484 /*
1485 * IDAAllocVectors
1486 *
1487 * This routine allocates the IDA vectors ewt, tempv1, tempv2, and
1488 * phi[0], ..., phi[maxord].
1489 * If all memory allocations are successful, IDAAllocVectors returns
1490 * SUNTRUE. Otherwise all allocated memory is freed and IDAAllocVectors
1491 * returns SUNFALSE.
1492 * This routine also sets the optional outputs lrw and liw, which are
1493 * (respectively) the lengths of the real and integer work spaces
1494 * allocated here.
1495 */
1496
IDAAllocVectors(IDAMem IDA_mem,N_Vector tmpl)1497 static booleantype IDAAllocVectors(IDAMem IDA_mem, N_Vector tmpl)
1498 {
1499 int i, j, maxcol;
1500
1501 /* Allocate ewt, ee, delta, yypredict, yppredict, savres, tempv1, tempv2, tempv3 */
1502
1503 IDA_mem->ida_ewt = N_VClone(tmpl);
1504 if (IDA_mem->ida_ewt == NULL) return(SUNFALSE);
1505
1506 IDA_mem->ida_ee = N_VClone(tmpl);
1507 if (IDA_mem->ida_ee == NULL) {
1508 N_VDestroy(IDA_mem->ida_ewt);
1509 return(SUNFALSE);
1510 }
1511
1512 IDA_mem->ida_delta = N_VClone(tmpl);
1513 if (IDA_mem->ida_delta == NULL) {
1514 N_VDestroy(IDA_mem->ida_ewt);
1515 N_VDestroy(IDA_mem->ida_ee);
1516 return(SUNFALSE);
1517 }
1518
1519 IDA_mem->ida_yypredict = N_VClone(tmpl);
1520 if (IDA_mem->ida_yypredict == NULL) {
1521 N_VDestroy(IDA_mem->ida_ewt);
1522 N_VDestroy(IDA_mem->ida_ee);
1523 N_VDestroy(IDA_mem->ida_delta);
1524 return(SUNFALSE);
1525 }
1526
1527 IDA_mem->ida_yppredict = N_VClone(tmpl);
1528 if (IDA_mem->ida_yppredict == NULL) {
1529 N_VDestroy(IDA_mem->ida_ewt);
1530 N_VDestroy(IDA_mem->ida_ee);
1531 N_VDestroy(IDA_mem->ida_delta);
1532 N_VDestroy(IDA_mem->ida_yypredict);
1533 return(SUNFALSE);
1534 }
1535
1536 IDA_mem->ida_savres = N_VClone(tmpl);
1537 if (IDA_mem->ida_savres == NULL) {
1538 N_VDestroy(IDA_mem->ida_ewt);
1539 N_VDestroy(IDA_mem->ida_ee);
1540 N_VDestroy(IDA_mem->ida_delta);
1541 N_VDestroy(IDA_mem->ida_yypredict);
1542 N_VDestroy(IDA_mem->ida_yppredict);
1543 return(SUNFALSE);
1544 }
1545
1546 IDA_mem->ida_tempv1 = N_VClone(tmpl);
1547 if (IDA_mem->ida_tempv1 == NULL) {
1548 N_VDestroy(IDA_mem->ida_ewt);
1549 N_VDestroy(IDA_mem->ida_ee);
1550 N_VDestroy(IDA_mem->ida_delta);
1551 N_VDestroy(IDA_mem->ida_yypredict);
1552 N_VDestroy(IDA_mem->ida_yppredict);
1553 N_VDestroy(IDA_mem->ida_savres);
1554 return(SUNFALSE);
1555 }
1556
1557 IDA_mem->ida_tempv2 = N_VClone(tmpl);
1558 if (IDA_mem->ida_tempv2 == NULL) {
1559 N_VDestroy(IDA_mem->ida_ewt);
1560 N_VDestroy(IDA_mem->ida_ee);
1561 N_VDestroy(IDA_mem->ida_delta);
1562 N_VDestroy(IDA_mem->ida_yypredict);
1563 N_VDestroy(IDA_mem->ida_yppredict);
1564 N_VDestroy(IDA_mem->ida_savres);
1565 N_VDestroy(IDA_mem->ida_tempv1);
1566 return(SUNFALSE);
1567 }
1568
1569 IDA_mem->ida_tempv3 = N_VClone(tmpl);
1570 if (IDA_mem->ida_tempv3 == NULL) {
1571 N_VDestroy(IDA_mem->ida_ewt);
1572 N_VDestroy(IDA_mem->ida_ee);
1573 N_VDestroy(IDA_mem->ida_delta);
1574 N_VDestroy(IDA_mem->ida_yypredict);
1575 N_VDestroy(IDA_mem->ida_yppredict);
1576 N_VDestroy(IDA_mem->ida_savres);
1577 N_VDestroy(IDA_mem->ida_tempv1);
1578 N_VDestroy(IDA_mem->ida_tempv2);
1579 return(SUNFALSE);
1580 }
1581
1582 /* Allocate phi[0] ... phi[maxord]. Make sure phi[2] and phi[3] are
1583 allocated (for use as temporary vectors), regardless of maxord. */
1584
1585 maxcol = SUNMAX(IDA_mem->ida_maxord,3);
1586 for (j=0; j <= maxcol; j++) {
1587 IDA_mem->ida_phi[j] = N_VClone(tmpl);
1588 if (IDA_mem->ida_phi[j] == NULL) {
1589 N_VDestroy(IDA_mem->ida_ewt);
1590 N_VDestroy(IDA_mem->ida_ee);
1591 N_VDestroy(IDA_mem->ida_delta);
1592 N_VDestroy(IDA_mem->ida_yypredict);
1593 N_VDestroy(IDA_mem->ida_yppredict);
1594 N_VDestroy(IDA_mem->ida_savres);
1595 N_VDestroy(IDA_mem->ida_tempv1);
1596 N_VDestroy(IDA_mem->ida_tempv2);
1597 N_VDestroy(IDA_mem->ida_tempv3);
1598 for (i=0; i < j; i++) N_VDestroy(IDA_mem->ida_phi[i]);
1599 return(SUNFALSE);
1600 }
1601 }
1602
1603 /* Update solver workspace lengths */
1604 IDA_mem->ida_lrw += (maxcol + 10)*IDA_mem->ida_lrw1;
1605 IDA_mem->ida_liw += (maxcol + 10)*IDA_mem->ida_liw1;
1606
1607 /* Store the value of maxord used here */
1608 IDA_mem->ida_maxord_alloc = IDA_mem->ida_maxord;
1609
1610 return(SUNTRUE);
1611 }
1612
1613 /*
1614 * IDAfreeVectors
1615 *
1616 * This routine frees the IDA vectors allocated for IDA.
1617 */
1618
IDAFreeVectors(IDAMem IDA_mem)1619 static void IDAFreeVectors(IDAMem IDA_mem)
1620 {
1621 int j, maxcol;
1622
1623 N_VDestroy(IDA_mem->ida_ewt); IDA_mem->ida_ewt = NULL;
1624 N_VDestroy(IDA_mem->ida_ee); IDA_mem->ida_ee = NULL;
1625 N_VDestroy(IDA_mem->ida_delta); IDA_mem->ida_delta = NULL;
1626 N_VDestroy(IDA_mem->ida_yypredict); IDA_mem->ida_yypredict = NULL;
1627 N_VDestroy(IDA_mem->ida_yppredict); IDA_mem->ida_yppredict = NULL;
1628 N_VDestroy(IDA_mem->ida_savres); IDA_mem->ida_savres = NULL;
1629 N_VDestroy(IDA_mem->ida_tempv1); IDA_mem->ida_tempv1 = NULL;
1630 N_VDestroy(IDA_mem->ida_tempv2); IDA_mem->ida_tempv2 = NULL;
1631 N_VDestroy(IDA_mem->ida_tempv3); IDA_mem->ida_tempv3 = NULL;
1632 maxcol = SUNMAX(IDA_mem->ida_maxord_alloc,3);
1633 for(j=0; j <= maxcol; j++) {
1634 N_VDestroy(IDA_mem->ida_phi[j]);
1635 IDA_mem->ida_phi[j] = NULL;
1636 }
1637
1638 IDA_mem->ida_lrw -= (maxcol + 10)*IDA_mem->ida_lrw1;
1639 IDA_mem->ida_liw -= (maxcol + 10)*IDA_mem->ida_liw1;
1640
1641 if (IDA_mem->ida_VatolMallocDone) {
1642 N_VDestroy(IDA_mem->ida_Vatol); IDA_mem->ida_Vatol = NULL;
1643 IDA_mem->ida_lrw -= IDA_mem->ida_lrw1;
1644 IDA_mem->ida_liw -= IDA_mem->ida_liw1;
1645 }
1646
1647 if (IDA_mem->ida_constraintsMallocDone) {
1648 N_VDestroy(IDA_mem->ida_constraints);
1649 IDA_mem->ida_constraints = NULL;
1650 IDA_mem->ida_lrw -= IDA_mem->ida_lrw1;
1651 IDA_mem->ida_liw -= IDA_mem->ida_liw1;
1652 }
1653
1654 if (IDA_mem->ida_idMallocDone) {
1655 N_VDestroy(IDA_mem->ida_id); IDA_mem->ida_id = NULL;
1656 IDA_mem->ida_lrw -= IDA_mem->ida_lrw1;
1657 IDA_mem->ida_liw -= IDA_mem->ida_liw1;
1658 }
1659
1660 }
1661
1662 /*
1663 * -----------------------------------------------------------------
1664 * Initial setup
1665 * -----------------------------------------------------------------
1666 */
1667
1668 /*
1669 * IDAInitialSetup
1670 *
1671 * This routine is called by IDASolve once at the first step.
1672 * It performs all checks on optional inputs and inputs to
1673 * IDAInit/IDAReInit that could not be done before.
1674 *
1675 * If no error is encountered, IDAInitialSetup returns IDA_SUCCESS.
1676 * Otherwise, it returns an error flag and reported to the error
1677 * handler function.
1678 */
1679
IDAInitialSetup(IDAMem IDA_mem)1680 int IDAInitialSetup(IDAMem IDA_mem)
1681 {
1682 booleantype conOK;
1683 int ier;
1684
1685 /* Test for more vector operations, depending on options */
1686 if (IDA_mem->ida_suppressalg)
1687 if (IDA_mem->ida_phi[0]->ops->nvwrmsnormmask == NULL) {
1688 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAInitialSetup",
1689 MSG_BAD_NVECTOR);
1690 return(IDA_ILL_INPUT);
1691 }
1692
1693 /* Test id vector for legality */
1694 if (IDA_mem->ida_suppressalg && (IDA_mem->ida_id==NULL)){
1695 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAInitialSetup",
1696 MSG_MISSING_ID);
1697 return(IDA_ILL_INPUT);
1698 }
1699
1700 /* Did the user specify tolerances? */
1701 if (IDA_mem->ida_itol == IDA_NN) {
1702 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAInitialSetup",
1703 MSG_NO_TOLS);
1704 return(IDA_ILL_INPUT);
1705 }
1706
1707 /* Set data for efun */
1708 if (IDA_mem->ida_user_efun) IDA_mem->ida_edata = IDA_mem->ida_user_data;
1709 else IDA_mem->ida_edata = IDA_mem;
1710
1711 /* Initial error weight vector */
1712 ier = IDA_mem->ida_efun(IDA_mem->ida_phi[0], IDA_mem->ida_ewt, IDA_mem->ida_edata);
1713 if (ier != 0) {
1714 if (IDA_mem->ida_itol == IDA_WF)
1715 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAInitialSetup",
1716 MSG_FAIL_EWT);
1717 else
1718 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAInitialSetup",
1719 MSG_BAD_EWT);
1720 return(IDA_ILL_INPUT);
1721 }
1722
1723 /* Check to see if y0 satisfies constraints. */
1724 if (IDA_mem->ida_constraintsSet) {
1725 conOK = N_VConstrMask(IDA_mem->ida_constraints, IDA_mem->ida_phi[0], IDA_mem->ida_tempv2);
1726 if (!conOK) {
1727 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDAInitialSetup",
1728 MSG_Y0_FAIL_CONSTR);
1729 return(IDA_ILL_INPUT);
1730 }
1731 }
1732
1733 /* Call linit function if it exists. */
1734 if (IDA_mem->ida_linit != NULL) {
1735 ier = IDA_mem->ida_linit(IDA_mem);
1736 if (ier != 0) {
1737 IDAProcessError(IDA_mem, IDA_LINIT_FAIL, "IDA", "IDAInitialSetup",
1738 MSG_LINIT_FAIL);
1739 return(IDA_LINIT_FAIL);
1740 }
1741 }
1742
1743 /* Initialize the nonlinear solver (must occur after linear solver is initialize) so
1744 * that lsetup and lsolve pointers have been set */
1745 ier = idaNlsInit(IDA_mem);
1746 if (ier != IDA_SUCCESS) {
1747 IDAProcessError(IDA_mem, IDA_NLS_INIT_FAIL, "IDA", "IDAInitialSetup",
1748 MSG_NLS_INIT_FAIL);
1749 return(IDA_NLS_INIT_FAIL);
1750 }
1751
1752 return(IDA_SUCCESS);
1753 }
1754
1755 /*
1756 * IDAEwtSet
1757 *
1758 * This routine is responsible for loading the error weight vector
1759 * ewt, according to itol, as follows:
1760 * (1) ewt[i] = 1 / (rtol * SUNRabs(ycur[i]) + atol), i=0,...,Neq-1
1761 * if itol = IDA_SS
1762 * (2) ewt[i] = 1 / (rtol * SUNRabs(ycur[i]) + atol[i]), i=0,...,Neq-1
1763 * if itol = IDA_SV
1764 *
1765 * IDAEwtSet returns 0 if ewt is successfully set as above to a
1766 * positive vector and -1 otherwise. In the latter case, ewt is
1767 * considered undefined.
1768 *
1769 * All the real work is done in the routines IDAEwtSetSS, IDAEwtSetSV.
1770 */
1771
IDAEwtSet(N_Vector ycur,N_Vector weight,void * data)1772 int IDAEwtSet(N_Vector ycur, N_Vector weight, void *data)
1773 {
1774 IDAMem IDA_mem;
1775 int flag = 0;
1776
1777 /* data points to IDA_mem here */
1778
1779 IDA_mem = (IDAMem) data;
1780
1781 switch(IDA_mem->ida_itol) {
1782 case IDA_SS:
1783 flag = IDAEwtSetSS(IDA_mem, ycur, weight);
1784 break;
1785 case IDA_SV:
1786 flag = IDAEwtSetSV(IDA_mem, ycur, weight);
1787 break;
1788 }
1789 return(flag);
1790 }
1791
1792 /*
1793 * IDAEwtSetSS
1794 *
1795 * This routine sets ewt as decribed above in the case itol=IDA_SS.
1796 * If the absolute tolerance is zero, it tests for non-positive components
1797 * before inverting. IDAEwtSetSS returns 0 if ewt is successfully set to a
1798 * positive vector and -1 otherwise. In the latter case, ewt is considered
1799 * undefined.
1800 */
1801
IDAEwtSetSS(IDAMem IDA_mem,N_Vector ycur,N_Vector weight)1802 static int IDAEwtSetSS(IDAMem IDA_mem, N_Vector ycur, N_Vector weight)
1803 {
1804 N_VAbs(ycur, IDA_mem->ida_tempv1);
1805 N_VScale(IDA_mem->ida_rtol, IDA_mem->ida_tempv1, IDA_mem->ida_tempv1);
1806 N_VAddConst(IDA_mem->ida_tempv1, IDA_mem->ida_Satol, IDA_mem->ida_tempv1);
1807 if (IDA_mem->ida_atolmin0) {
1808 if (N_VMin(IDA_mem->ida_tempv1) <= ZERO) return(-1);
1809 }
1810 N_VInv(IDA_mem->ida_tempv1, weight);
1811 return(0);
1812 }
1813
1814 /*
1815 * IDAEwtSetSV
1816 *
1817 * This routine sets ewt as decribed above in the case itol=IDA_SV.
1818 * If the absolute tolerance is zero, it tests for non-positive components
1819 * before inverting. IDAEwtSetSV returns 0 if ewt is successfully set to a
1820 * positive vector and -1 otherwise. In the latter case, ewt is considered
1821 * undefined.
1822 */
1823
IDAEwtSetSV(IDAMem IDA_mem,N_Vector ycur,N_Vector weight)1824 static int IDAEwtSetSV(IDAMem IDA_mem, N_Vector ycur, N_Vector weight)
1825 {
1826 N_VAbs(ycur, IDA_mem->ida_tempv1);
1827 N_VLinearSum(IDA_mem->ida_rtol, IDA_mem->ida_tempv1,
1828 ONE, IDA_mem->ida_Vatol, IDA_mem->ida_tempv1);
1829 if (IDA_mem->ida_atolmin0) {
1830 if (N_VMin(IDA_mem->ida_tempv1) <= ZERO) return(-1);
1831 }
1832 N_VInv(IDA_mem->ida_tempv1, weight);
1833 return(0);
1834 }
1835
1836 /*
1837 * -----------------------------------------------------------------
1838 * Stopping tests
1839 * -----------------------------------------------------------------
1840 */
1841
1842 /*
1843 * IDAStopTest1
1844 *
1845 * This routine tests for stop conditions before taking a step.
1846 * The tests depend on the value of itask.
1847 * The variable tretlast is the previously returned value of tret.
1848 *
1849 * The return values are:
1850 * CONTINUE_STEPS if no stop conditions were found
1851 * IDA_SUCCESS for a normal return to the user
1852 * IDA_TSTOP_RETURN for a tstop-reached return to the user
1853 * IDA_ILL_INPUT for an illegal-input return to the user
1854 *
1855 * In the tstop cases, this routine may adjust the stepsize hh to cause
1856 * the next step to reach tstop exactly.
1857 */
1858
IDAStopTest1(IDAMem IDA_mem,realtype tout,realtype * tret,N_Vector yret,N_Vector ypret,int itask)1859 static int IDAStopTest1(IDAMem IDA_mem, realtype tout, realtype *tret,
1860 N_Vector yret, N_Vector ypret, int itask)
1861 {
1862 int ier;
1863 realtype troundoff;
1864
1865 switch (itask) {
1866
1867 case IDA_NORMAL:
1868
1869 if (IDA_mem->ida_tstopset) {
1870 /* Test for tn past tstop, tn = tretlast, tn past tout, tn near tstop. */
1871 if ( (IDA_mem->ida_tn - IDA_mem->ida_tstop)*IDA_mem->ida_hh > ZERO) {
1872 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve",
1873 MSG_BAD_TSTOP, IDA_mem->ida_tstop, IDA_mem->ida_tn);
1874 return(IDA_ILL_INPUT);
1875 }
1876 }
1877
1878 /* Test for tout = tretlast, and for tn past tout. */
1879 if (tout == IDA_mem->ida_tretlast) {
1880 *tret = IDA_mem->ida_tretlast = tout;
1881 return(IDA_SUCCESS);
1882 }
1883 if ((IDA_mem->ida_tn - tout)*IDA_mem->ida_hh >= ZERO) {
1884 ier = IDAGetSolution(IDA_mem, tout, yret, ypret);
1885 if (ier != IDA_SUCCESS) {
1886 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve", MSG_BAD_TOUT, tout);
1887 return(IDA_ILL_INPUT);
1888 }
1889 *tret = IDA_mem->ida_tretlast = tout;
1890 return(IDA_SUCCESS);
1891 }
1892
1893 if (IDA_mem->ida_tstopset) {
1894 troundoff = HUNDRED * IDA_mem->ida_uround * (SUNRabs(IDA_mem->ida_tn) + SUNRabs(IDA_mem->ida_hh));
1895 if (SUNRabs(IDA_mem->ida_tn - IDA_mem->ida_tstop) <= troundoff) {
1896 ier = IDAGetSolution(IDA_mem, IDA_mem->ida_tstop, yret, ypret);
1897 if (ier != IDA_SUCCESS) {
1898 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve",
1899 MSG_BAD_TSTOP, IDA_mem->ida_tstop, IDA_mem->ida_tn);
1900 return(IDA_ILL_INPUT);
1901 }
1902 *tret = IDA_mem->ida_tretlast = IDA_mem->ida_tstop;
1903 IDA_mem->ida_tstopset = SUNFALSE;
1904 return(IDA_TSTOP_RETURN);
1905 }
1906 if ((IDA_mem->ida_tn + IDA_mem->ida_hh - IDA_mem->ida_tstop)*IDA_mem->ida_hh > ZERO)
1907 IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn)*(ONE - FOUR * IDA_mem->ida_uround);
1908 }
1909
1910 return(CONTINUE_STEPS);
1911
1912 case IDA_ONE_STEP:
1913
1914 if (IDA_mem->ida_tstopset) {
1915 /* Test for tn past tstop, tn past tretlast, and tn near tstop. */
1916 if ((IDA_mem->ida_tn - IDA_mem->ida_tstop)*IDA_mem->ida_hh > ZERO) {
1917 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve",
1918 MSG_BAD_TSTOP, IDA_mem->ida_tstop, IDA_mem->ida_tn);
1919 return(IDA_ILL_INPUT);
1920 }
1921 }
1922
1923 /* Test for tn past tretlast. */
1924 if ((IDA_mem->ida_tn - IDA_mem->ida_tretlast)*IDA_mem->ida_hh > ZERO) {
1925 ier = IDAGetSolution(IDA_mem, IDA_mem->ida_tn, yret, ypret);
1926 *tret = IDA_mem->ida_tretlast = IDA_mem->ida_tn;
1927 return(IDA_SUCCESS);
1928 }
1929
1930 if (IDA_mem->ida_tstopset) {
1931 troundoff = HUNDRED * IDA_mem->ida_uround * (SUNRabs(IDA_mem->ida_tn) + SUNRabs(IDA_mem->ida_hh));
1932 if (SUNRabs(IDA_mem->ida_tn - IDA_mem->ida_tstop) <= troundoff) {
1933 ier = IDAGetSolution(IDA_mem, IDA_mem->ida_tstop, yret, ypret);
1934 if (ier != IDA_SUCCESS) {
1935 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDA", "IDASolve",
1936 MSG_BAD_TSTOP, IDA_mem->ida_tstop, IDA_mem->ida_tn);
1937 return(IDA_ILL_INPUT);
1938 }
1939 *tret = IDA_mem->ida_tretlast = IDA_mem->ida_tstop;
1940 IDA_mem->ida_tstopset = SUNFALSE;
1941 return(IDA_TSTOP_RETURN);
1942 }
1943 if ((IDA_mem->ida_tn + IDA_mem->ida_hh - IDA_mem->ida_tstop)*IDA_mem->ida_hh > ZERO)
1944 IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn)*(ONE - FOUR * IDA_mem->ida_uround);
1945 }
1946
1947 return(CONTINUE_STEPS);
1948
1949 }
1950 return(IDA_ILL_INPUT); /* This return should never happen. */
1951 }
1952
1953 /*
1954 * IDAStopTest2
1955 *
1956 * This routine tests for stop conditions after taking a step.
1957 * The tests depend on the value of itask.
1958 *
1959 * The return values are:
1960 * CONTINUE_STEPS if no stop conditions were found
1961 * IDA_SUCCESS for a normal return to the user
1962 * IDA_TSTOP_RETURN for a tstop-reached return to the user
1963 * IDA_ILL_INPUT for an illegal-input return to the user
1964 *
1965 * In the two cases with tstop, this routine may reset the stepsize hh
1966 * to cause the next step to reach tstop exactly.
1967 *
1968 * In the two cases with ONE_STEP mode, no interpolation to tn is needed
1969 * because yret and ypret already contain the current y and y' values.
1970 *
1971 * Note: No test is made for an error return from IDAGetSolution here,
1972 * because the same test was made prior to the step.
1973 */
1974
IDAStopTest2(IDAMem IDA_mem,realtype tout,realtype * tret,N_Vector yret,N_Vector ypret,int itask)1975 static int IDAStopTest2(IDAMem IDA_mem, realtype tout, realtype *tret,
1976 N_Vector yret, N_Vector ypret, int itask)
1977 {
1978 /* int ier; */
1979 realtype troundoff;
1980
1981 switch (itask) {
1982
1983 case IDA_NORMAL:
1984
1985 /* Test for tn past tout. */
1986 if ((IDA_mem->ida_tn - tout)*IDA_mem->ida_hh >= ZERO) {
1987 /* ier = */ IDAGetSolution(IDA_mem, tout, yret, ypret);
1988 *tret = IDA_mem->ida_tretlast = tout;
1989 return(IDA_SUCCESS);
1990 }
1991
1992 if (IDA_mem->ida_tstopset) {
1993 /* Test for tn at tstop and for tn near tstop */
1994 troundoff = HUNDRED * IDA_mem->ida_uround * (SUNRabs(IDA_mem->ida_tn) + SUNRabs(IDA_mem->ida_hh));
1995 if (SUNRabs(IDA_mem->ida_tn - IDA_mem->ida_tstop) <= troundoff) {
1996 /* ier = */ IDAGetSolution(IDA_mem, IDA_mem->ida_tstop, yret, ypret);
1997 *tret = IDA_mem->ida_tretlast = IDA_mem->ida_tstop;
1998 IDA_mem->ida_tstopset = SUNFALSE;
1999 return(IDA_TSTOP_RETURN);
2000 }
2001 if ((IDA_mem->ida_tn + IDA_mem->ida_hh - IDA_mem->ida_tstop)*IDA_mem->ida_hh > ZERO)
2002 IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn)*(ONE - FOUR * IDA_mem->ida_uround);
2003 }
2004
2005 return(CONTINUE_STEPS);
2006
2007 case IDA_ONE_STEP:
2008
2009 if (IDA_mem->ida_tstopset) {
2010 /* Test for tn at tstop and for tn near tstop */
2011 troundoff = HUNDRED * IDA_mem->ida_uround * (SUNRabs(IDA_mem->ida_tn) + SUNRabs(IDA_mem->ida_hh));
2012 if (SUNRabs(IDA_mem->ida_tn - IDA_mem->ida_tstop) <= troundoff) {
2013 /* ier = */ IDAGetSolution(IDA_mem, IDA_mem->ida_tstop, yret, ypret);
2014 *tret = IDA_mem->ida_tretlast = IDA_mem->ida_tstop;
2015 IDA_mem->ida_tstopset = SUNFALSE;
2016 return(IDA_TSTOP_RETURN);
2017 }
2018 if ((IDA_mem->ida_tn + IDA_mem->ida_hh - IDA_mem->ida_tstop)*IDA_mem->ida_hh > ZERO)
2019 IDA_mem->ida_hh = (IDA_mem->ida_tstop - IDA_mem->ida_tn)*(ONE - FOUR * IDA_mem->ida_uround);
2020 }
2021
2022 *tret = IDA_mem->ida_tretlast = IDA_mem->ida_tn;
2023 return(IDA_SUCCESS);
2024
2025 }
2026 return IDA_ILL_INPUT; /* This return should never happen. */
2027 }
2028
2029 /*
2030 * -----------------------------------------------------------------
2031 * Error handler
2032 * -----------------------------------------------------------------
2033 */
2034
2035 /*
2036 * IDAHandleFailure
2037 *
2038 * This routine prints error messages for all cases of failure by
2039 * IDAStep. It returns to IDASolve the value that it is to return to
2040 * the user.
2041 */
2042
IDAHandleFailure(IDAMem IDA_mem,int sflag)2043 static int IDAHandleFailure(IDAMem IDA_mem, int sflag)
2044 {
2045 /* Depending on sflag, print error message and return error flag */
2046 switch (sflag) {
2047
2048 case IDA_ERR_FAIL:
2049 IDAProcessError(IDA_mem, IDA_ERR_FAIL, "IDA", "IDASolve",
2050 MSG_ERR_FAILS, IDA_mem->ida_tn, IDA_mem->ida_hh);
2051 return(IDA_ERR_FAIL);
2052
2053 case IDA_CONV_FAIL:
2054 IDAProcessError(IDA_mem, IDA_CONV_FAIL, "IDA", "IDASolve",
2055 MSG_CONV_FAILS, IDA_mem->ida_tn, IDA_mem->ida_hh);
2056 return(IDA_CONV_FAIL);
2057
2058 case IDA_LSETUP_FAIL:
2059 IDAProcessError(IDA_mem, IDA_LSETUP_FAIL, "IDA", "IDASolve",
2060 MSG_SETUP_FAILED, IDA_mem->ida_tn);
2061 return(IDA_LSETUP_FAIL);
2062
2063 case IDA_LSOLVE_FAIL:
2064 IDAProcessError(IDA_mem, IDA_LSOLVE_FAIL, "IDA", "IDASolve",
2065 MSG_SOLVE_FAILED, IDA_mem->ida_tn);
2066 return(IDA_LSOLVE_FAIL);
2067
2068 case IDA_REP_RES_ERR:
2069 IDAProcessError(IDA_mem, IDA_REP_RES_ERR, "IDA", "IDASolve",
2070 MSG_REP_RES_ERR, IDA_mem->ida_tn);
2071 return(IDA_REP_RES_ERR);
2072
2073 case IDA_RES_FAIL:
2074 IDAProcessError(IDA_mem, IDA_RES_FAIL, "IDA", "IDASolve",
2075 MSG_RES_NONRECOV, IDA_mem->ida_tn);
2076 return(IDA_RES_FAIL);
2077
2078 case IDA_CONSTR_FAIL:
2079 IDAProcessError(IDA_mem, IDA_CONSTR_FAIL, "IDA", "IDASolve",
2080 MSG_FAILED_CONSTR, IDA_mem->ida_tn);
2081 return(IDA_CONSTR_FAIL);
2082
2083 case IDA_MEM_NULL:
2084 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDASolve",
2085 MSG_NO_MEM);
2086 return(IDA_MEM_NULL);
2087
2088 case SUN_NLS_MEM_NULL:
2089 IDAProcessError(IDA_mem, IDA_MEM_NULL, "IDA", "IDASolve",
2090 MSG_NLS_INPUT_NULL, IDA_mem->ida_tn);
2091 return(IDA_MEM_NULL);
2092
2093 case IDA_NLS_SETUP_FAIL:
2094 IDAProcessError(IDA_mem, IDA_NLS_SETUP_FAIL, "IDA", "IDASolve",
2095 MSG_NLS_SETUP_FAILED, IDA_mem->ida_tn);
2096 return(IDA_NLS_SETUP_FAIL);
2097 case IDA_NLS_FAIL:
2098 IDAProcessError(IDA_mem, IDA_NLS_FAIL, "IDA", "IDASolve",
2099 MSG_NLS_FAIL, IDA_mem->ida_tn);
2100 return(IDA_NLS_FAIL);
2101 }
2102
2103 /* This return should never happen */
2104 IDAProcessError(IDA_mem, IDA_UNRECOGNIZED_ERROR, "IDA", "IDASolve",
2105 "IDA encountered an unrecognized error. Please report this to the Sundials developers at sundials-users@llnl.gov");
2106 return (IDA_UNRECOGNIZED_ERROR);
2107 }
2108
2109 /*
2110 * -----------------------------------------------------------------
2111 * Main IDAStep function
2112 * -----------------------------------------------------------------
2113 */
2114
2115 /*
2116 * IDAStep
2117 *
2118 * This routine performs one internal IDA step, from tn to tn + hh.
2119 * It calls other routines to do all the work.
2120 *
2121 * It solves a system of differential/algebraic equations of the form
2122 * F(t,y,y') = 0, for one step. In IDA, tt is used for t,
2123 * yy is used for y, and yp is used for y'. The function F is supplied as 'res'
2124 * by the user.
2125 *
2126 * The methods used are modified divided difference, fixed leading
2127 * coefficient forms of backward differentiation formulas.
2128 * The code adjusts the stepsize and order to control the local error per step.
2129 *
2130 * The main operations done here are as follows:
2131 * * initialize various quantities;
2132 * * setting of multistep method coefficients;
2133 * * solution of the nonlinear system for yy at t = tn + hh;
2134 * * deciding on order reduction and testing the local error;
2135 * * attempting to recover from failure in nonlinear solver or error test;
2136 * * resetting stepsize and order for the next step.
2137 * * updating phi and other state data if successful;
2138 *
2139 * On a failure in the nonlinear system solution or error test, the
2140 * step may be reattempted, depending on the nature of the failure.
2141 *
2142 * Variables or arrays (all in the IDAMem structure) used in IDAStep are:
2143 *
2144 * tt -- Independent variable.
2145 * yy -- Solution vector at tt.
2146 * yp -- Derivative of solution vector after successful stelp.
2147 * res -- User-supplied function to evaluate the residual. See the
2148 * description given in file ida.h .
2149 * lsetup -- Routine to prepare for the linear solver call. It may either
2150 * save or recalculate quantities used by lsolve. (Optional)
2151 * lsolve -- Routine to solve a linear system. A prior call to lsetup
2152 * may be required.
2153 * hh -- Appropriate step size for next step.
2154 * ewt -- Vector of weights used in all convergence tests.
2155 * phi -- Array of divided differences used by IDAStep. This array is composed
2156 * of (maxord+1) nvectors (each of size Neq). (maxord+1) is the maximum
2157 * order for the problem, maxord, plus 1.
2158 *
2159 * Return values are:
2160 * IDA_SUCCESS IDA_RES_FAIL LSETUP_ERROR_NONRECVR
2161 * IDA_LSOLVE_FAIL IDA_ERR_FAIL
2162 * IDA_CONSTR_FAIL IDA_CONV_FAIL
2163 * IDA_REP_RES_ERR
2164 */
2165
IDAStep(IDAMem IDA_mem)2166 static int IDAStep(IDAMem IDA_mem)
2167 {
2168 realtype saved_t, ck;
2169 realtype err_k, err_km1;
2170 int ncf, nef;
2171 int nflag, kflag;
2172
2173 saved_t = IDA_mem->ida_tn;
2174 ncf = nef = 0;
2175
2176 if (IDA_mem->ida_nst == ZERO){
2177 IDA_mem->ida_kk = 1;
2178 IDA_mem->ida_kused = 0;
2179 IDA_mem->ida_hused = ZERO;
2180 IDA_mem->ida_psi[0] = IDA_mem->ida_hh;
2181 IDA_mem->ida_cj = ONE/IDA_mem->ida_hh;
2182 IDA_mem->ida_phase = 0;
2183 IDA_mem->ida_ns = 0;
2184 }
2185
2186 /* To prevent 'unintialized variable' warnings */
2187 err_k = ZERO;
2188 err_km1 = ZERO;
2189
2190 /* Looping point for attempts to take a step */
2191
2192 for(;;) {
2193
2194 /*-----------------------
2195 Set method coefficients
2196 -----------------------*/
2197
2198 IDASetCoeffs(IDA_mem, &ck);
2199
2200 kflag = IDA_SUCCESS;
2201
2202 /*----------------------------------------------------
2203 If tn is past tstop (by roundoff), reset it to tstop.
2204 -----------------------------------------------------*/
2205
2206 IDA_mem->ida_tn = IDA_mem->ida_tn + IDA_mem->ida_hh;
2207 if (IDA_mem->ida_tstopset) {
2208 if ((IDA_mem->ida_tn - IDA_mem->ida_tstop)*IDA_mem->ida_hh > ZERO)
2209 IDA_mem->ida_tn = IDA_mem->ida_tstop;
2210 }
2211
2212 /*-----------------------
2213 Advance state variables
2214 -----------------------*/
2215
2216 /* Compute predicted values for yy and yp */
2217 IDAPredict(IDA_mem);
2218
2219 /* Nonlinear system solution */
2220 nflag = IDANls(IDA_mem);
2221
2222 /* If NLS was successful, perform error test */
2223 if (nflag == IDA_SUCCESS)
2224 nflag = IDATestError(IDA_mem, ck, &err_k, &err_km1);
2225
2226 /* Test for convergence or error test failures */
2227 if (nflag != IDA_SUCCESS) {
2228
2229 /* restore and decide what to do */
2230 IDARestore(IDA_mem, saved_t);
2231 kflag = IDAHandleNFlag(IDA_mem, nflag, err_k, err_km1,
2232 &(IDA_mem->ida_ncfn), &ncf,
2233 &(IDA_mem->ida_netf), &nef);
2234
2235 /* exit on nonrecoverable failure */
2236 if (kflag != PREDICT_AGAIN) return(kflag);
2237
2238 /* recoverable error; predict again */
2239 if(IDA_mem->ida_nst==0) IDAReset(IDA_mem);
2240 continue;
2241
2242 }
2243
2244 /* kflag == IDA_SUCCESS */
2245 break;
2246
2247 } /* end loop */
2248
2249 /* Nonlinear system solve and error test were both successful;
2250 update data, and consider change of step and/or order */
2251
2252 IDACompleteStep(IDA_mem, err_k, err_km1);
2253
2254 /*
2255 Rescale ee vector to be the estimated local error
2256 Notes:
2257 (1) altering the value of ee is permissible since
2258 it will be overwritten by
2259 IDASolve()->IDAStep()->IDANls()
2260 before it is needed again
2261 (2) the value of ee is only valid if IDAHandleNFlag()
2262 returns either PREDICT_AGAIN or IDA_SUCCESS
2263 */
2264
2265 N_VScale(ck, IDA_mem->ida_ee, IDA_mem->ida_ee);
2266
2267 return(IDA_SUCCESS);
2268 }
2269
2270 /*
2271 * IDASetCoeffs
2272 *
2273 * This routine computes the coefficients relevant to the current step.
2274 * The counter ns counts the number of consecutive steps taken at
2275 * constant stepsize h and order k, up to a maximum of k + 2.
2276 * Then the first ns components of beta will be one, and on a step
2277 * with ns = k + 2, the coefficients alpha, etc. need not be reset here.
2278 * Also, IDACompleteStep prohibits an order increase until ns = k + 2.
2279 */
2280
IDASetCoeffs(IDAMem IDA_mem,realtype * ck)2281 static void IDASetCoeffs(IDAMem IDA_mem, realtype *ck)
2282 {
2283 int i;
2284 realtype temp1, temp2, alpha0, alphas;
2285
2286 /* Set coefficients for the current stepsize h */
2287
2288 if ( (IDA_mem->ida_hh != IDA_mem->ida_hused) ||
2289 (IDA_mem->ida_kk != IDA_mem->ida_kused) )
2290 IDA_mem->ida_ns = 0;
2291 IDA_mem->ida_ns = SUNMIN(IDA_mem->ida_ns+1, IDA_mem->ida_kused+2);
2292 if (IDA_mem->ida_kk + 1 >= IDA_mem->ida_ns) {
2293 IDA_mem->ida_beta[0] = ONE;
2294 IDA_mem->ida_alpha[0] = ONE;
2295 temp1 = IDA_mem->ida_hh;
2296 IDA_mem->ida_gamma[0] = ZERO;
2297 IDA_mem->ida_sigma[0] = ONE;
2298 for(i=1; i<=IDA_mem->ida_kk; i++) {
2299 temp2 = IDA_mem->ida_psi[i-1];
2300 IDA_mem->ida_psi[i-1] = temp1;
2301 IDA_mem->ida_beta[i] = IDA_mem->ida_beta[i-1] * IDA_mem->ida_psi[i-1] / temp2;
2302 temp1 = temp2 + IDA_mem->ida_hh;
2303 IDA_mem->ida_alpha[i] = IDA_mem->ida_hh / temp1;
2304 IDA_mem->ida_sigma[i] = i * IDA_mem->ida_sigma[i-1] * IDA_mem->ida_alpha[i];
2305 IDA_mem->ida_gamma[i] = IDA_mem->ida_gamma[i-1] + IDA_mem->ida_alpha[i-1] / IDA_mem->ida_hh;
2306 }
2307 IDA_mem->ida_psi[IDA_mem->ida_kk] = temp1;
2308 }
2309 /* compute alphas, alpha0 */
2310 alphas = ZERO;
2311 alpha0 = ZERO;
2312 for(i=0; i<IDA_mem->ida_kk; i++) {
2313 alphas = alphas - ONE/(i+1);
2314 alpha0 = alpha0 - IDA_mem->ida_alpha[i];
2315 }
2316
2317 /* compute leading coefficient cj */
2318 IDA_mem->ida_cjlast = IDA_mem->ida_cj;
2319 IDA_mem->ida_cj = -alphas/IDA_mem->ida_hh;
2320
2321 /* compute variable stepsize error coefficient ck */
2322
2323 *ck = SUNRabs(IDA_mem->ida_alpha[IDA_mem->ida_kk] + alphas - alpha0);
2324 *ck = SUNMAX(*ck, IDA_mem->ida_alpha[IDA_mem->ida_kk]);
2325
2326 /* change phi to phi-star */
2327
2328 /* Scale i=IDA_mem->ida_ns to i<=IDA_mem->ida_kk */
2329 if (IDA_mem->ida_ns <= IDA_mem->ida_kk) {
2330 (void) N_VScaleVectorArray(IDA_mem->ida_kk - IDA_mem->ida_ns + 1,
2331 IDA_mem->ida_beta+IDA_mem->ida_ns,
2332 IDA_mem->ida_phi+IDA_mem->ida_ns,
2333 IDA_mem->ida_phi+IDA_mem->ida_ns);
2334 }
2335
2336 }
2337
2338 /*
2339 * -----------------------------------------------------------------
2340 * Nonlinear solver functions
2341 * -----------------------------------------------------------------
2342 */
2343
2344 /*
2345 * IDANls
2346 *
2347 * This routine attempts to solve the nonlinear system using the linear
2348 * solver specified. NOTE: this routine uses N_Vector ee as the scratch
2349 * vector tempv3 passed to lsetup.
2350 *
2351 * Possible return values:
2352 *
2353 * IDA_SUCCESS
2354 *
2355 * IDA_RES_RECVR IDA_RES_FAIL
2356 * IDA_LSETUP_RECVR IDA_LSETUP_FAIL
2357 * IDA_LSOLVE_RECVR IDA_LSOLVE_FAIL
2358 *
2359 * IDA_CONSTR_RECVR
2360 * SUN_NLS_CONV_RECVR
2361 * IDA_MEM_NULL
2362 */
2363
IDANls(IDAMem IDA_mem)2364 static int IDANls(IDAMem IDA_mem)
2365 {
2366 int retval;
2367 booleantype constraintsPassed, callLSetup;
2368 realtype temp1, temp2, vnorm;
2369 N_Vector mm, tmp;
2370
2371 callLSetup = SUNFALSE;
2372
2373 /* Initialize if the first time called */
2374
2375 if (IDA_mem->ida_nst == 0){
2376 IDA_mem->ida_cjold = IDA_mem->ida_cj;
2377 IDA_mem->ida_ss = TWENTY;
2378 if (IDA_mem->ida_lsetup) callLSetup = SUNTRUE;
2379 }
2380
2381 /* Decide if lsetup is to be called */
2382
2383 if (IDA_mem->ida_lsetup) {
2384 IDA_mem->ida_cjratio = IDA_mem->ida_cj / IDA_mem->ida_cjold;
2385 temp1 = (ONE - XRATE) / (ONE + XRATE);
2386 temp2 = ONE/temp1;
2387 if (IDA_mem->ida_cjratio < temp1 || IDA_mem->ida_cjratio > temp2) callLSetup = SUNTRUE;
2388 if (IDA_mem->ida_cj != IDA_mem->ida_cjlast) IDA_mem->ida_ss = HUNDRED;
2389 }
2390
2391 /* initial guess for the correction to the predictor */
2392 N_VConst(ZERO, IDA_mem->ida_ee);
2393
2394 /* call nonlinear solver setup if it exists */
2395 if ((IDA_mem->NLS)->ops->setup) {
2396 retval = SUNNonlinSolSetup(IDA_mem->NLS, IDA_mem->ida_ee, IDA_mem);
2397 if (retval < 0) return(IDA_NLS_SETUP_FAIL);
2398 if (retval > 0) return(IDA_NLS_SETUP_RECVR);
2399 }
2400
2401 /* solve the nonlinear system */
2402 retval = SUNNonlinSolSolve(IDA_mem->NLS,
2403 IDA_mem->ida_yypredict, IDA_mem->ida_ee,
2404 IDA_mem->ida_ewt, IDA_mem->ida_epsNewt,
2405 callLSetup, IDA_mem);
2406
2407 /* update yy and yp based on the final correction from the nonlinear solver */
2408 N_VLinearSum(ONE, IDA_mem->ida_yypredict, ONE, IDA_mem->ida_ee, IDA_mem->ida_yy);
2409 N_VLinearSum(ONE, IDA_mem->ida_yppredict, IDA_mem->ida_cj, IDA_mem->ida_ee, IDA_mem->ida_yp);
2410
2411 /* return if nonlinear solver failed */
2412 if (retval != IDA_SUCCESS) return(retval);
2413
2414 /* If otherwise successful, check and enforce inequality constraints. */
2415
2416 if (IDA_mem->ida_constraintsSet) {
2417
2418 /* shortcut names for temporary work vectors */
2419 mm = IDA_mem->ida_tempv2;
2420 tmp = IDA_mem->ida_tempv1;
2421
2422 /* Get mask vector mm, set where constraints failed */
2423 constraintsPassed = N_VConstrMask(IDA_mem->ida_constraints,
2424 IDA_mem->ida_yy, mm);
2425 if (constraintsPassed) return(IDA_SUCCESS);
2426
2427 /* Constraints not met */
2428
2429 /* Compute correction to satisfy constraints */
2430 N_VCompare(ONEPT5, IDA_mem->ida_constraints, tmp); /* a[i] =1 when |c[i]| = 2 */
2431 N_VProd(tmp, IDA_mem->ida_constraints, tmp); /* a * c */
2432 N_VDiv(tmp, IDA_mem->ida_ewt, tmp); /* a * c * wt */
2433 N_VLinearSum(ONE, IDA_mem->ida_yy, -PT1, tmp, tmp); /* y - 0.1 * a * c * wt */
2434 N_VProd(tmp, mm, tmp); /* v = mm*(y-.1*a*c*wt) */
2435
2436 vnorm = IDAWrmsNorm(IDA_mem, tmp, IDA_mem->ida_ewt, SUNFALSE); /* ||v|| */
2437
2438 /* If vector v of constraint corrections is small in norm, correct and
2439 accept this step */
2440 if (vnorm <= IDA_mem->ida_epsNewt) {
2441 N_VLinearSum(ONE, IDA_mem->ida_ee,
2442 -ONE, tmp, IDA_mem->ida_ee); /* ee <- ee - v */
2443 return(IDA_SUCCESS);
2444 }
2445
2446 /* Constraints correction is too large, reduce h by computing rr = h'/h */
2447 N_VLinearSum(ONE, IDA_mem->ida_phi[0], -ONE, IDA_mem->ida_yy, tmp);
2448 N_VProd(mm, tmp, tmp);
2449 IDA_mem->ida_rr = PT9*N_VMinQuotient(IDA_mem->ida_phi[0], tmp);
2450 IDA_mem->ida_rr = SUNMAX(IDA_mem->ida_rr, PT1);
2451
2452 /* Reattempt step with new step size */
2453 return(IDA_CONSTR_RECVR);
2454 }
2455
2456 return(IDA_SUCCESS);
2457 }
2458
2459
2460 /*
2461 * IDAPredict
2462 *
2463 * This routine predicts the new values for vectors yy and yp.
2464 */
2465
IDAPredict(IDAMem IDA_mem)2466 static void IDAPredict(IDAMem IDA_mem)
2467 {
2468 int j;
2469
2470 for(j=0; j<=IDA_mem->ida_kk; j++)
2471 IDA_mem->ida_cvals[j] = ONE;
2472
2473 (void) N_VLinearCombination(IDA_mem->ida_kk+1, IDA_mem->ida_cvals,
2474 IDA_mem->ida_phi, IDA_mem->ida_yypredict);
2475
2476 (void) N_VLinearCombination(IDA_mem->ida_kk, IDA_mem->ida_gamma+1,
2477 IDA_mem->ida_phi+1, IDA_mem->ida_yppredict);
2478 }
2479
2480 /*
2481 * -----------------------------------------------------------------
2482 * Error test
2483 * -----------------------------------------------------------------
2484 */
2485
2486 /*
2487 * IDATestError
2488 *
2489 * This routine estimates errors at orders k, k-1, k-2, decides
2490 * whether or not to suggest an order decrease, and performs
2491 * the local error test.
2492 *
2493 * IDATestError returns either IDA_SUCCESS or ERROR_TEST_FAIL.
2494 */
2495
IDATestError(IDAMem IDA_mem,realtype ck,realtype * err_k,realtype * err_km1)2496 static int IDATestError(IDAMem IDA_mem, realtype ck,
2497 realtype *err_k, realtype *err_km1)
2498 {
2499 realtype err_km2; /* estimated error at k-2 */
2500 realtype enorm_k, enorm_km1, enorm_km2; /* error norms */
2501 realtype terr_k, terr_km1, terr_km2; /* local truncation error norms */
2502
2503 /* Compute error for order k. */
2504 enorm_k = IDAWrmsNorm(IDA_mem, IDA_mem->ida_ee, IDA_mem->ida_ewt, IDA_mem->ida_suppressalg);
2505 *err_k = IDA_mem->ida_sigma[IDA_mem->ida_kk] * enorm_k;
2506 terr_k = (IDA_mem->ida_kk + 1) * (*err_k);
2507
2508 IDA_mem->ida_knew = IDA_mem->ida_kk;
2509
2510 if ( IDA_mem->ida_kk > 1 ) {
2511
2512 /* Compute error at order k-1 */
2513 N_VLinearSum(ONE, IDA_mem->ida_phi[IDA_mem->ida_kk], ONE, IDA_mem->ida_ee, IDA_mem->ida_delta);
2514 enorm_km1 = IDAWrmsNorm(IDA_mem, IDA_mem->ida_delta,
2515 IDA_mem->ida_ewt, IDA_mem->ida_suppressalg);
2516 *err_km1 = IDA_mem->ida_sigma[IDA_mem->ida_kk - 1] * enorm_km1;
2517 terr_km1 = IDA_mem->ida_kk * (*err_km1);
2518
2519 if ( IDA_mem->ida_kk > 2 ) {
2520
2521 /* Compute error at order k-2 */
2522 N_VLinearSum(ONE, IDA_mem->ida_phi[IDA_mem->ida_kk - 1], ONE,
2523 IDA_mem->ida_delta, IDA_mem->ida_delta);
2524 enorm_km2 = IDAWrmsNorm(IDA_mem, IDA_mem->ida_delta,
2525 IDA_mem->ida_ewt, IDA_mem->ida_suppressalg);
2526 err_km2 = IDA_mem->ida_sigma[IDA_mem->ida_kk - 2] * enorm_km2;
2527 terr_km2 = (IDA_mem->ida_kk - 1) * err_km2;
2528
2529 /* Decrease order if errors are reduced */
2530 if (SUNMAX(terr_km1, terr_km2) <= terr_k)
2531 IDA_mem->ida_knew = IDA_mem->ida_kk - 1;
2532
2533 } else {
2534
2535 /* Decrease order to 1 if errors are reduced by at least 1/2 */
2536 if (terr_km1 <= (HALF * terr_k) )
2537 IDA_mem->ida_knew = IDA_mem->ida_kk - 1;
2538
2539 }
2540
2541 }
2542
2543 /* Perform error test */
2544 if (ck * enorm_k > ONE) return(ERROR_TEST_FAIL);
2545 else return(IDA_SUCCESS);
2546 }
2547
2548 /*
2549 * IDARestore
2550 *
2551 * This routine restores tn, psi, and phi in the event of a failure.
2552 * It changes back phi-star to phi (changed in IDASetCoeffs)
2553 */
2554
IDARestore(IDAMem IDA_mem,realtype saved_t)2555 static void IDARestore(IDAMem IDA_mem, realtype saved_t)
2556 {
2557 int j;
2558
2559 IDA_mem->ida_tn = saved_t;
2560
2561 for (j = 1; j <= IDA_mem->ida_kk; j++)
2562 IDA_mem->ida_psi[j-1] = IDA_mem->ida_psi[j] - IDA_mem->ida_hh;
2563
2564 if (IDA_mem->ida_ns <= IDA_mem->ida_kk) {
2565
2566 for (j = IDA_mem->ida_ns; j <= IDA_mem->ida_kk; j++)
2567 IDA_mem->ida_cvals[j-IDA_mem->ida_ns] = ONE/IDA_mem->ida_beta[j];
2568
2569 (void) N_VScaleVectorArray(IDA_mem->ida_kk - IDA_mem->ida_ns + 1,
2570 IDA_mem->ida_cvals,
2571 IDA_mem->ida_phi+IDA_mem->ida_ns,
2572 IDA_mem->ida_phi+IDA_mem->ida_ns);
2573 }
2574
2575 }
2576
2577 /*
2578 * -----------------------------------------------------------------
2579 * Handler for convergence and/or error test failures
2580 * -----------------------------------------------------------------
2581 */
2582
2583 /*
2584 * IDAHandleNFlag
2585 *
2586 * This routine handles failures indicated by the input variable nflag.
2587 * Positive values indicate various recoverable failures while negative
2588 * values indicate nonrecoverable failures. This routine adjusts the
2589 * step size for recoverable failures.
2590 *
2591 * Possible nflag values (input):
2592 *
2593 * --convergence failures--
2594 * IDA_RES_RECVR > 0
2595 * IDA_LSOLVE_RECVR > 0
2596 * IDA_CONSTR_RECVR > 0
2597 * SUN_NLS_CONV_RECVR > 0
2598 * IDA_RES_FAIL < 0
2599 * IDA_LSOLVE_FAIL < 0
2600 * IDA_LSETUP_FAIL < 0
2601 *
2602 * --error test failure--
2603 * ERROR_TEST_FAIL > 0
2604 *
2605 * Possible kflag values (output):
2606 *
2607 * --recoverable--
2608 * PREDICT_AGAIN
2609 *
2610 * --nonrecoverable--
2611 * IDA_CONSTR_FAIL
2612 * IDA_REP_RES_ERR
2613 * IDA_ERR_FAIL
2614 * IDA_CONV_FAIL
2615 * IDA_RES_FAIL
2616 * IDA_LSETUP_FAIL
2617 * IDA_LSOLVE_FAIL
2618 */
2619
IDAHandleNFlag(IDAMem IDA_mem,int nflag,realtype err_k,realtype err_km1,long int * ncfnPtr,int * ncfPtr,long int * netfPtr,int * nefPtr)2620 static int IDAHandleNFlag(IDAMem IDA_mem, int nflag, realtype err_k, realtype err_km1,
2621 long int *ncfnPtr, int *ncfPtr, long int *netfPtr, int *nefPtr)
2622 {
2623 realtype err_knew;
2624
2625 IDA_mem->ida_phase = 1;
2626
2627 if (nflag != ERROR_TEST_FAIL) {
2628
2629 /*-----------------------
2630 Nonlinear solver failed
2631 -----------------------*/
2632
2633 (*ncfPtr)++; /* local counter for convergence failures */
2634 (*ncfnPtr)++; /* global counter for convergence failures */
2635
2636 if (nflag < 0) { /* nonrecoverable failure */
2637
2638 if (nflag == IDA_LSOLVE_FAIL) return(IDA_LSOLVE_FAIL);
2639 else if (nflag == IDA_LSETUP_FAIL) return(IDA_LSETUP_FAIL);
2640 else if (nflag == IDA_RES_FAIL) return(IDA_RES_FAIL);
2641 else return(IDA_NLS_FAIL);
2642
2643 } else { /* recoverable failure */
2644
2645 /* Reduce step size for a new prediction
2646 Note that if nflag=IDA_CONSTR_RECVR then rr was already set in IDANls */
2647 if (nflag != IDA_CONSTR_RECVR) IDA_mem->ida_rr = QUARTER;
2648 IDA_mem->ida_hh *= IDA_mem->ida_rr;
2649
2650 /* Test if there were too many convergence failures */
2651 if (*ncfPtr < IDA_mem->ida_maxncf) return(PREDICT_AGAIN);
2652 else if (nflag == IDA_RES_RECVR) return(IDA_REP_RES_ERR);
2653 else if (nflag == IDA_CONSTR_RECVR) return(IDA_CONSTR_FAIL);
2654 else return(IDA_CONV_FAIL);
2655 }
2656
2657 } else {
2658
2659 /*-----------------
2660 Error Test failed
2661 -----------------*/
2662
2663 (*nefPtr)++; /* local counter for error test failures */
2664 (*netfPtr)++; /* global counter for error test failures */
2665
2666 if (*nefPtr == 1) {
2667
2668 /* On first error test failure, keep current order or lower order by one.
2669 Compute new stepsize based on differences of the solution. */
2670
2671 err_knew = (IDA_mem->ida_kk == IDA_mem->ida_knew) ? err_k : err_km1;
2672
2673 IDA_mem->ida_kk = IDA_mem->ida_knew;
2674 IDA_mem->ida_rr = PT9 * SUNRpowerR( TWO * err_knew + PT0001, -ONE/(IDA_mem->ida_kk + 1) );
2675 IDA_mem->ida_rr = SUNMAX(QUARTER, SUNMIN(PT9,IDA_mem->ida_rr));
2676 IDA_mem->ida_hh *= IDA_mem->ida_rr;
2677 return(PREDICT_AGAIN);
2678
2679 } else if (*nefPtr == 2) {
2680
2681 /* On second error test failure, use current order or decrease order by one.
2682 Reduce stepsize by factor of 1/4. */
2683
2684 IDA_mem->ida_kk = IDA_mem->ida_knew;
2685 IDA_mem->ida_rr = QUARTER;
2686 IDA_mem->ida_hh *= IDA_mem->ida_rr;
2687 return(PREDICT_AGAIN);
2688
2689 } else if (*nefPtr < IDA_mem->ida_maxnef) {
2690
2691 /* On third and subsequent error test failures, set order to 1.
2692 Reduce stepsize by factor of 1/4. */
2693 IDA_mem->ida_kk = 1;
2694 IDA_mem->ida_rr = QUARTER;
2695 IDA_mem->ida_hh *= IDA_mem->ida_rr;
2696 return(PREDICT_AGAIN);
2697
2698 } else {
2699
2700 /* Too many error test failures */
2701 return(IDA_ERR_FAIL);
2702
2703 }
2704
2705 }
2706
2707 }
2708
2709 /*
2710 * IDAReset
2711 *
2712 * This routine is called only if we need to predict again at the
2713 * very first step. In such a case, reset phi[1] and psi[0].
2714 */
2715
IDAReset(IDAMem IDA_mem)2716 static void IDAReset(IDAMem IDA_mem)
2717 {
2718 IDA_mem->ida_psi[0] = IDA_mem->ida_hh;
2719
2720 N_VScale(IDA_mem->ida_rr, IDA_mem->ida_phi[1], IDA_mem->ida_phi[1]);
2721 }
2722
2723 /*
2724 * -----------------------------------------------------------------
2725 * Function called after a successful step
2726 * -----------------------------------------------------------------
2727 */
2728
2729 /*
2730 * IDACompleteStep
2731 *
2732 * This routine completes a successful step. It increments nst,
2733 * saves the stepsize and order used, makes the final selection of
2734 * stepsize and order for the next step, and updates the phi array.
2735 */
2736
IDACompleteStep(IDAMem IDA_mem,realtype err_k,realtype err_km1)2737 static void IDACompleteStep(IDAMem IDA_mem, realtype err_k, realtype err_km1)
2738 {
2739 int j, kdiff, action;
2740 realtype terr_k, terr_km1, terr_kp1;
2741 realtype err_knew, err_kp1;
2742 realtype enorm, tmp, hnew;
2743
2744 IDA_mem->ida_nst++;
2745 kdiff = IDA_mem->ida_kk - IDA_mem->ida_kused;
2746 IDA_mem->ida_kused = IDA_mem->ida_kk;
2747 IDA_mem->ida_hused = IDA_mem->ida_hh;
2748
2749 if ( (IDA_mem->ida_knew == IDA_mem->ida_kk - 1) ||
2750 (IDA_mem->ida_kk == IDA_mem->ida_maxord) )
2751 IDA_mem->ida_phase = 1;
2752
2753 /* For the first few steps, until either a step fails, or the order is
2754 reduced, or the order reaches its maximum, we raise the order and double
2755 the stepsize. During these steps, phase = 0. Thereafter, phase = 1, and
2756 stepsize and order are set by the usual local error algorithm.
2757
2758 Note that, after the first step, the order is not increased, as not all
2759 of the neccessary information is available yet. */
2760
2761 if (IDA_mem->ida_phase == 0) {
2762
2763 if(IDA_mem->ida_nst > 1) {
2764 IDA_mem->ida_kk++;
2765 hnew = TWO * IDA_mem->ida_hh;
2766 if( (tmp = SUNRabs(hnew) * IDA_mem->ida_hmax_inv) > ONE )
2767 hnew /= tmp;
2768 IDA_mem->ida_hh = hnew;
2769 }
2770
2771 } else {
2772
2773 action = UNSET;
2774
2775 /* Set action = LOWER/MAINTAIN/RAISE to specify order decision */
2776
2777 if (IDA_mem->ida_knew == IDA_mem->ida_kk - 1) {action = LOWER; goto takeaction;}
2778 if (IDA_mem->ida_kk == IDA_mem->ida_maxord) {action = MAINTAIN; goto takeaction;}
2779 if ( (IDA_mem->ida_kk + 1 >= IDA_mem->ida_ns ) ||
2780 (kdiff == 1)) {action = MAINTAIN; goto takeaction;}
2781
2782 /* Estimate the error at order k+1, unless already decided to
2783 reduce order, or already using maximum order, or stepsize has not
2784 been constant, or order was just raised. */
2785
2786 N_VLinearSum(ONE, IDA_mem->ida_ee, -ONE,
2787 IDA_mem->ida_phi[IDA_mem->ida_kk + 1], IDA_mem->ida_tempv1);
2788 enorm = IDAWrmsNorm(IDA_mem, IDA_mem->ida_tempv1, IDA_mem->ida_ewt,
2789 IDA_mem->ida_suppressalg);
2790 err_kp1= enorm/(IDA_mem->ida_kk + 2);
2791
2792 /* Choose among orders k-1, k, k+1 using local truncation error norms. */
2793
2794 terr_k = (IDA_mem->ida_kk + 1) * err_k;
2795 terr_kp1 = (IDA_mem->ida_kk + 2) * err_kp1;
2796
2797 if (IDA_mem->ida_kk == 1) {
2798 if (terr_kp1 >= HALF * terr_k) {action = MAINTAIN; goto takeaction;}
2799 else {action = RAISE; goto takeaction;}
2800 } else {
2801 terr_km1 = IDA_mem->ida_kk * err_km1;
2802 if (terr_km1 <= SUNMIN(terr_k, terr_kp1)) {action = LOWER; goto takeaction;}
2803 else if (terr_kp1 >= terr_k) {action = MAINTAIN; goto takeaction;}
2804 else {action = RAISE; goto takeaction;}
2805 }
2806
2807 takeaction:
2808
2809 /* Set the estimated error norm and, on change of order, reset kk. */
2810 if (action == RAISE) { IDA_mem->ida_kk++; err_knew = err_kp1; }
2811 else if (action == LOWER) { IDA_mem->ida_kk--; err_knew = err_km1; }
2812 else { err_knew = err_k; }
2813
2814 /* Compute rr = tentative ratio hnew/hh from error norm estimate.
2815 Reduce hh if rr <= 1, double hh if rr >= 2, else leave hh as is.
2816 If hh is reduced, hnew/hh is restricted to be between .5 and .9. */
2817
2818 hnew = IDA_mem->ida_hh;
2819 IDA_mem->ida_rr = SUNRpowerR( TWO * err_knew + PT0001, -ONE/(IDA_mem->ida_kk + 1) );
2820
2821 if (IDA_mem->ida_rr >= TWO) {
2822 hnew = TWO * IDA_mem->ida_hh;
2823 if( (tmp = SUNRabs(hnew) * IDA_mem->ida_hmax_inv) > ONE )
2824 hnew /= tmp;
2825 } else if (IDA_mem->ida_rr <= ONE ) {
2826 IDA_mem->ida_rr = SUNMAX(HALF, SUNMIN(PT9,IDA_mem->ida_rr));
2827 hnew = IDA_mem->ida_hh * IDA_mem->ida_rr;
2828 }
2829
2830 IDA_mem->ida_hh = hnew;
2831
2832 } /* end of phase if block */
2833
2834 /* Save ee for possible order increase on next step */
2835 if (IDA_mem->ida_kused < IDA_mem->ida_maxord) {
2836 N_VScale(ONE, IDA_mem->ida_ee, IDA_mem->ida_phi[IDA_mem->ida_kused + 1]);
2837 }
2838
2839 /* Update phi arrays */
2840
2841 /* To update phi arrays compute X += Z where */
2842 /* X = [ phi[kused], phi[kused-1], phi[kused-2], ... phi[1] ] */
2843 /* Z = [ ee, phi[kused], phi[kused-1], ... phi[0] ] */
2844
2845 IDA_mem->ida_Zvecs[0] = IDA_mem->ida_ee;
2846 IDA_mem->ida_Xvecs[0] = IDA_mem->ida_phi[IDA_mem->ida_kused];
2847 for (j=1; j<=IDA_mem->ida_kused; j++) {
2848 IDA_mem->ida_Zvecs[j] = IDA_mem->ida_phi[IDA_mem->ida_kused-j+1];
2849 IDA_mem->ida_Xvecs[j] = IDA_mem->ida_phi[IDA_mem->ida_kused-j];
2850 }
2851
2852 (void) N_VLinearSumVectorArray(IDA_mem->ida_kused+1,
2853 ONE, IDA_mem->ida_Xvecs,
2854 ONE, IDA_mem->ida_Zvecs,
2855 IDA_mem->ida_Xvecs);
2856 }
2857
2858 /*
2859 * -----------------------------------------------------------------
2860 * Interpolated output
2861 * -----------------------------------------------------------------
2862 */
2863
2864 /*
2865 * IDAGetSolution
2866 *
2867 * This routine evaluates y(t) and y'(t) as the value and derivative of
2868 * the interpolating polynomial at the independent variable t, and stores
2869 * the results in the vectors yret and ypret. It uses the current
2870 * independent variable value, tn, and the method order last used, kused.
2871 * This function is called by IDASolve with t = tout, t = tn, or t = tstop.
2872 *
2873 * If kused = 0 (no step has been taken), or if t = tn, then the order used
2874 * here is taken to be 1, giving yret = phi[0], ypret = phi[1]/psi[0].
2875 *
2876 * The return values are:
2877 * IDA_SUCCESS if t is legal, or
2878 * IDA_BAD_T if t is not within the interval of the last step taken.
2879 */
2880
IDAGetSolution(void * ida_mem,realtype t,N_Vector yret,N_Vector ypret)2881 int IDAGetSolution(void *ida_mem, realtype t, N_Vector yret, N_Vector ypret)
2882 {
2883 IDAMem IDA_mem;
2884 realtype tfuzz, tp, delt, c, d, gam;
2885 int j, kord, retval;
2886
2887 if (ida_mem == NULL) {
2888 IDAProcessError(NULL, IDA_MEM_NULL, "IDA", "IDAGetSolution", MSG_NO_MEM);
2889 return (IDA_MEM_NULL);
2890 }
2891 IDA_mem = (IDAMem) ida_mem;
2892
2893 /* Check t for legality. Here tn - hused is t_{n-1}. */
2894
2895 tfuzz = HUNDRED * IDA_mem->ida_uround * (SUNRabs(IDA_mem->ida_tn) + SUNRabs(IDA_mem->ida_hh));
2896 if (IDA_mem->ida_hh < ZERO) tfuzz = - tfuzz;
2897 tp = IDA_mem->ida_tn - IDA_mem->ida_hused - tfuzz;
2898 if ((t - tp)*IDA_mem->ida_hh < ZERO) {
2899 IDAProcessError(IDA_mem, IDA_BAD_T, "IDA", "IDAGetSolution", MSG_BAD_T, t,
2900 IDA_mem->ida_tn-IDA_mem->ida_hused, IDA_mem->ida_tn);
2901 return(IDA_BAD_T);
2902 }
2903
2904 /* Initialize kord = (kused or 1). */
2905
2906 kord = IDA_mem->ida_kused;
2907 if (IDA_mem->ida_kused == 0) kord = 1;
2908
2909 /* Accumulate multiples of columns phi[j] into yret and ypret. */
2910
2911 delt = t - IDA_mem->ida_tn;
2912 c = ONE; d = ZERO;
2913 gam = delt / IDA_mem->ida_psi[0];
2914
2915 IDA_mem->ida_cvals[0] = c;
2916 for (j=1; j <= kord; j++) {
2917 d = d*gam + c / IDA_mem->ida_psi[j-1];
2918 c = c*gam;
2919 gam = (delt + IDA_mem->ida_psi[j-1]) / IDA_mem->ida_psi[j];
2920
2921 IDA_mem->ida_cvals[j] = c;
2922 IDA_mem->ida_dvals[j-1] = d;
2923 }
2924
2925 retval = N_VLinearCombination(kord+1, IDA_mem->ida_cvals,
2926 IDA_mem->ida_phi, yret);
2927 if (retval != IDA_SUCCESS) return(IDA_VECTOROP_ERR);
2928
2929 retval = N_VLinearCombination(kord, IDA_mem->ida_dvals,
2930 IDA_mem->ida_phi+1, ypret);
2931 if (retval != IDA_SUCCESS) return(IDA_VECTOROP_ERR);
2932
2933 return(IDA_SUCCESS);
2934 }
2935
2936 /*
2937 * -----------------------------------------------------------------
2938 * Norm function
2939 * -----------------------------------------------------------------
2940 */
2941
2942 /*
2943 * IDAWrmsNorm
2944 *
2945 * Returns the WRMS norm of vector x with weights w.
2946 * If mask = SUNTRUE, the weight vector w is masked by id, i.e.,
2947 * nrm = N_VWrmsNormMask(x,w,id);
2948 * Otherwise,
2949 * nrm = N_VWrmsNorm(x,w);
2950 *
2951 * mask = SUNFALSE when the call is made from the nonlinear solver.
2952 * mask = suppressalg otherwise.
2953 */
2954
IDAWrmsNorm(IDAMem IDA_mem,N_Vector x,N_Vector w,booleantype mask)2955 realtype IDAWrmsNorm(IDAMem IDA_mem, N_Vector x, N_Vector w,
2956 booleantype mask)
2957 {
2958 realtype nrm;
2959
2960 if (mask) nrm = N_VWrmsNormMask(x, w, IDA_mem->ida_id);
2961 else nrm = N_VWrmsNorm(x, w);
2962
2963 return(nrm);
2964 }
2965
2966 /*
2967 * -----------------------------------------------------------------
2968 * Functions for rootfinding
2969 * -----------------------------------------------------------------
2970 */
2971
2972 /*
2973 * IDARcheck1
2974 *
2975 * This routine completes the initialization of rootfinding memory
2976 * information, and checks whether g has a zero both at and very near
2977 * the initial point of the IVP.
2978 *
2979 * This routine returns an int equal to:
2980 * IDA_RTFUNC_FAIL < 0 if the g function failed, or
2981 * IDA_SUCCESS = 0 otherwise.
2982 */
2983
IDARcheck1(IDAMem IDA_mem)2984 static int IDARcheck1(IDAMem IDA_mem)
2985 {
2986 int i, retval;
2987 realtype smallh, hratio, tplus;
2988 booleantype zroot;
2989
2990 for (i = 0; i < IDA_mem->ida_nrtfn; i++)
2991 IDA_mem->ida_iroots[i] = 0;
2992 IDA_mem->ida_tlo = IDA_mem->ida_tn;
2993 IDA_mem->ida_ttol = ((SUNRabs(IDA_mem->ida_tn) + SUNRabs(IDA_mem->ida_hh)) *
2994 IDA_mem->ida_uround * HUNDRED);
2995
2996 /* Evaluate g at initial t and check for zero values. */
2997 retval = IDA_mem->ida_gfun(IDA_mem->ida_tlo, IDA_mem->ida_phi[0], IDA_mem->ida_phi[1],
2998 IDA_mem->ida_glo, IDA_mem->ida_user_data);
2999 IDA_mem->ida_nge = 1;
3000 if (retval != 0) return(IDA_RTFUNC_FAIL);
3001
3002 zroot = SUNFALSE;
3003 for (i = 0; i < IDA_mem->ida_nrtfn; i++) {
3004 if (SUNRabs(IDA_mem->ida_glo[i]) == ZERO) {
3005 zroot = SUNTRUE;
3006 IDA_mem->ida_gactive[i] = SUNFALSE;
3007 }
3008 }
3009 if (!zroot) return(IDA_SUCCESS);
3010
3011 /* Some g_i is zero at t0; look at g at t0+(small increment). */
3012 hratio = SUNMAX(IDA_mem->ida_ttol / SUNRabs(IDA_mem->ida_hh), PT1);
3013 smallh = hratio * IDA_mem->ida_hh;
3014 tplus = IDA_mem->ida_tlo + smallh;
3015 N_VLinearSum(ONE, IDA_mem->ida_phi[0], smallh, IDA_mem->ida_phi[1], IDA_mem->ida_yy);
3016 retval = IDA_mem->ida_gfun(tplus, IDA_mem->ida_yy, IDA_mem->ida_phi[1],
3017 IDA_mem->ida_ghi, IDA_mem->ida_user_data);
3018 IDA_mem->ida_nge++;
3019 if (retval != 0) return(IDA_RTFUNC_FAIL);
3020
3021 /* We check now only the components of g which were exactly 0.0 at t0
3022 * to see if we can 'activate' them. */
3023 for (i = 0; i < IDA_mem->ida_nrtfn; i++) {
3024 if (!IDA_mem->ida_gactive[i] && SUNRabs(IDA_mem->ida_ghi[i]) != ZERO) {
3025 IDA_mem->ida_gactive[i] = SUNTRUE;
3026 IDA_mem->ida_glo[i] = IDA_mem->ida_ghi[i];
3027 }
3028 }
3029 return(IDA_SUCCESS);
3030 }
3031
3032 /*
3033 * IDARcheck2
3034 *
3035 * This routine checks for exact zeros of g at the last root found,
3036 * if the last return was a root. It then checks for a close pair of
3037 * zeros (an error condition), and for a new root at a nearby point.
3038 * The array glo = g(tlo) at the left endpoint of the search interval
3039 * is adjusted if necessary to assure that all g_i are nonzero
3040 * there, before returning to do a root search in the interval.
3041 *
3042 * On entry, tlo = tretlast is the last value of tret returned by
3043 * IDASolve. This may be the previous tn, the previous tout value,
3044 * or the last root location.
3045 *
3046 * This routine returns an int equal to:
3047 * IDA_RTFUNC_FAIL < 0 if the g function failed, or
3048 * CLOSERT = 3 if a close pair of zeros was found, or
3049 * RTFOUND = 1 if a new zero of g was found near tlo, or
3050 * IDA_SUCCESS = 0 otherwise.
3051 */
3052
IDARcheck2(IDAMem IDA_mem)3053 static int IDARcheck2(IDAMem IDA_mem)
3054 {
3055 int i, retval;
3056 realtype smallh, hratio, tplus;
3057 booleantype zroot;
3058
3059 if (IDA_mem->ida_irfnd == 0) return(IDA_SUCCESS);
3060
3061 (void) IDAGetSolution(IDA_mem, IDA_mem->ida_tlo, IDA_mem->ida_yy, IDA_mem->ida_yp);
3062 retval = IDA_mem->ida_gfun(IDA_mem->ida_tlo, IDA_mem->ida_yy, IDA_mem->ida_yp,
3063 IDA_mem->ida_glo, IDA_mem->ida_user_data);
3064 IDA_mem->ida_nge++;
3065 if (retval != 0) return(IDA_RTFUNC_FAIL);
3066
3067 zroot = SUNFALSE;
3068 for (i = 0; i < IDA_mem->ida_nrtfn; i++)
3069 IDA_mem->ida_iroots[i] = 0;
3070 for (i = 0; i < IDA_mem->ida_nrtfn; i++) {
3071 if (!IDA_mem->ida_gactive[i]) continue;
3072 if (SUNRabs(IDA_mem->ida_glo[i]) == ZERO) {
3073 zroot = SUNTRUE;
3074 IDA_mem->ida_iroots[i] = 1;
3075 }
3076 }
3077 if (!zroot) return(IDA_SUCCESS);
3078
3079 /* One or more g_i has a zero at tlo. Check g at tlo+smallh. */
3080 IDA_mem->ida_ttol = ((SUNRabs(IDA_mem->ida_tn) + SUNRabs(IDA_mem->ida_hh)) *
3081 IDA_mem->ida_uround * HUNDRED);
3082 smallh = (IDA_mem->ida_hh > ZERO) ? IDA_mem->ida_ttol : -IDA_mem->ida_ttol;
3083 tplus = IDA_mem->ida_tlo + smallh;
3084 if ( (tplus - IDA_mem->ida_tn)*IDA_mem->ida_hh >= ZERO) {
3085 hratio = smallh/IDA_mem->ida_hh;
3086 N_VLinearSum(ONE, IDA_mem->ida_yy,
3087 hratio, IDA_mem->ida_phi[1], IDA_mem->ida_yy);
3088 } else {
3089 (void) IDAGetSolution(IDA_mem, tplus, IDA_mem->ida_yy, IDA_mem->ida_yp);
3090 }
3091 retval = IDA_mem->ida_gfun(tplus, IDA_mem->ida_yy, IDA_mem->ida_yp,
3092 IDA_mem->ida_ghi, IDA_mem->ida_user_data);
3093 IDA_mem->ida_nge++;
3094 if (retval != 0) return(IDA_RTFUNC_FAIL);
3095
3096 /* Check for close roots (error return), for a new zero at tlo+smallh,
3097 and for a g_i that changed from zero to nonzero. */
3098 zroot = SUNFALSE;
3099 for (i = 0; i < IDA_mem->ida_nrtfn; i++) {
3100 if (!IDA_mem->ida_gactive[i]) continue;
3101 if (SUNRabs(IDA_mem->ida_ghi[i]) == ZERO) {
3102 if (IDA_mem->ida_iroots[i] == 1) return(CLOSERT);
3103 zroot = SUNTRUE;
3104 IDA_mem->ida_iroots[i] = 1;
3105 } else {
3106 if (IDA_mem->ida_iroots[i] == 1)
3107 IDA_mem->ida_glo[i] = IDA_mem->ida_ghi[i];
3108 }
3109 }
3110 if (zroot) return(RTFOUND);
3111 return(IDA_SUCCESS);
3112 }
3113
3114 /*
3115 * IDARcheck3
3116 *
3117 * This routine interfaces to IDARootfind to look for a root of g
3118 * between tlo and either tn or tout, whichever comes first.
3119 * Only roots beyond tlo in the direction of integration are sought.
3120 *
3121 * This routine returns an int equal to:
3122 * IDA_RTFUNC_FAIL < 0 if the g function failed, or
3123 * RTFOUND = 1 if a root of g was found, or
3124 * IDA_SUCCESS = 0 otherwise.
3125 */
3126
IDARcheck3(IDAMem IDA_mem)3127 static int IDARcheck3(IDAMem IDA_mem)
3128 {
3129 int i, ier, retval;
3130
3131 /* Set thi = tn or tout, whichever comes first. */
3132 if (IDA_mem->ida_taskc == IDA_ONE_STEP) IDA_mem->ida_thi = IDA_mem->ida_tn;
3133 if (IDA_mem->ida_taskc == IDA_NORMAL) {
3134 IDA_mem->ida_thi = ((IDA_mem->ida_toutc - IDA_mem->ida_tn)*IDA_mem->ida_hh >= ZERO)
3135 ? IDA_mem->ida_tn : IDA_mem->ida_toutc;
3136 }
3137
3138 /* Get y and y' at thi. */
3139 (void) IDAGetSolution(IDA_mem, IDA_mem->ida_thi, IDA_mem->ida_yy, IDA_mem->ida_yp);
3140
3141
3142 /* Set ghi = g(thi) and call IDARootfind to search (tlo,thi) for roots. */
3143 retval = IDA_mem->ida_gfun(IDA_mem->ida_thi, IDA_mem->ida_yy,
3144 IDA_mem->ida_yp, IDA_mem->ida_ghi,
3145 IDA_mem->ida_user_data);
3146 IDA_mem->ida_nge++;
3147 if (retval != 0) return(IDA_RTFUNC_FAIL);
3148
3149 IDA_mem->ida_ttol = ((SUNRabs(IDA_mem->ida_tn) + SUNRabs(IDA_mem->ida_hh)) *
3150 IDA_mem->ida_uround * HUNDRED);
3151 ier = IDARootfind(IDA_mem);
3152 if (ier == IDA_RTFUNC_FAIL) return(IDA_RTFUNC_FAIL);
3153 for(i=0; i<IDA_mem->ida_nrtfn; i++) {
3154 if(!IDA_mem->ida_gactive[i] && IDA_mem->ida_grout[i] != ZERO)
3155 IDA_mem->ida_gactive[i] = SUNTRUE;
3156 }
3157 IDA_mem->ida_tlo = IDA_mem->ida_trout;
3158 for (i = 0; i < IDA_mem->ida_nrtfn; i++)
3159 IDA_mem->ida_glo[i] = IDA_mem->ida_grout[i];
3160
3161 /* If no root found, return IDA_SUCCESS. */
3162 if (ier == IDA_SUCCESS) return(IDA_SUCCESS);
3163
3164 /* If a root was found, interpolate to get y(trout) and return. */
3165 (void) IDAGetSolution(IDA_mem, IDA_mem->ida_trout, IDA_mem->ida_yy, IDA_mem->ida_yp);
3166 return(RTFOUND);
3167 }
3168
3169 /*
3170 * IDARootfind
3171 *
3172 * This routine solves for a root of g(t) between tlo and thi, if
3173 * one exists. Only roots of odd multiplicity (i.e. with a change
3174 * of sign in one of the g_i), or exact zeros, are found.
3175 * Here the sign of tlo - thi is arbitrary, but if multiple roots
3176 * are found, the one closest to tlo is returned.
3177 *
3178 * The method used is the Illinois algorithm, a modified secant method.
3179 * Reference: Kathie L. Hiebert and Lawrence F. Shampine, Implicitly
3180 * Defined Output Points for Solutions of ODEs, Sandia National
3181 * Laboratory Report SAND80-0180, February 1980.
3182 *
3183 * This routine uses the following parameters for communication:
3184 *
3185 * nrtfn = number of functions g_i, or number of components of
3186 * the vector-valued function g(t). Input only.
3187 *
3188 * gfun = user-defined function for g(t). Its form is
3189 * (void) gfun(t, y, yp, gt, user_data)
3190 *
3191 * rootdir = in array specifying the direction of zero-crossings.
3192 * If rootdir[i] > 0, search for roots of g_i only if
3193 * g_i is increasing; if rootdir[i] < 0, search for
3194 * roots of g_i only if g_i is decreasing; otherwise
3195 * always search for roots of g_i.
3196 *
3197 * gactive = array specifying whether a component of g should
3198 * or should not be monitored. gactive[i] is initially
3199 * set to SUNTRUE for all i=0,...,nrtfn-1, but it may be
3200 * reset to SUNFALSE if at the first step g[i] is 0.0
3201 * both at the I.C. and at a small perturbation of them.
3202 * gactive[i] is then set back on SUNTRUE only after the
3203 * corresponding g function moves away from 0.0.
3204 *
3205 * nge = cumulative counter for gfun calls.
3206 *
3207 * ttol = a convergence tolerance for trout. Input only.
3208 * When a root at trout is found, it is located only to
3209 * within a tolerance of ttol. Typically, ttol should
3210 * be set to a value on the order of
3211 * 100 * UROUND * max (SUNRabs(tlo), SUNRabs(thi))
3212 * where UROUND is the unit roundoff of the machine.
3213 *
3214 * tlo, thi = endpoints of the interval in which roots are sought.
3215 * On input, these must be distinct, but tlo - thi may
3216 * be of either sign. The direction of integration is
3217 * assumed to be from tlo to thi. On return, tlo and thi
3218 * are the endpoints of the final relevant interval.
3219 *
3220 * glo, ghi = arrays of length nrtfn containing the vectors g(tlo)
3221 * and g(thi) respectively. Input and output. On input,
3222 * none of the glo[i] should be zero.
3223 *
3224 * trout = root location, if a root was found, or thi if not.
3225 * Output only. If a root was found other than an exact
3226 * zero of g, trout is the endpoint thi of the final
3227 * interval bracketing the root, with size at most ttol.
3228 *
3229 * grout = array of length nrtfn containing g(trout) on return.
3230 *
3231 * iroots = int array of length nrtfn with root information.
3232 * Output only. If a root was found, iroots indicates
3233 * which components g_i have a root at trout. For
3234 * i = 0, ..., nrtfn-1, iroots[i] = 1 if g_i has a root
3235 * and g_i is increasing, iroots[i] = -1 if g_i has a
3236 * root and g_i is decreasing, and iroots[i] = 0 if g_i
3237 * has no roots or g_i varies in the direction opposite
3238 * to that indicated by rootdir[i].
3239 *
3240 * This routine returns an int equal to:
3241 * IDA_RTFUNC_FAIL < 0 if the g function failed, or
3242 * RTFOUND = 1 if a root of g was found, or
3243 * IDA_SUCCESS = 0 otherwise.
3244 *
3245 */
3246
IDARootfind(IDAMem IDA_mem)3247 static int IDARootfind(IDAMem IDA_mem)
3248 {
3249 realtype alph, tmid, gfrac, maxfrac, fracint, fracsub;
3250 int i, retval, imax, side, sideprev;
3251 booleantype zroot, sgnchg;
3252
3253 imax = 0;
3254
3255 /* First check for change in sign in ghi or for a zero in ghi. */
3256 maxfrac = ZERO;
3257 zroot = SUNFALSE;
3258 sgnchg = SUNFALSE;
3259 for (i = 0; i < IDA_mem->ida_nrtfn; i++) {
3260 if(!IDA_mem->ida_gactive[i]) continue;
3261 if (SUNRabs(IDA_mem->ida_ghi[i]) == ZERO) {
3262 if(IDA_mem->ida_rootdir[i] * IDA_mem->ida_glo[i] <= ZERO) {
3263 zroot = SUNTRUE;
3264 }
3265 } else {
3266 if ( (IDA_mem->ida_glo[i] * IDA_mem->ida_ghi[i] < ZERO) &&
3267 (IDA_mem->ida_rootdir[i] * IDA_mem->ida_glo[i] <= ZERO) ) {
3268 gfrac = SUNRabs(IDA_mem->ida_ghi[i] / (IDA_mem->ida_ghi[i] - IDA_mem->ida_glo[i]));
3269 if (gfrac > maxfrac) {
3270 sgnchg = SUNTRUE;
3271 maxfrac = gfrac;
3272 imax = i;
3273 }
3274 }
3275 }
3276 }
3277
3278 /* If no sign change was found, reset trout and grout. Then return
3279 IDA_SUCCESS if no zero was found, or set iroots and return RTFOUND. */
3280 if (!sgnchg) {
3281 IDA_mem->ida_trout = IDA_mem->ida_thi;
3282 for (i = 0; i < IDA_mem->ida_nrtfn; i++)
3283 IDA_mem->ida_grout[i] = IDA_mem->ida_ghi[i];
3284 if (!zroot) return(IDA_SUCCESS);
3285 for (i = 0; i < IDA_mem->ida_nrtfn; i++) {
3286 IDA_mem->ida_iroots[i] = 0;
3287 if(!IDA_mem->ida_gactive[i]) continue;
3288 if ( (SUNRabs(IDA_mem->ida_ghi[i]) == ZERO) &&
3289 (IDA_mem->ida_rootdir[i] * IDA_mem->ida_glo[i] <= ZERO) )
3290 IDA_mem->ida_iroots[i] = IDA_mem->ida_glo[i] > 0 ? -1:1;
3291 }
3292 return(RTFOUND);
3293 }
3294
3295 /* Initialize alph to avoid compiler warning */
3296 alph = ONE;
3297
3298 /* A sign change was found. Loop to locate nearest root. */
3299
3300 side = 0; sideprev = -1;
3301 for(;;) { /* Looping point */
3302
3303 /* If interval size is already less than tolerance ttol, break. */
3304 if (SUNRabs(IDA_mem->ida_thi - IDA_mem->ida_tlo) <= IDA_mem->ida_ttol)
3305 break;
3306
3307 /* Set weight alph.
3308 On the first two passes, set alph = 1. Thereafter, reset alph
3309 according to the side (low vs high) of the subinterval in which
3310 the sign change was found in the previous two passes.
3311 If the sides were opposite, set alph = 1.
3312 If the sides were the same, then double alph (if high side),
3313 or halve alph (if low side).
3314 The next guess tmid is the secant method value if alph = 1, but
3315 is closer to tlo if alph < 1, and closer to thi if alph > 1. */
3316
3317 if (sideprev == side) {
3318 alph = (side == 2) ? alph*TWO : alph*HALF;
3319 } else {
3320 alph = ONE;
3321 }
3322
3323 /* Set next root approximation tmid and get g(tmid).
3324 If tmid is too close to tlo or thi, adjust it inward,
3325 by a fractional distance that is between 0.1 and 0.5. */
3326 tmid = IDA_mem->ida_thi - (IDA_mem->ida_thi - IDA_mem->ida_tlo) *
3327 IDA_mem->ida_ghi[imax] / (IDA_mem->ida_ghi[imax] - alph*IDA_mem->ida_glo[imax]);
3328 if (SUNRabs(tmid - IDA_mem->ida_tlo) < HALF * IDA_mem->ida_ttol) {
3329 fracint = SUNRabs(IDA_mem->ida_thi - IDA_mem->ida_tlo) / IDA_mem->ida_ttol;
3330 fracsub = (fracint > FIVE) ? PT1 : HALF/fracint;
3331 tmid = IDA_mem->ida_tlo + fracsub*(IDA_mem->ida_thi - IDA_mem->ida_tlo);
3332 }
3333 if (SUNRabs(IDA_mem->ida_thi - tmid) < HALF * IDA_mem->ida_ttol) {
3334 fracint = SUNRabs(IDA_mem->ida_thi - IDA_mem->ida_tlo) / IDA_mem->ida_ttol;
3335 fracsub = (fracint > FIVE) ? PT1 : HALF/fracint;
3336 tmid = IDA_mem->ida_thi - fracsub*(IDA_mem->ida_thi - IDA_mem->ida_tlo);
3337 }
3338
3339 (void) IDAGetSolution(IDA_mem, tmid, IDA_mem->ida_yy, IDA_mem->ida_yp);
3340 retval = IDA_mem->ida_gfun(tmid, IDA_mem->ida_yy, IDA_mem->ida_yp,
3341 IDA_mem->ida_grout, IDA_mem->ida_user_data);
3342 IDA_mem->ida_nge++;
3343 if (retval != 0) return(IDA_RTFUNC_FAIL);
3344
3345 /* Check to see in which subinterval g changes sign, and reset imax.
3346 Set side = 1 if sign change is on low side, or 2 if on high side. */
3347 maxfrac = ZERO;
3348 zroot = SUNFALSE;
3349 sgnchg = SUNFALSE;
3350 sideprev = side;
3351 for (i = 0; i < IDA_mem->ida_nrtfn; i++) {
3352 if(!IDA_mem->ida_gactive[i]) continue;
3353 if (SUNRabs(IDA_mem->ida_grout[i]) == ZERO) {
3354 if(IDA_mem->ida_rootdir[i] * IDA_mem->ida_glo[i] <= ZERO)
3355 zroot = SUNTRUE;
3356 } else {
3357 if ( (IDA_mem->ida_glo[i] * IDA_mem->ida_grout[i] < ZERO) &&
3358 (IDA_mem->ida_rootdir[i] * IDA_mem->ida_glo[i] <= ZERO) ) {
3359 gfrac = SUNRabs(IDA_mem->ida_grout[i] /
3360 (IDA_mem->ida_grout[i] - IDA_mem->ida_glo[i]));
3361 if (gfrac > maxfrac) {
3362 sgnchg = SUNTRUE;
3363 maxfrac = gfrac;
3364 imax = i;
3365 }
3366 }
3367 }
3368 }
3369 if (sgnchg) {
3370 /* Sign change found in (tlo,tmid); replace thi with tmid. */
3371 IDA_mem->ida_thi = tmid;
3372 for (i = 0; i < IDA_mem->ida_nrtfn; i++)
3373 IDA_mem->ida_ghi[i] = IDA_mem->ida_grout[i];
3374 side = 1;
3375 /* Stop at root thi if converged; otherwise loop. */
3376 if (SUNRabs(IDA_mem->ida_thi - IDA_mem->ida_tlo) <= IDA_mem->ida_ttol)
3377 break;
3378 continue; /* Return to looping point. */
3379 }
3380
3381 if (zroot) {
3382 /* No sign change in (tlo,tmid), but g = 0 at tmid; return root tmid. */
3383 IDA_mem->ida_thi = tmid;
3384 for (i = 0; i < IDA_mem->ida_nrtfn; i++)
3385 IDA_mem->ida_ghi[i] = IDA_mem->ida_grout[i];
3386 break;
3387 }
3388
3389 /* No sign change in (tlo,tmid), and no zero at tmid.
3390 Sign change must be in (tmid,thi). Replace tlo with tmid. */
3391 IDA_mem->ida_tlo = tmid;
3392 for (i = 0; i < IDA_mem->ida_nrtfn; i++)
3393 IDA_mem->ida_glo[i] = IDA_mem->ida_grout[i];
3394 side = 2;
3395 /* Stop at root thi if converged; otherwise loop back. */
3396 if (SUNRabs(IDA_mem->ida_thi - IDA_mem->ida_tlo) <= IDA_mem->ida_ttol)
3397 break;
3398
3399 } /* End of root-search loop */
3400
3401 /* Reset trout and grout, set iroots, and return RTFOUND. */
3402 IDA_mem->ida_trout = IDA_mem->ida_thi;
3403 for (i = 0; i < IDA_mem->ida_nrtfn; i++) {
3404 IDA_mem->ida_grout[i] = IDA_mem->ida_ghi[i];
3405 IDA_mem->ida_iroots[i] = 0;
3406 if(!IDA_mem->ida_gactive[i]) continue;
3407 if ( (SUNRabs(IDA_mem->ida_ghi[i]) == ZERO) &&
3408 (IDA_mem->ida_rootdir[i] * IDA_mem->ida_glo[i] <= ZERO) )
3409 IDA_mem->ida_iroots[i] = IDA_mem->ida_glo[i] > 0 ? -1:1;
3410 if ( (IDA_mem->ida_glo[i] * IDA_mem->ida_ghi[i] < ZERO) &&
3411 (IDA_mem->ida_rootdir[i] * IDA_mem->ida_glo[i] <= ZERO) )
3412 IDA_mem->ida_iroots[i] = IDA_mem->ida_glo[i] > 0 ? -1:1;
3413 }
3414 return(RTFOUND);
3415 }
3416
3417 /*
3418 * =================================================================
3419 * IDA error message handling functions
3420 * =================================================================
3421 */
3422
3423 /*
3424 * IDAProcessError is a high level error handling function.
3425 * - If ida_mem==NULL it prints the error message to stderr.
3426 * - Otherwise, it sets up and calls the error handling function
3427 * pointed to by ida_ehfun.
3428 */
3429
IDAProcessError(IDAMem IDA_mem,int error_code,const char * module,const char * fname,const char * msgfmt,...)3430 void IDAProcessError(IDAMem IDA_mem,
3431 int error_code, const char *module, const char *fname,
3432 const char *msgfmt, ...)
3433 {
3434 va_list ap;
3435 char msg[256];
3436
3437 /* Initialize the argument pointer variable
3438 (msgfmt is the last required argument to IDAProcessError) */
3439
3440 va_start(ap, msgfmt);
3441
3442 /* Compose the message */
3443
3444 vsprintf(msg, msgfmt, ap);
3445
3446 if (IDA_mem == NULL) { /* We write to stderr */
3447 #ifndef NO_FPRINTF_OUTPUT
3448 fprintf(stderr, "\n[%s ERROR] %s\n ", module, fname);
3449 fprintf(stderr, "%s\n\n", msg);
3450 #endif
3451
3452 } else { /* We can call ehfun */
3453 IDA_mem->ida_ehfun(error_code, module, fname, msg, IDA_mem->ida_eh_data);
3454 }
3455
3456 /* Finalize argument processing */
3457 va_end(ap);
3458
3459 return;
3460 }
3461
3462 /* IDAErrHandler is the default error handling function.
3463 It sends the error message to the stream pointed to by ida_errfp */
3464
IDAErrHandler(int error_code,const char * module,const char * function,char * msg,void * data)3465 void IDAErrHandler(int error_code, const char *module,
3466 const char *function, char *msg, void *data)
3467 {
3468 IDAMem IDA_mem;
3469 char err_type[10];
3470
3471 /* data points to IDA_mem here */
3472
3473 IDA_mem = (IDAMem) data;
3474
3475 if (error_code == IDA_WARNING)
3476 sprintf(err_type,"WARNING");
3477 else
3478 sprintf(err_type,"ERROR");
3479
3480 #ifndef NO_FPRINTF_OUTPUT
3481 if (IDA_mem->ida_errfp != NULL) {
3482 fprintf(IDA_mem->ida_errfp,"\n[%s %s] %s\n",module,err_type,function);
3483 fprintf(IDA_mem->ida_errfp," %s\n\n",msg);
3484 }
3485 #endif
3486
3487 return;
3488 }
3489