1 /*-----------------------------------------------------------------
2 * Programmer(s): Daniel R. Reynolds @ SMU
3 * Alan C. Hindmarsh and Radu Serban @ LLNL
4 *-----------------------------------------------------------------
5 * SUNDIALS Copyright Start
6 * Copyright (c) 2002-2021, Lawrence Livermore National Security
7 * and Southern Methodist University.
8 * All rights reserved.
9 *
10 * See the top-level LICENSE and NOTICE files for details.
11 *
12 * SPDX-License-Identifier: BSD-3-Clause
13 * SUNDIALS Copyright End
14 *-----------------------------------------------------------------
15 * Implementation file for IDAS' linear solver interface
16 *-----------------------------------------------------------------*/
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <string.h>
21
22 #include "idas_impl.h"
23 #include "idas_ls_impl.h"
24 #include <sundials/sundials_math.h>
25 #include <sundials/sundials_linearsolver.h>
26 #include <sunmatrix/sunmatrix_band.h>
27 #include <sunmatrix/sunmatrix_dense.h>
28 #include <sunmatrix/sunmatrix_sparse.h>
29
30 /* constants */
31 #define MAX_ITERS 3 /* max. number of attempts to recover in DQ J*v */
32 #define ZERO RCONST(0.0)
33 #define PT25 RCONST(0.25)
34 #define PT05 RCONST(0.05)
35 #define PT9 RCONST(0.9)
36 #define ONE RCONST(1.0)
37 #define TWO RCONST(2.0)
38
39
40 /*=================================================================
41 PRIVATE FUNCTION PROTOTYPES
42 =================================================================*/
43
44 static int idaLsJacBWrapper(realtype tt, realtype c_jB, N_Vector yyB,
45 N_Vector ypB, N_Vector rBr, SUNMatrix JacB,
46 void *ida_mem, N_Vector tmp1B,
47 N_Vector tmp2B, N_Vector tmp3B);
48 static int idaLsJacBSWrapper(realtype tt, realtype c_jB, N_Vector yyB,
49 N_Vector ypB, N_Vector rBr, SUNMatrix JacB,
50 void *ida_mem, N_Vector tmp1B,
51 N_Vector tmp2B, N_Vector tmp3B);
52
53 static int idaLsPrecSetupB(realtype tt, N_Vector yyB,
54 N_Vector ypB, N_Vector rrB,
55 realtype c_jB, void *idaadj_mem);
56 static int idaLsPrecSetupBS(realtype tt, N_Vector yyB,
57 N_Vector ypB, N_Vector rrB,
58 realtype c_jB, void *idaadj_mem);
59
60 static int idaLsPrecSolveB(realtype tt, N_Vector yyB,
61 N_Vector ypB, N_Vector rrB,
62 N_Vector rvecB, N_Vector zvecB,
63 realtype c_jB, realtype deltaB,
64 void *idaadj_mem);
65 static int idaLsPrecSolveBS(realtype tt, N_Vector yyB,
66 N_Vector ypB, N_Vector rrB,
67 N_Vector rvecB, N_Vector zvecB,
68 realtype c_jB, realtype deltaB,
69 void *idaadj_mem);
70
71 static int idaLsJacTimesSetupB(realtype tt, N_Vector yyB,
72 N_Vector ypB, N_Vector rrB,
73 realtype c_jB, void *idaadj_mem);
74 static int idaLsJacTimesSetupBS(realtype tt, N_Vector yyB,
75 N_Vector ypB, N_Vector rrB,
76 realtype c_jB, void *idaadj_mem);
77
78 static int idaLsJacTimesVecB(realtype tt, N_Vector yyB,
79 N_Vector ypB, N_Vector rrB,
80 N_Vector vB, N_Vector JvB,
81 realtype c_jB, void *idaadj_mem,
82 N_Vector tmp1B, N_Vector tmp2B);
83 static int idaLsJacTimesVecBS(realtype tt, N_Vector yyB,
84 N_Vector ypB, N_Vector rrB,
85 N_Vector vB, N_Vector JvB,
86 realtype c_jB, void *idaadj_mem,
87 N_Vector tmp1B, N_Vector tmp2B);
88
89
90 /*================================================================
91 PART I - forward problems
92 ================================================================*/
93
94
95 /*---------------------------------------------------------------
96 IDASetLinearSolver specifies the linear solver
97 ---------------------------------------------------------------*/
IDASetLinearSolver(void * ida_mem,SUNLinearSolver LS,SUNMatrix A)98 int IDASetLinearSolver(void *ida_mem, SUNLinearSolver LS, SUNMatrix A)
99 {
100 IDAMem IDA_mem;
101 IDALsMem idals_mem;
102 int retval, LSType;
103 booleantype iterative; /* is the solver iterative? */
104 booleantype matrixbased; /* is a matrix structure used? */
105
106 /* Return immediately if any input is NULL */
107 if (ida_mem == NULL) {
108 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASLS",
109 "IDASetLinearSolver", MSG_LS_IDAMEM_NULL);
110 return(IDALS_MEM_NULL);
111 }
112 if (LS == NULL) {
113 IDAProcessError(NULL, IDALS_ILL_INPUT, "IDASLS",
114 "IDASetLinearSolver",
115 "LS must be non-NULL");
116 return(IDALS_ILL_INPUT);
117 }
118 IDA_mem = (IDAMem) ida_mem;
119
120 /* Test if solver is compatible with LS interface */
121 if ( (LS->ops->gettype == NULL) || (LS->ops->solve == NULL) ) {
122 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS",
123 "IDASetLinearSolver",
124 "LS object is missing a required operation");
125 return(IDALS_ILL_INPUT);
126 }
127
128 /* Retrieve the LS type */
129 LSType = SUNLinSolGetType(LS);
130
131 /* Set flags based on LS type */
132 iterative = (LSType != SUNLINEARSOLVER_DIRECT);
133 matrixbased = (LSType != SUNLINEARSOLVER_ITERATIVE);
134
135 /* Test if vector is compatible with LS interface */
136 if (IDA_mem->ida_tempv1->ops->nvconst == NULL ||
137 IDA_mem->ida_tempv1->ops->nvwrmsnorm == NULL) {
138 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS",
139 "IDASetLinearSolver", MSG_LS_BAD_NVECTOR);
140 return(IDALS_ILL_INPUT);
141 }
142
143 /* Check for compatible LS type, matrix and "atimes" support */
144 if (iterative) {
145
146 if (IDA_mem->ida_tempv1->ops->nvgetlength == NULL) {
147 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS",
148 "IDASetLinearSolver", MSG_LS_BAD_NVECTOR);
149 return(IDALS_ILL_INPUT);
150 }
151
152 if (LS->ops->resid == NULL || LS->ops->numiters == NULL) {
153 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS", "IDASetLinearSolver",
154 "Iterative LS object requires 'resid' and 'numiters' routines");
155 return(IDALS_ILL_INPUT);
156 }
157
158 if (!matrixbased && LS->ops->setatimes == NULL) {
159 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS", "IDASetLinearSolver",
160 "Incompatible inputs: iterative LS must support ATimes routine");
161 return(IDALS_ILL_INPUT);
162 }
163
164 if (matrixbased && A == NULL) {
165 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS", "IDASetLinearSolver",
166 "Incompatible inputs: matrix-iterative LS requires non-NULL matrix");
167 return(IDALS_ILL_INPUT);
168 }
169
170 } else if (A == NULL) {
171
172 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS", "IDASetLinearSolver",
173 "Incompatible inputs: direct LS requires non-NULL matrix");
174 return(IDALS_ILL_INPUT);
175
176 }
177
178 /* free any existing system solver attached to IDA */
179 if (IDA_mem->ida_lfree) IDA_mem->ida_lfree(IDA_mem);
180
181 /* Set four main system linear solver function fields in IDA_mem */
182 IDA_mem->ida_linit = idaLsInitialize;
183 IDA_mem->ida_lsetup = idaLsSetup;
184 IDA_mem->ida_lsolve = idaLsSolve;
185 IDA_mem->ida_lfree = idaLsFree;
186
187 /* Set ida_lperf if using an iterative SUNLinearSolver object */
188 IDA_mem->ida_lperf = (iterative) ? idaLsPerf : NULL;
189
190 /* Allocate memory for IDALsMemRec */
191 idals_mem = NULL;
192 idals_mem = (IDALsMem) malloc(sizeof(struct IDALsMemRec));
193 if (idals_mem == NULL) {
194 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASLS",
195 "IDASetLinearSolver", MSG_LS_MEM_FAIL);
196 return(IDALS_MEM_FAIL);
197 }
198 memset(idals_mem, 0, sizeof(struct IDALsMemRec));
199
200 /* set SUNLinearSolver pointer */
201 idals_mem->LS = LS;
202
203 /* Linear solver type information */
204 idals_mem->iterative = iterative;
205 idals_mem->matrixbased = matrixbased;
206
207 /* Set defaults for Jacobian-related fields */
208 idals_mem->J = A;
209 if (A != NULL) {
210 idals_mem->jacDQ = SUNTRUE;
211 idals_mem->jac = idaLsDQJac;
212 idals_mem->J_data = IDA_mem;
213 } else {
214 idals_mem->jacDQ = SUNFALSE;
215 idals_mem->jac = NULL;
216 idals_mem->J_data = NULL;
217 }
218 idals_mem->jtimesDQ = SUNTRUE;
219 idals_mem->jtsetup = NULL;
220 idals_mem->jtimes = idaLsDQJtimes;
221 idals_mem->jt_res = IDA_mem->ida_res;
222 idals_mem->jt_data = IDA_mem;
223
224 /* Set defaults for preconditioner-related fields */
225 idals_mem->pset = NULL;
226 idals_mem->psolve = NULL;
227 idals_mem->pfree = NULL;
228 idals_mem->pdata = IDA_mem->ida_user_data;
229
230 /* Initialize counters */
231 idaLsInitializeCounters(idals_mem);
232
233 /* Set default values for the rest of the Ls parameters */
234 idals_mem->eplifac = PT05;
235 idals_mem->dqincfac = ONE;
236 idals_mem->last_flag = IDALS_SUCCESS;
237
238 /* If LS supports ATimes, attach IDALs routine */
239 if (LS->ops->setatimes) {
240 retval = SUNLinSolSetATimes(LS, IDA_mem, idaLsATimes);
241 if (retval != SUNLS_SUCCESS) {
242 IDAProcessError(IDA_mem, IDALS_SUNLS_FAIL, "IDASLS",
243 "IDASetLinearSolver",
244 "Error in calling SUNLinSolSetATimes");
245 free(idals_mem); idals_mem = NULL;
246 return(IDALS_SUNLS_FAIL);
247 }
248 }
249
250 /* If LS supports preconditioning, initialize pset/psol to NULL */
251 if (LS->ops->setpreconditioner) {
252 retval = SUNLinSolSetPreconditioner(LS, IDA_mem, NULL, NULL);
253 if (retval != SUNLS_SUCCESS) {
254 IDAProcessError(IDA_mem, IDALS_SUNLS_FAIL, "IDASLS",
255 "IDASetLinearSolver",
256 "Error in calling SUNLinSolSetPreconditioner");
257 free(idals_mem); idals_mem = NULL;
258 return(IDALS_SUNLS_FAIL);
259 }
260 }
261
262 /* Allocate memory for ytemp, yptemp and x */
263 idals_mem->ytemp = N_VClone(IDA_mem->ida_tempv1);
264 if (idals_mem->ytemp == NULL) {
265 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASLS",
266 "IDASetLinearSolver", MSG_LS_MEM_FAIL);
267 free(idals_mem); idals_mem = NULL;
268 return(IDALS_MEM_FAIL);
269 }
270
271 idals_mem->yptemp = N_VClone(IDA_mem->ida_tempv1);
272 if (idals_mem->yptemp == NULL) {
273 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASLS",
274 "IDASetLinearSolver", MSG_LS_MEM_FAIL);
275 N_VDestroy(idals_mem->ytemp);
276 free(idals_mem); idals_mem = NULL;
277 return(IDALS_MEM_FAIL);
278 }
279
280 idals_mem->x = N_VClone(IDA_mem->ida_tempv1);
281 if (idals_mem->x == NULL) {
282 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASLS",
283 "IDASetLinearSolver", MSG_LS_MEM_FAIL);
284 N_VDestroy(idals_mem->ytemp);
285 N_VDestroy(idals_mem->yptemp);
286 free(idals_mem); idals_mem = NULL;
287 return(IDALS_MEM_FAIL);
288 }
289
290 /* For iterative LS, compute sqrtN */
291 if (iterative)
292 idals_mem->nrmfac = SUNRsqrt( N_VGetLength(idals_mem->ytemp) );
293
294 /* For matrix-based LS, enable soltuion scaling */
295 if (matrixbased)
296 idals_mem->scalesol = SUNTRUE;
297 else
298 idals_mem->scalesol = SUNFALSE;
299
300 /* Attach linear solver memory to integrator memory */
301 IDA_mem->ida_lmem = idals_mem;
302
303 return(IDALS_SUCCESS);
304 }
305
306
307 /*===============================================================
308 Optional input/output routines
309 ===============================================================*/
310
311
312 /* IDASetJacFn specifies the Jacobian function */
IDASetJacFn(void * ida_mem,IDALsJacFn jac)313 int IDASetJacFn(void *ida_mem, IDALsJacFn jac)
314 {
315 IDAMem IDA_mem;
316 IDALsMem idals_mem;
317 int retval;
318
319 /* access IDALsMem structure */
320 retval = idaLs_AccessLMem(ida_mem, "IDALsSetJacFn",
321 &IDA_mem, &idals_mem);
322 if (retval != IDALS_SUCCESS) return(retval);
323
324 /* return with failure if jac cannot be used */
325 if ((jac != NULL) && (idals_mem->J == NULL)) {
326 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS", "IDASetJacFn",
327 "Jacobian routine cannot be supplied for NULL SUNMatrix");
328 return(IDALS_ILL_INPUT);
329 }
330
331 /* set Jacobian routine pointer, and update relevant flags */
332 if (jac != NULL) {
333 idals_mem->jacDQ = SUNFALSE;
334 idals_mem->jac = jac;
335 idals_mem->J_data = IDA_mem->ida_user_data;
336 } else {
337 idals_mem->jacDQ = SUNTRUE;
338 idals_mem->jac = idaLsDQJac;
339 idals_mem->J_data = IDA_mem;
340 }
341
342 return(IDALS_SUCCESS);
343 }
344
345
346 /* IDASetEpsLin specifies the nonlinear -> linear tolerance scale factor */
IDASetEpsLin(void * ida_mem,realtype eplifac)347 int IDASetEpsLin(void *ida_mem, realtype eplifac)
348 {
349 IDAMem IDA_mem;
350 IDALsMem idals_mem;
351 int retval;
352
353 /* access IDALsMem structure */
354 retval = idaLs_AccessLMem(ida_mem, "IDASetEpsLin",
355 &IDA_mem, &idals_mem);
356 if (retval != IDALS_SUCCESS) return(retval);
357
358 /* Check for legal eplifac */
359 if (eplifac < ZERO) {
360 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS",
361 "IDASetEpsLin", MSG_LS_NEG_EPLIFAC);
362 return(IDALS_ILL_INPUT);
363 }
364
365 idals_mem->eplifac = (eplifac == ZERO) ? PT05 : eplifac;
366
367 return(IDALS_SUCCESS);
368 }
369
370
371 /* IDASetWRMSNormFactor sets or computes the factor to use when converting from
372 the integrator tolerance to the linear solver tolerance (WRMS to L2 norm). */
IDASetLSNormFactor(void * ida_mem,realtype nrmfac)373 int IDASetLSNormFactor(void *ida_mem, realtype nrmfac)
374 {
375 IDAMem IDA_mem;
376 IDALsMem idals_mem;
377 int retval;
378
379 /* access IDALsMem structure */
380 retval = idaLs_AccessLMem(ida_mem, "IDASetLSNormFactor",
381 &IDA_mem, &idals_mem);
382 if (retval != IDALS_SUCCESS) return(retval);
383
384 if (nrmfac > ZERO) {
385 /* user-provided factor */
386 idals_mem->nrmfac = nrmfac;
387 } else if (nrmfac < ZERO) {
388 /* compute factor for WRMS norm with dot product */
389 N_VConst(ONE, idals_mem->ytemp);
390 idals_mem->nrmfac = SUNRsqrt(N_VDotProd(idals_mem->ytemp, idals_mem->ytemp));
391 } else {
392 /* compute default factor for WRMS norm from vector legnth */
393 idals_mem->nrmfac = SUNRsqrt(N_VGetLength(idals_mem->ytemp));
394 }
395
396 return(IDALS_SUCCESS);
397 }
398
399
400 /* IDASetLinearSolutionScaling enables or disables scaling the linear solver
401 solution to account for changes in cj. */
IDASetLinearSolutionScaling(void * ida_mem,booleantype onoff)402 int IDASetLinearSolutionScaling(void *ida_mem, booleantype onoff)
403 {
404 IDAMem IDA_mem;
405 IDALsMem idals_mem;
406 int retval;
407
408 /* access IDALsMem structure */
409 retval = idaLs_AccessLMem(ida_mem, "IDASetLinearSolutionScaling",
410 &IDA_mem, &idals_mem);
411 if (retval != IDALS_SUCCESS) return(retval);
412
413 /* check for valid solver type */
414 if (!(idals_mem->matrixbased)) return(IDALS_ILL_INPUT);
415
416 /* set solution scaling flag */
417 idals_mem->scalesol = onoff;
418
419 return(IDALS_SUCCESS);
420 }
421
422
423 /* IDASetIncrementFactor specifies increment factor for DQ approximations to Jv */
IDASetIncrementFactor(void * ida_mem,realtype dqincfac)424 int IDASetIncrementFactor(void *ida_mem, realtype dqincfac)
425 {
426 IDAMem IDA_mem;
427 IDALsMem idals_mem;
428 int retval;
429
430 /* access IDALsMem structure */
431 retval = idaLs_AccessLMem(ida_mem, "IDASetIncrementFactor",
432 &IDA_mem, &idals_mem);
433 if (retval != IDALS_SUCCESS) return(retval);
434
435 /* Check for legal dqincfac */
436 if (dqincfac <= ZERO) {
437 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS",
438 "IDASetIncrementFactor", MSG_LS_NEG_DQINCFAC);
439 return(IDALS_ILL_INPUT);
440 }
441
442 idals_mem->dqincfac = dqincfac;
443
444 return(IDALS_SUCCESS);
445 }
446
447
448 /* IDASetPreconditioner specifies the user-supplied psetup and psolve routines */
IDASetPreconditioner(void * ida_mem,IDALsPrecSetupFn psetup,IDALsPrecSolveFn psolve)449 int IDASetPreconditioner(void *ida_mem,
450 IDALsPrecSetupFn psetup,
451 IDALsPrecSolveFn psolve)
452 {
453 IDAMem IDA_mem;
454 IDALsMem idals_mem;
455 PSetupFn idals_psetup;
456 PSolveFn idals_psolve;
457 int retval;
458
459 /* access IDALsMem structure */
460 retval = idaLs_AccessLMem(ida_mem, "IDASetPreconditioner",
461 &IDA_mem, &idals_mem);
462 if (retval != IDALS_SUCCESS) return(retval);
463
464 /* store function pointers for user-supplied routines in IDALs interface */
465 idals_mem->pset = psetup;
466 idals_mem->psolve = psolve;
467
468 /* issue error if LS object does not allow user-supplied preconditioning */
469 if (idals_mem->LS->ops->setpreconditioner == NULL) {
470 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS",
471 "IDASetPreconditioner",
472 "SUNLinearSolver object does not support user-supplied preconditioning");
473 return(IDALS_ILL_INPUT);
474 }
475
476 /* notify iterative linear solver to call IDALs interface routines */
477 idals_psetup = (psetup == NULL) ? NULL : idaLsPSetup;
478 idals_psolve = (psolve == NULL) ? NULL : idaLsPSolve;
479 retval = SUNLinSolSetPreconditioner(idals_mem->LS, IDA_mem,
480 idals_psetup, idals_psolve);
481 if (retval != SUNLS_SUCCESS) {
482 IDAProcessError(IDA_mem, IDALS_SUNLS_FAIL, "IDASLS",
483 "IDASetPreconditioner",
484 "Error in calling SUNLinSolSetPreconditioner");
485 return(IDALS_SUNLS_FAIL);
486 }
487
488 return(IDALS_SUCCESS);
489 }
490
491
492 /* IDASetJacTimes specifies the user-supplied Jacobian-vector product
493 setup and multiply routines */
IDASetJacTimes(void * ida_mem,IDALsJacTimesSetupFn jtsetup,IDALsJacTimesVecFn jtimes)494 int IDASetJacTimes(void *ida_mem, IDALsJacTimesSetupFn jtsetup,
495 IDALsJacTimesVecFn jtimes)
496 {
497 IDAMem IDA_mem;
498 IDALsMem idals_mem;
499 int retval;
500
501 /* access IDALsMem structure */
502 retval = idaLs_AccessLMem(ida_mem, "IDASetJacTimes",
503 &IDA_mem, &idals_mem);
504 if (retval != IDALS_SUCCESS) return(retval);
505
506 /* issue error if LS object does not allow user-supplied ATimes */
507 if (idals_mem->LS->ops->setatimes == NULL) {
508 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS",
509 "IDASetJacTimes",
510 "SUNLinearSolver object does not support user-supplied ATimes routine");
511 return(IDALS_ILL_INPUT);
512 }
513
514 /* store function pointers for user-supplied routines in IDALs
515 interface (NULL jtimes implies use of DQ default) */
516 if (jtimes != NULL) {
517 idals_mem->jtimesDQ = SUNFALSE;
518 idals_mem->jtsetup = jtsetup;
519 idals_mem->jtimes = jtimes;
520 idals_mem->jt_data = IDA_mem->ida_user_data;
521 } else {
522 idals_mem->jtimesDQ = SUNTRUE;
523 idals_mem->jtsetup = NULL;
524 idals_mem->jtimes = idaLsDQJtimes;
525 idals_mem->jt_res = IDA_mem->ida_res;
526 idals_mem->jt_data = IDA_mem;
527 }
528
529 return(IDALS_SUCCESS);
530 }
531
532
533 /* IDASetJacTimesResFn specifies an alternative user-supplied DAE residual
534 function to use in the internal finite difference Jacobian-vector
535 product */
IDASetJacTimesResFn(void * ida_mem,IDAResFn jtimesResFn)536 int IDASetJacTimesResFn(void *ida_mem, IDAResFn jtimesResFn)
537 {
538 IDAMem IDA_mem;
539 IDALsMem idals_mem;
540 int retval;
541
542 /* access IDALsMem structure */
543 retval = idaLs_AccessLMem(ida_mem, "IDASetJacTimesResFn",
544 &IDA_mem, &idals_mem);
545 if (retval != IDALS_SUCCESS) return(retval);
546
547 /* check if using internal finite difference approximation */
548 if (!(idals_mem->jtimesDQ)) {
549 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS", "IDASetJacTimesResFn",
550 "Internal finite-difference Jacobian-vector product is disabled.");
551 return(IDALS_ILL_INPUT);
552 }
553
554 /* store function pointers for Res function (NULL implies use DAE Res) */
555 if (jtimesResFn != NULL)
556 idals_mem->jt_res = jtimesResFn;
557 else
558 idals_mem->jt_res = IDA_mem->ida_res;
559
560 return(IDALS_SUCCESS);
561 }
562
563
564 /* IDAGetLinWorkSpace returns the length of workspace allocated
565 for the IDALS linear solver interface */
IDAGetLinWorkSpace(void * ida_mem,long int * lenrwLS,long int * leniwLS)566 int IDAGetLinWorkSpace(void *ida_mem, long int *lenrwLS,
567 long int *leniwLS)
568 {
569 IDAMem IDA_mem;
570 IDALsMem idals_mem;
571 sunindextype lrw1, liw1;
572 long int lrw, liw;
573 int retval;
574
575 /* access IDALsMem structure */
576 retval = idaLs_AccessLMem(ida_mem, "IDAGetLinWorkSpace",
577 &IDA_mem, &idals_mem);
578 if (retval != IDALS_SUCCESS) return(retval);
579
580 /* start with fixed sizes plus vector/matrix pointers */
581 *lenrwLS = 3;
582 *leniwLS = 34;
583
584 /* add N_Vector sizes */
585 if (IDA_mem->ida_tempv1->ops->nvspace) {
586 N_VSpace(IDA_mem->ida_tempv1, &lrw1, &liw1);
587 *lenrwLS += 3*lrw1;
588 *leniwLS += 3*liw1;
589 }
590
591 /* add LS sizes */
592 if (idals_mem->LS->ops->space) {
593 retval = SUNLinSolSpace(idals_mem->LS, &lrw, &liw);
594 if (retval == 0) {
595 *lenrwLS += lrw;
596 *leniwLS += liw;
597 }
598 }
599
600 return(IDALS_SUCCESS);
601 }
602
603
604 /* IDAGetNumJacEvals returns the number of Jacobian evaluations */
IDAGetNumJacEvals(void * ida_mem,long int * njevals)605 int IDAGetNumJacEvals(void *ida_mem, long int *njevals)
606 {
607 IDAMem IDA_mem;
608 IDALsMem idals_mem;
609 int retval;
610
611 /* access IDALsMem structure; store output and return */
612 retval = idaLs_AccessLMem(ida_mem, "IDAGetNumJacEvals",
613 &IDA_mem, &idals_mem);
614 if (retval != IDALS_SUCCESS) return(retval);
615 *njevals = idals_mem->nje;
616 return(IDALS_SUCCESS);
617 }
618
619
620 /* IDAGetNumPrecEvals returns the number of preconditioner evaluations */
IDAGetNumPrecEvals(void * ida_mem,long int * npevals)621 int IDAGetNumPrecEvals(void *ida_mem, long int *npevals)
622 {
623 IDAMem IDA_mem;
624 IDALsMem idals_mem;
625 int retval;
626
627 /* access IDALsMem structure; store output and return */
628 retval = idaLs_AccessLMem(ida_mem, "IDAGetNumPrecEvals",
629 &IDA_mem, &idals_mem);
630 if (retval != IDALS_SUCCESS) return(retval);
631 *npevals = idals_mem->npe;
632 return(IDALS_SUCCESS);
633 }
634
635
636 /* IDAGetNumPrecSolves returns the number of preconditioner solves */
IDAGetNumPrecSolves(void * ida_mem,long int * npsolves)637 int IDAGetNumPrecSolves(void *ida_mem, long int *npsolves)
638 {
639 IDAMem IDA_mem;
640 IDALsMem idals_mem;
641 int retval;
642
643 /* access IDALsMem structure; store output and return */
644 retval = idaLs_AccessLMem(ida_mem, "IDAGetNumPrecSolves",
645 &IDA_mem, &idals_mem);
646 if (retval != IDALS_SUCCESS) return(retval);
647 *npsolves = idals_mem->nps;
648 return(IDALS_SUCCESS);
649 }
650
651
652 /* IDAGetNumLinIters returns the number of linear iterations */
IDAGetNumLinIters(void * ida_mem,long int * nliters)653 int IDAGetNumLinIters(void *ida_mem, long int *nliters)
654 {
655 IDAMem IDA_mem;
656 IDALsMem idals_mem;
657 int retval;
658
659 /* access IDALsMem structure; store output and return */
660 retval = idaLs_AccessLMem(ida_mem, "IDAGetNumLinIters",
661 &IDA_mem, &idals_mem);
662 if (retval != IDALS_SUCCESS) return(retval);
663 *nliters = idals_mem->nli;
664 return(IDALS_SUCCESS);
665 }
666
667
668 /* IDAGetNumLinConvFails returns the number of linear convergence failures */
IDAGetNumLinConvFails(void * ida_mem,long int * nlcfails)669 int IDAGetNumLinConvFails(void *ida_mem, long int *nlcfails)
670 {
671 IDAMem IDA_mem;
672 IDALsMem idals_mem;
673 int retval;
674
675 /* access IDALsMem structure; store output and return */
676 retval = idaLs_AccessLMem(ida_mem, "IDAGetNumLinConvFails",
677 &IDA_mem, &idals_mem);
678 if (retval != IDALS_SUCCESS) return(retval);
679 *nlcfails = idals_mem->ncfl;
680 return(IDALS_SUCCESS);
681 }
682
683
684 /* IDAGetNumJTSetupEvals returns the number of calls to the
685 user-supplied Jacobian-vector product setup routine */
IDAGetNumJTSetupEvals(void * ida_mem,long int * njtsetups)686 int IDAGetNumJTSetupEvals(void *ida_mem, long int *njtsetups)
687 {
688 IDAMem IDA_mem;
689 IDALsMem idals_mem;
690 int retval;
691
692 /* access IDALsMem structure; store output and return */
693 retval = idaLs_AccessLMem(ida_mem, "IDAGetNumJTSetupEvals",
694 &IDA_mem, &idals_mem);
695 if (retval != IDALS_SUCCESS) return(retval);
696 *njtsetups = idals_mem->njtsetup;
697 return(IDALS_SUCCESS);
698 }
699
700
701 /* IDAGetNumJtimesEvals returns the number of calls to the
702 Jacobian-vector product multiply routine */
IDAGetNumJtimesEvals(void * ida_mem,long int * njvevals)703 int IDAGetNumJtimesEvals(void *ida_mem, long int *njvevals)
704 {
705 IDAMem IDA_mem;
706 IDALsMem idals_mem;
707 int retval;
708
709 /* access IDALsMem structure; store output and return */
710 retval = idaLs_AccessLMem(ida_mem, "IDAGetNumJtimesEvals",
711 &IDA_mem, &idals_mem);
712 if (retval != IDALS_SUCCESS) return(retval);
713 *njvevals = idals_mem->njtimes;
714 return(IDALS_SUCCESS);
715 }
716
717
718 /* IDAGetNumLinResEvals returns the number of calls to the DAE
719 residual needed for the DQ Jacobian approximation or J*v
720 product approximation */
IDAGetNumLinResEvals(void * ida_mem,long int * nrevalsLS)721 int IDAGetNumLinResEvals(void *ida_mem, long int *nrevalsLS)
722 {
723 IDAMem IDA_mem;
724 IDALsMem idals_mem;
725 int retval;
726
727 /* access IDALsMem structure; store output and return */
728 retval = idaLs_AccessLMem(ida_mem, "IDAGetNumLinResEvals",
729 &IDA_mem, &idals_mem);
730 if (retval != IDALS_SUCCESS) return(retval);
731 *nrevalsLS = idals_mem->nreDQ;
732 return(IDALS_SUCCESS);
733 }
734
735
736 /* IDAGetLastLinFlag returns the last flag set in a IDALS function */
IDAGetLastLinFlag(void * ida_mem,long int * flag)737 int IDAGetLastLinFlag(void *ida_mem, long int *flag)
738 {
739 IDAMem IDA_mem;
740 IDALsMem idals_mem;
741 int retval;
742
743 /* access IDALsMem structure; store output and return */
744 retval = idaLs_AccessLMem(ida_mem, "IDAGetLastLinFlag",
745 &IDA_mem, &idals_mem);
746 if (retval != IDALS_SUCCESS) return(retval);
747 *flag = idals_mem->last_flag;
748 return(IDALS_SUCCESS);
749 }
750
751
752 /* IDAGetLinReturnFlagName translates from the integer error code
753 returned by an IDALs routine to the corresponding string
754 equivalent for that flag */
IDAGetLinReturnFlagName(long int flag)755 char *IDAGetLinReturnFlagName(long int flag)
756 {
757 char *name = (char *)malloc(30*sizeof(char));
758
759 switch(flag) {
760 case IDALS_SUCCESS:
761 sprintf(name,"IDALS_SUCCESS");
762 break;
763 case IDALS_MEM_NULL:
764 sprintf(name,"IDALS_MEM_NULL");
765 break;
766 case IDALS_LMEM_NULL:
767 sprintf(name,"IDALS_LMEM_NULL");
768 break;
769 case IDALS_ILL_INPUT:
770 sprintf(name,"IDALS_ILL_INPUT");
771 break;
772 case IDALS_MEM_FAIL:
773 sprintf(name,"IDALS_MEM_FAIL");
774 break;
775 case IDALS_PMEM_NULL:
776 sprintf(name,"IDALS_PMEM_NULL");
777 break;
778 case IDALS_JACFUNC_UNRECVR:
779 sprintf(name,"IDALS_JACFUNC_UNRECVR");
780 break;
781 case IDALS_JACFUNC_RECVR:
782 sprintf(name,"IDALS_JACFUNC_RECVR");
783 break;
784 case IDALS_SUNMAT_FAIL:
785 sprintf(name,"IDALS_SUNMAT_FAIL");
786 break;
787 case IDALS_SUNLS_FAIL:
788 sprintf(name,"IDALS_SUNLS_FAIL");
789 break;
790 default:
791 sprintf(name,"NONE");
792 }
793
794 return(name);
795 }
796
797
798 /*===============================================================
799 IDASLS Private functions
800 ===============================================================*/
801
802 /*---------------------------------------------------------------
803 idaLsATimes:
804
805 This routine generates the matrix-vector product z = Jv, where
806 J is the system Jacobian, by calling either the user provided
807 routine or the internal DQ routine. The return value is
808 the same as the value returned by jtimes --
809 0 if successful, nonzero otherwise.
810 ---------------------------------------------------------------*/
idaLsATimes(void * ida_mem,N_Vector v,N_Vector z)811 int idaLsATimes(void *ida_mem, N_Vector v, N_Vector z)
812 {
813 IDAMem IDA_mem;
814 IDALsMem idals_mem;
815 int retval;
816
817 /* access IDALsMem structure */
818 retval = idaLs_AccessLMem(ida_mem, "idaLsATimes",
819 &IDA_mem, &idals_mem);
820 if (retval != IDALS_SUCCESS) return(retval);
821
822 /* call Jacobian-times-vector product routine
823 (either user-supplied or internal DQ) */
824 retval = idals_mem->jtimes(IDA_mem->ida_tn, idals_mem->ycur,
825 idals_mem->ypcur, idals_mem->rcur,
826 v, z, IDA_mem->ida_cj,
827 idals_mem->jt_data, idals_mem->ytemp,
828 idals_mem->yptemp);
829 idals_mem->njtimes++;
830 return(retval);
831 }
832
833
834 /*---------------------------------------------------------------
835 idaLsPSetup:
836
837 This routine interfaces between the generic iterative linear
838 solvers and the user's psetup routine. It passes to psetup all
839 required state information from ida_mem. Its return value
840 is the same as that returned by psetup. Note that the generic
841 iterative linear solvers guarantee that idaLsPSetup will only
842 be called in the case that the user's psetup routine is non-NULL.
843 ---------------------------------------------------------------*/
idaLsPSetup(void * ida_mem)844 int idaLsPSetup(void *ida_mem)
845 {
846 IDAMem IDA_mem;
847 IDALsMem idals_mem;
848 int retval;
849
850 /* access IDALsMem structure */
851 retval = idaLs_AccessLMem(ida_mem, "idaLsPSetup",
852 &IDA_mem, &idals_mem);
853 if (retval != IDALS_SUCCESS) return(retval);
854
855 /* Call user pset routine to update preconditioner and possibly
856 reset jcur (pass !jbad as update suggestion) */
857 retval = idals_mem->pset(IDA_mem->ida_tn, idals_mem->ycur,
858 idals_mem->ypcur, idals_mem->rcur,
859 IDA_mem->ida_cj, idals_mem->pdata);
860 idals_mem->npe++;
861 return(retval);
862 }
863
864
865 /*---------------------------------------------------------------
866 idaLsPSolve:
867
868 This routine interfaces between the generic SUNLinSolSolve
869 routine and the user's psolve routine. It passes to psolve all
870 required state information from ida_mem. Its return value is
871 the same as that returned by psolve. Note that the generic
872 SUNLinSol solver guarantees that IDASilsPSolve will not be
873 called in the case in which preconditioning is not done. This
874 is the only case in which the user's psolve routine is allowed
875 to be NULL.
876 ---------------------------------------------------------------*/
idaLsPSolve(void * ida_mem,N_Vector r,N_Vector z,realtype tol,int lr)877 int idaLsPSolve(void *ida_mem, N_Vector r, N_Vector z, realtype tol, int lr)
878 {
879 IDAMem IDA_mem;
880 IDALsMem idals_mem;
881 int retval;
882
883 /* access IDALsMem structure */
884 retval = idaLs_AccessLMem(ida_mem, "idaLsPSolve",
885 &IDA_mem, &idals_mem);
886 if (retval != IDALS_SUCCESS) return(retval);
887
888 /* call the user-supplied psolve routine, and accumulate count */
889 retval = idals_mem->psolve(IDA_mem->ida_tn, idals_mem->ycur,
890 idals_mem->ypcur, idals_mem->rcur,
891 r, z, IDA_mem->ida_cj, tol,
892 idals_mem->pdata);
893 idals_mem->nps++;
894 return(retval);
895 }
896
897
898 /*---------------------------------------------------------------
899 idaLsDQJac:
900
901 This routine is a wrapper for the Dense and Band
902 implementations of the difference quotient Jacobian
903 approximation routines.
904 ---------------------------------------------------------------*/
idaLsDQJac(realtype t,realtype c_j,N_Vector y,N_Vector yp,N_Vector r,SUNMatrix Jac,void * ida_mem,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)905 int idaLsDQJac(realtype t, realtype c_j, N_Vector y, N_Vector yp,
906 N_Vector r, SUNMatrix Jac, void *ida_mem,
907 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
908 {
909 int retval;
910 IDAMem IDA_mem;
911 IDA_mem = (IDAMem) ida_mem;
912
913 /* access IDAMem structure */
914 if (ida_mem == NULL) {
915 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASLS",
916 "idaLsDQJac", MSG_LS_IDAMEM_NULL);
917 return(IDALS_MEM_NULL);
918 }
919
920 /* verify that Jac is non-NULL */
921 if (Jac == NULL) {
922 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDASLS",
923 "idaLsDQJac", MSG_LS_LMEM_NULL);
924 return(IDALS_LMEM_NULL);
925 }
926
927 /* Verify that N_Vector supports required operations */
928 if (IDA_mem->ida_tempv1->ops->nvcloneempty == NULL ||
929 IDA_mem->ida_tempv1->ops->nvlinearsum == NULL ||
930 IDA_mem->ida_tempv1->ops->nvdestroy == NULL ||
931 IDA_mem->ida_tempv1->ops->nvscale == NULL ||
932 IDA_mem->ida_tempv1->ops->nvgetarraypointer == NULL ||
933 IDA_mem->ida_tempv1->ops->nvsetarraypointer == NULL) {
934 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS",
935 "idaLsDQJac", MSG_LS_BAD_NVECTOR);
936 return(IDALS_ILL_INPUT);
937 }
938
939 /* Call the matrix-structure-specific DQ approximation routine */
940 if (SUNMatGetID(Jac) == SUNMATRIX_DENSE) {
941 retval = idaLsDenseDQJac(t, c_j, y, yp, r, Jac, IDA_mem, tmp1);
942 } else if (SUNMatGetID(Jac) == SUNMATRIX_BAND) {
943 retval = idaLsBandDQJac(t, c_j, y, yp, r, Jac, IDA_mem, tmp1, tmp2, tmp3);
944 } else {
945 IDAProcessError(IDA_mem, IDA_ILL_INPUT, "IDASLS",
946 "idaLsDQJac",
947 "unrecognized matrix type for idaLsDQJac");
948 retval = IDA_ILL_INPUT;
949 }
950 return(retval);
951 }
952
953
954 /*---------------------------------------------------------------
955 idaLsDenseDQJac
956
957 This routine generates a dense difference quotient approximation
958 to the Jacobian F_y + c_j*F_y'. It assumes a dense SUNmatrix
959 input (stored column-wise, and that elements within each column
960 are contiguous). The address of the jth column of J is obtained
961 via the function SUNDenseMatrix_Column() and this pointer is
962 associated with an N_Vector using the
963 N_VGetArrayPointer/N_VSetArrayPointer functions. Finally, the
964 actual computation of the jth column of the Jacobian is
965 done with a call to N_VLinearSum.
966 ---------------------------------------------------------------*/
idaLsDenseDQJac(realtype tt,realtype c_j,N_Vector yy,N_Vector yp,N_Vector rr,SUNMatrix Jac,IDAMem IDA_mem,N_Vector tmp1)967 int idaLsDenseDQJac(realtype tt, realtype c_j, N_Vector yy,
968 N_Vector yp, N_Vector rr, SUNMatrix Jac,
969 IDAMem IDA_mem, N_Vector tmp1)
970 {
971 realtype inc, inc_inv, yj, ypj, srur, conj;
972 realtype *y_data, *yp_data, *ewt_data, *cns_data = NULL;
973 N_Vector rtemp, jthCol;
974 sunindextype j, N;
975 IDALsMem idals_mem;
976 int retval = 0;
977
978 /* access LsMem interface structure */
979 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
980
981 /* access matrix dimension */
982 N = SUNDenseMatrix_Columns(Jac);
983
984 /* Rename work vectors for readibility */
985 rtemp = tmp1;
986
987 /* Create an empty vector for matrix column calculations */
988 jthCol = N_VCloneEmpty(tmp1);
989
990 /* Obtain pointers to the data for ewt, yy, yp. */
991 ewt_data = N_VGetArrayPointer(IDA_mem->ida_ewt);
992 y_data = N_VGetArrayPointer(yy);
993 yp_data = N_VGetArrayPointer(yp);
994 if(IDA_mem->ida_constraintsSet)
995 cns_data = N_VGetArrayPointer(IDA_mem->ida_constraints);
996
997 srur = SUNRsqrt(IDA_mem->ida_uround);
998
999 for (j=0; j < N; j++) {
1000
1001 /* Generate the jth col of J(tt,yy,yp) as delta(F)/delta(y_j). */
1002
1003 /* Set data address of jthCol, and save y_j and yp_j values. */
1004 N_VSetArrayPointer(SUNDenseMatrix_Column(Jac,j), jthCol);
1005 yj = y_data[j];
1006 ypj = yp_data[j];
1007
1008 /* Set increment inc to y_j based on sqrt(uround)*abs(y_j), with
1009 adjustments using yp_j and ewt_j if this is small, and a further
1010 adjustment to give it the same sign as hh*yp_j. */
1011
1012 inc = SUNMAX( srur * SUNMAX( SUNRabs(yj), SUNRabs(IDA_mem->ida_hh*ypj) ),
1013 ONE/ewt_data[j] );
1014
1015 if (IDA_mem->ida_hh*ypj < ZERO) inc = -inc;
1016 inc = (yj + inc) - yj;
1017
1018 /* Adjust sign(inc) again if y_j has an inequality constraint. */
1019 if (IDA_mem->ida_constraintsSet) {
1020 conj = cns_data[j];
1021 if (SUNRabs(conj) == ONE) {if((yj+inc)*conj < ZERO) inc = -inc;}
1022 else if (SUNRabs(conj) == TWO) {if((yj+inc)*conj <= ZERO) inc = -inc;}
1023 }
1024
1025 /* Increment y_j and yp_j, call res, and break on error return. */
1026 y_data[j] += inc;
1027 yp_data[j] += c_j*inc;
1028
1029 retval = idals_mem->jt_res(tt, yy, yp, rtemp, IDA_mem->ida_user_data);
1030 idals_mem->nreDQ++;
1031 if (retval != 0) break;
1032
1033 /* Construct difference quotient in jthCol */
1034 inc_inv = ONE/inc;
1035 N_VLinearSum(inc_inv, rtemp, -inc_inv, rr, jthCol);
1036
1037 /* reset y_j, yp_j */
1038 y_data[j] = yj;
1039 yp_data[j] = ypj;
1040 }
1041
1042 /* Destroy jthCol vector */
1043 N_VSetArrayPointer(NULL, jthCol); /* SHOULDN'T BE NEEDED */
1044 N_VDestroy(jthCol);
1045
1046 return(retval);
1047 }
1048
1049
1050 /*---------------------------------------------------------------
1051 idaLsBandDQJac
1052
1053 This routine generates a banded difference quotient approximation
1054 JJ to the DAE system Jacobian J. It assumes a band SUNMatrix
1055 input (stored column-wise, and that elements within each column
1056 are contiguous). This makes it possible to get the address
1057 of a column of JJ via the function SUNBandMatrix_Column(). The
1058 columns of the Jacobian are constructed using mupper + mlower + 1
1059 calls to the res routine, and appropriate differencing.
1060 The return value is either IDABAND_SUCCESS = 0, or the nonzero
1061 value returned by the res routine, if any.
1062 ---------------------------------------------------------------*/
idaLsBandDQJac(realtype tt,realtype c_j,N_Vector yy,N_Vector yp,N_Vector rr,SUNMatrix Jac,IDAMem IDA_mem,N_Vector tmp1,N_Vector tmp2,N_Vector tmp3)1063 int idaLsBandDQJac(realtype tt, realtype c_j, N_Vector yy,
1064 N_Vector yp, N_Vector rr, SUNMatrix Jac,
1065 IDAMem IDA_mem, N_Vector tmp1, N_Vector tmp2,
1066 N_Vector tmp3)
1067 {
1068 realtype inc, inc_inv, yj, ypj, srur, conj, ewtj;
1069 realtype *y_data, *yp_data, *ewt_data, *cns_data = NULL;
1070 realtype *ytemp_data, *yptemp_data, *rtemp_data, *r_data, *col_j;
1071 N_Vector rtemp, ytemp, yptemp;
1072 sunindextype i, j, i1, i2, width, ngroups, group;
1073 sunindextype N, mupper, mlower;
1074 IDALsMem idals_mem;
1075 int retval = 0;
1076
1077 /* access LsMem interface structure */
1078 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
1079
1080 /* access matrix dimensions */
1081 N = SUNBandMatrix_Columns(Jac);
1082 mupper = SUNBandMatrix_UpperBandwidth(Jac);
1083 mlower = SUNBandMatrix_LowerBandwidth(Jac);
1084
1085 /* Rename work vectors for use as temporary values of r, y and yp */
1086 rtemp = tmp1;
1087 ytemp = tmp2;
1088 yptemp= tmp3;
1089
1090 /* Obtain pointers to the data for all eight vectors used. */
1091 ewt_data = N_VGetArrayPointer(IDA_mem->ida_ewt);
1092 r_data = N_VGetArrayPointer(rr);
1093 y_data = N_VGetArrayPointer(yy);
1094 yp_data = N_VGetArrayPointer(yp);
1095 rtemp_data = N_VGetArrayPointer(rtemp);
1096 ytemp_data = N_VGetArrayPointer(ytemp);
1097 yptemp_data = N_VGetArrayPointer(yptemp);
1098 if (IDA_mem->ida_constraintsSet)
1099 cns_data = N_VGetArrayPointer(IDA_mem->ida_constraints);
1100
1101 /* Initialize ytemp and yptemp. */
1102 N_VScale(ONE, yy, ytemp);
1103 N_VScale(ONE, yp, yptemp);
1104
1105 /* Compute miscellaneous values for the Jacobian computation. */
1106 srur = SUNRsqrt(IDA_mem->ida_uround);
1107 width = mlower + mupper + 1;
1108 ngroups = SUNMIN(width, N);
1109
1110 /* Loop over column groups. */
1111 for (group=1; group <= ngroups; group++) {
1112
1113 /* Increment all yy[j] and yp[j] for j in this group. */
1114 for (j=group-1; j<N; j+=width) {
1115 yj = y_data[j];
1116 ypj = yp_data[j];
1117 ewtj = ewt_data[j];
1118
1119 /* Set increment inc to yj based on sqrt(uround)*abs(yj), with
1120 adjustments using ypj and ewtj if this is small, and a further
1121 adjustment to give it the same sign as hh*ypj. */
1122 inc = SUNMAX( srur * SUNMAX( SUNRabs(yj), SUNRabs(IDA_mem->ida_hh*ypj) ),
1123 ONE/ewtj );
1124 if (IDA_mem->ida_hh*ypj < ZERO) inc = -inc;
1125 inc = (yj + inc) - yj;
1126
1127 /* Adjust sign(inc) again if yj has an inequality constraint. */
1128 if (IDA_mem->ida_constraintsSet) {
1129 conj = cns_data[j];
1130 if (SUNRabs(conj) == ONE) {if((yj+inc)*conj < ZERO) inc = -inc;}
1131 else if (SUNRabs(conj) == TWO) {if((yj+inc)*conj <= ZERO) inc = -inc;}
1132 }
1133
1134 /* Increment yj and ypj. */
1135 ytemp_data[j] += inc;
1136 yptemp_data[j] += IDA_mem->ida_cj*inc;
1137 }
1138
1139 /* Call res routine with incremented arguments. */
1140 retval = IDA_mem->ida_res(tt, ytemp, yptemp, rtemp, IDA_mem->ida_user_data);
1141 idals_mem->nreDQ++;
1142 if (retval != 0) break;
1143
1144 /* Loop over the indices j in this group again. */
1145 for (j=group-1; j<N; j+=width) {
1146
1147 /* Reset ytemp and yptemp components that were perturbed. */
1148 yj = ytemp_data[j] = y_data[j];
1149 ypj = yptemp_data[j] = yp_data[j];
1150 col_j = SUNBandMatrix_Column(Jac, j);
1151 ewtj = ewt_data[j];
1152
1153 /* Set increment inc exactly as above. */
1154 inc = SUNMAX( srur * SUNMAX( SUNRabs(yj), SUNRabs(IDA_mem->ida_hh*ypj) ),
1155 ONE/ewtj );
1156 if (IDA_mem->ida_hh*ypj < ZERO) inc = -inc;
1157 inc = (yj + inc) - yj;
1158 if (IDA_mem->ida_constraintsSet) {
1159 conj = cns_data[j];
1160 if (SUNRabs(conj) == ONE) {if((yj+inc)*conj < ZERO) inc = -inc;}
1161 else if (SUNRabs(conj) == TWO) {if((yj+inc)*conj <= ZERO) inc = -inc;}
1162 }
1163
1164 /* Load the difference quotient Jacobian elements for column j */
1165 inc_inv = ONE/inc;
1166 i1 = SUNMAX(0, j-mupper);
1167 i2 = SUNMIN(j+mlower,N-1);
1168 for (i=i1; i<=i2; i++)
1169 SM_COLUMN_ELEMENT_B(col_j,i,j) = inc_inv * (rtemp_data[i]-r_data[i]);
1170 }
1171 }
1172
1173 return(retval);
1174 }
1175
1176
1177 /*---------------------------------------------------------------
1178 idaLsDQJtimes
1179
1180 This routine generates a difference quotient approximation to
1181 the matrix-vector product z = Jv, where J is the system
1182 Jacobian. The approximation is
1183 Jv = [F(t,y1,yp1) - F(t,y,yp)]/sigma,
1184 where
1185 y1 = y + sigma*v, yp1 = yp + cj*sigma*v,
1186 sigma = sqrt(Neq)*dqincfac.
1187 The return value from the call to res is saved in order to set
1188 the return flag from idaLsSolve.
1189 ---------------------------------------------------------------*/
idaLsDQJtimes(realtype tt,N_Vector yy,N_Vector yp,N_Vector rr,N_Vector v,N_Vector Jv,realtype c_j,void * ida_mem,N_Vector work1,N_Vector work2)1190 int idaLsDQJtimes(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr,
1191 N_Vector v, N_Vector Jv, realtype c_j,
1192 void *ida_mem, N_Vector work1, N_Vector work2)
1193 {
1194 IDAMem IDA_mem;
1195 IDALsMem idals_mem;
1196 N_Vector y_tmp, yp_tmp;
1197 realtype sig, siginv;
1198 int iter, retval;
1199 SUNLinearSolver_ID LSID;
1200
1201 /* access IDALsMem structure */
1202 retval = idaLs_AccessLMem(ida_mem, "idaLsDQJtimes",
1203 &IDA_mem, &idals_mem);
1204 if (retval != IDALS_SUCCESS) return(retval);
1205
1206 LSID = SUNLinSolGetID(idals_mem->LS);
1207 if (LSID == SUNLINEARSOLVER_SPGMR || LSID == SUNLINEARSOLVER_SPFGMR)
1208 sig = idals_mem->nrmfac * idals_mem->dqincfac;
1209 else
1210 sig = idals_mem->dqincfac / N_VWrmsNorm(v, IDA_mem->ida_ewt);
1211
1212 /* Rename work1 and work2 for readibility */
1213 y_tmp = work1;
1214 yp_tmp = work2;
1215
1216 for (iter=0; iter<MAX_ITERS; iter++) {
1217
1218 /* Set y_tmp = yy + sig*v, yp_tmp = yp + cj*sig*v. */
1219 N_VLinearSum(sig, v, ONE, yy, y_tmp);
1220 N_VLinearSum(c_j*sig, v, ONE, yp, yp_tmp);
1221
1222 /* Call res for Jv = F(t, y_tmp, yp_tmp), and return if it failed. */
1223 retval = IDA_mem->ida_res(tt, y_tmp, yp_tmp, Jv, IDA_mem->ida_user_data);
1224 idals_mem->nreDQ++;
1225 if (retval == 0) break;
1226 if (retval < 0) return(-1);
1227
1228 sig *= PT25;
1229 }
1230
1231 if (retval > 0) return(+1);
1232
1233 /* Set Jv to [Jv - rr]/sig and return. */
1234 siginv = ONE/sig;
1235 N_VLinearSum(siginv, Jv, -siginv, rr, Jv);
1236
1237 return(0);
1238 }
1239
1240
1241 /*---------------------------------------------------------------
1242 idaLsInitialize
1243
1244 This routine performs remaining initializations specific
1245 to the iterative linear solver interface (and solver itself)
1246 ---------------------------------------------------------------*/
idaLsInitialize(IDAMem IDA_mem)1247 int idaLsInitialize(IDAMem IDA_mem)
1248 {
1249 IDALsMem idals_mem;
1250 int retval;
1251
1252 /* access IDALsMem structure */
1253 if (IDA_mem->ida_lmem == NULL) {
1254 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDASLS",
1255 "idaLsInitialize", MSG_LS_LMEM_NULL);
1256 return(IDALS_LMEM_NULL);
1257 }
1258 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
1259
1260
1261 /* Test for valid combinations of matrix & Jacobian routines: */
1262 if (idals_mem->J == NULL) {
1263
1264 /* If SUNMatrix A is NULL: ensure 'jac' function pointer is NULL */
1265 idals_mem->jacDQ = SUNFALSE;
1266 idals_mem->jac = NULL;
1267 idals_mem->J_data = NULL;
1268
1269 } else if (idals_mem->jacDQ) {
1270
1271 /* If J is non-NULL, and 'jac' is not user-supplied:
1272 - if J is dense or band, ensure that our DQ approx. is used
1273 - otherwise => error */
1274 retval = 0;
1275 if (idals_mem->J->ops->getid) {
1276
1277 if ( (SUNMatGetID(idals_mem->J) == SUNMATRIX_DENSE) ||
1278 (SUNMatGetID(idals_mem->J) == SUNMATRIX_BAND) ) {
1279 idals_mem->jac = idaLsDQJac;
1280 idals_mem->J_data = IDA_mem;
1281 } else {
1282 retval++;
1283 }
1284
1285 } else {
1286 retval++;
1287 }
1288 if (retval) {
1289 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS", "idaLsInitialize",
1290 "No Jacobian constructor available for SUNMatrix type");
1291 idals_mem->last_flag = IDALS_ILL_INPUT;
1292 return(IDALS_ILL_INPUT);
1293 }
1294
1295 } else {
1296
1297 /* If J is non-NULL, and 'jac' is user-supplied,
1298 reset J_data pointer (just in case) */
1299 idals_mem->J_data = IDA_mem->ida_user_data;
1300 }
1301
1302 /* reset counters */
1303 idaLsInitializeCounters(idals_mem);
1304
1305 /* Set Jacobian-related fields, based on jtimesDQ */
1306 if (idals_mem->jtimesDQ) {
1307 idals_mem->jtsetup = NULL;
1308 idals_mem->jtimes = idaLsDQJtimes;
1309 idals_mem->jt_data = IDA_mem;
1310 } else {
1311 idals_mem->jt_data = IDA_mem->ida_user_data;
1312 }
1313
1314 /* if J is NULL and psetup is not present, then idaLsSetup does
1315 not need to be called, so set the lsetup function to NULL */
1316 if ( (idals_mem->J == NULL) && (idals_mem->pset == NULL) )
1317 IDA_mem->ida_lsetup = NULL;
1318
1319 /* Call LS initialize routine */
1320 idals_mem->last_flag = SUNLinSolInitialize(idals_mem->LS);
1321 return(idals_mem->last_flag);
1322 }
1323
1324
1325 /*---------------------------------------------------------------
1326 idaLsSetup
1327
1328 This calls the Jacobian evaluation routine (if using a SUNMatrix
1329 object), updates counters, and calls the LS 'setup' routine to
1330 prepare for subsequent calls to the LS 'solve' routine.
1331 ---------------------------------------------------------------*/
idaLsSetup(IDAMem IDA_mem,N_Vector y,N_Vector yp,N_Vector r,N_Vector vt1,N_Vector vt2,N_Vector vt3)1332 int idaLsSetup(IDAMem IDA_mem, N_Vector y, N_Vector yp, N_Vector r,
1333 N_Vector vt1, N_Vector vt2, N_Vector vt3)
1334 {
1335 IDALsMem idals_mem;
1336 int retval;
1337
1338 /* access IDALsMem structure */
1339 if (IDA_mem->ida_lmem == NULL) {
1340 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDASLS",
1341 "idaLsSetup", MSG_LS_LMEM_NULL);
1342 return(IDALS_LMEM_NULL);
1343 }
1344 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
1345
1346 /* Set IDALs N_Vector pointers to inputs */
1347 idals_mem->ycur = y;
1348 idals_mem->ypcur = yp;
1349 idals_mem->rcur = r;
1350
1351 /* recompute if J if it is non-NULL */
1352 if (idals_mem->J) {
1353
1354 /* Increment nje counter. */
1355 idals_mem->nje++;
1356
1357 /* Clear the linear system matrix if necessary */
1358 if (SUNLinSolGetType(idals_mem->LS) == SUNLINEARSOLVER_DIRECT) {
1359 retval = SUNMatZero(idals_mem->J);
1360 if (retval != 0) {
1361 IDAProcessError(IDA_mem, IDALS_SUNMAT_FAIL, "IDASLS",
1362 "idaLsSetup", MSG_LS_MATZERO_FAILED);
1363 idals_mem->last_flag = IDALS_SUNMAT_FAIL;
1364 return(idals_mem->last_flag);
1365 }
1366 }
1367
1368 /* Call Jacobian routine */
1369 retval = idals_mem->jac(IDA_mem->ida_tn, IDA_mem->ida_cj, y,
1370 yp, r, idals_mem->J,
1371 idals_mem->J_data, vt1, vt2, vt3);
1372 if (retval < 0) {
1373 IDAProcessError(IDA_mem, IDALS_JACFUNC_UNRECVR, "IDASLS",
1374 "idaLsSetup", MSG_LS_JACFUNC_FAILED);
1375 idals_mem->last_flag = IDALS_JACFUNC_UNRECVR;
1376 return(-1);
1377 }
1378 if (retval > 0) {
1379 idals_mem->last_flag = IDALS_JACFUNC_RECVR;
1380 return(1);
1381 }
1382
1383 }
1384
1385 /* Call LS setup routine -- the LS will call idaLsPSetup if applicable */
1386 idals_mem->last_flag = SUNLinSolSetup(idals_mem->LS, idals_mem->J);
1387 return(idals_mem->last_flag);
1388 }
1389
1390
1391 /*---------------------------------------------------------------
1392 idaLsSolve
1393
1394 This routine interfaces between IDA and the generic
1395 SUNLinearSolver object LS, by setting the appropriate tolerance
1396 and scaling vectors, calling the solver, accumulating
1397 statistics from the solve for use/reporting by IDA, and scaling
1398 the result if using a non-NULL SUNMatrix and cjratio does not
1399 equal one.
1400 ---------------------------------------------------------------*/
idaLsSolve(IDAMem IDA_mem,N_Vector b,N_Vector weight,N_Vector ycur,N_Vector ypcur,N_Vector rescur)1401 int idaLsSolve(IDAMem IDA_mem, N_Vector b, N_Vector weight,
1402 N_Vector ycur, N_Vector ypcur, N_Vector rescur)
1403 {
1404 IDALsMem idals_mem;
1405 int nli_inc, retval;
1406 realtype tol, w_mean;
1407
1408 /* access IDALsMem structure */
1409 if (IDA_mem->ida_lmem == NULL) {
1410 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDASLS",
1411 "idaLsSolve", MSG_LS_LMEM_NULL);
1412 return(IDALS_LMEM_NULL);
1413 }
1414 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
1415
1416 /* If the linear solver is iterative: set convergence test constant tol,
1417 in terms of the Newton convergence test constant epsNewt and safety
1418 factors. The factor nrmfac assures that the convergence test is
1419 applied to the WRMS norm of the residual vector, rather than the
1420 weighted L2 norm. */
1421 if (idals_mem->iterative) {
1422 tol = idals_mem->nrmfac * idals_mem->eplifac * IDA_mem->ida_epsNewt;
1423 } else {
1424 tol = ZERO;
1425 }
1426
1427 /* Set vectors ycur, ypcur and rcur for use by the Atimes and
1428 Psolve interface routines */
1429 idals_mem->ycur = ycur;
1430 idals_mem->ypcur = ypcur;
1431 idals_mem->rcur = rescur;
1432
1433 /* Set scaling vectors for LS to use (if applicable) */
1434 if (idals_mem->LS->ops->setscalingvectors) {
1435 retval = SUNLinSolSetScalingVectors(idals_mem->LS, weight, weight);
1436 if (retval != SUNLS_SUCCESS) {
1437 IDAProcessError(IDA_mem, IDALS_SUNLS_FAIL, "IDASLS", "idaLsSolve",
1438 "Error in calling SUNLinSolSetScalingVectors");
1439 idals_mem->last_flag = IDALS_SUNLS_FAIL;
1440 return(idals_mem->last_flag);
1441 }
1442
1443 /* If solver is iterative and does not support scaling vectors, update the
1444 tolerance in an attempt to account for weight vector. We make the
1445 following assumptions:
1446 1. w_i = w_mean, for i=0,...,n-1 (i.e. the weights are homogeneous)
1447 2. the linear solver uses a basic 2-norm to measure convergence
1448 Hence (using the notation from sunlinsol_spgmr.h, with S = diag(w)),
1449 || bbar - Abar xbar ||_2 < tol
1450 <=> || S b - S A x ||_2 < tol
1451 <=> || S (b - A x) ||_2 < tol
1452 <=> \sum_{i=0}^{n-1} (w_i (b - A x)_i)^2 < tol^2
1453 <=> w_mean^2 \sum_{i=0}^{n-1} (b - A x_i)^2 < tol^2
1454 <=> \sum_{i=0}^{n-1} (b - A x_i)^2 < tol^2 / w_mean^2
1455 <=> || b - A x ||_2 < tol / w_mean
1456 So we compute w_mean = ||w||_RMS and scale the desired tolerance accordingly. */
1457 } else if (idals_mem->iterative) {
1458
1459 N_VConst(ONE, idals_mem->x);
1460 w_mean = N_VWrmsNorm(weight, idals_mem->x);
1461 tol /= w_mean;
1462
1463 }
1464
1465 /* Set initial guess x = 0 to LS */
1466 N_VConst(ZERO, idals_mem->x);
1467
1468 /* If a user-provided jtsetup routine is supplied, call that here */
1469 if (idals_mem->jtsetup) {
1470 idals_mem->last_flag = idals_mem->jtsetup(IDA_mem->ida_tn, ycur, ypcur, rescur,
1471 IDA_mem->ida_cj, idals_mem->jt_data);
1472 idals_mem->njtsetup++;
1473 if (idals_mem->last_flag != 0) {
1474 IDAProcessError(IDA_mem, retval, "IDASLS",
1475 "idaLsSolve", MSG_LS_JTSETUP_FAILED);
1476 return(idals_mem->last_flag);
1477 }
1478 }
1479
1480 /* Call solver */
1481 retval = SUNLinSolSolve(idals_mem->LS, idals_mem->J,
1482 idals_mem->x, b, tol);
1483
1484 /* Copy appropriate result to b (depending on solver type) */
1485 if (idals_mem->iterative) {
1486
1487 /* Retrieve solver statistics */
1488 nli_inc = SUNLinSolNumIters(idals_mem->LS);
1489
1490 /* Copy x (or preconditioned residual vector if no iterations required) to b */
1491 if (nli_inc == 0) N_VScale(ONE, SUNLinSolResid(idals_mem->LS), b);
1492 else N_VScale(ONE, idals_mem->x, b);
1493
1494 /* Increment nli counter */
1495 idals_mem->nli += nli_inc;
1496
1497 } else {
1498
1499 /* Copy x to b */
1500 N_VScale(ONE, idals_mem->x, b);
1501
1502 }
1503
1504 /* If using a direct or matrix-iterative solver, scale the correction to
1505 account for change in cj */
1506 if (idals_mem->scalesol && (IDA_mem->ida_cjratio != ONE))
1507 N_VScale(TWO/(ONE + IDA_mem->ida_cjratio), b, b);
1508
1509 /* Increment ncfl counter */
1510 if (retval != SUNLS_SUCCESS) idals_mem->ncfl++;
1511
1512 /* Interpret solver return value */
1513 idals_mem->last_flag = retval;
1514
1515 switch(retval) {
1516
1517 case SUNLS_SUCCESS:
1518 return(0);
1519 break;
1520 case SUNLS_RES_REDUCED:
1521 case SUNLS_CONV_FAIL:
1522 case SUNLS_PSOLVE_FAIL_REC:
1523 case SUNLS_PACKAGE_FAIL_REC:
1524 case SUNLS_QRFACT_FAIL:
1525 case SUNLS_LUFACT_FAIL:
1526 return(1);
1527 break;
1528 case SUNLS_MEM_NULL:
1529 case SUNLS_ILL_INPUT:
1530 case SUNLS_MEM_FAIL:
1531 case SUNLS_GS_FAIL:
1532 case SUNLS_QRSOL_FAIL:
1533 return(-1);
1534 break;
1535 case SUNLS_PACKAGE_FAIL_UNREC:
1536 IDAProcessError(IDA_mem, SUNLS_PACKAGE_FAIL_UNREC, "IDASLS",
1537 "idaLsSolve",
1538 "Failure in SUNLinSol external package");
1539 return(-1);
1540 break;
1541 case SUNLS_PSOLVE_FAIL_UNREC:
1542 IDAProcessError(IDA_mem, SUNLS_PSOLVE_FAIL_UNREC, "IDASLS",
1543 "idaLsSolve", MSG_LS_PSOLVE_FAILED);
1544 return(-1);
1545 break;
1546 }
1547
1548 return(0);
1549 }
1550
1551
1552 /*---------------------------------------------------------------
1553 idaLsPerf: accumulates performance statistics information
1554 for IDA
1555 ---------------------------------------------------------------*/
idaLsPerf(IDAMem IDA_mem,int perftask)1556 int idaLsPerf(IDAMem IDA_mem, int perftask)
1557 {
1558 IDALsMem idals_mem;
1559 realtype rcfn, rcfl;
1560 long int nstd, nnid;
1561 booleantype lcfn, lcfl;
1562
1563 /* access IDALsMem structure */
1564 if (IDA_mem->ida_lmem == NULL) {
1565 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDASLS",
1566 "idaLsPerf", MSG_LS_LMEM_NULL);
1567 return(IDALS_LMEM_NULL);
1568 }
1569 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
1570
1571 /* when perftask == 0, store current performance statistics */
1572 if (perftask == 0) {
1573 idals_mem->nst0 = IDA_mem->ida_nst;
1574 idals_mem->nni0 = IDA_mem->ida_nni;
1575 idals_mem->ncfn0 = IDA_mem->ida_ncfn;
1576 idals_mem->ncfl0 = idals_mem->ncfl;
1577 idals_mem->nwarn = 0;
1578 return(0);
1579 }
1580
1581 /* Compute statistics since last call
1582
1583 Note: the performance monitor that checked whether the average
1584 number of linear iterations was too close to maxl has been
1585 removed, since the 'maxl' value is no longer owned by the
1586 IDALs interface.
1587 */
1588 nstd = IDA_mem->ida_nst - idals_mem->nst0;
1589 nnid = IDA_mem->ida_nni - idals_mem->nni0;
1590 if (nstd == 0 || nnid == 0) return(0);
1591
1592 rcfn = ((realtype) (IDA_mem->ida_ncfn - idals_mem->ncfn0)) / ((realtype) nstd);
1593 rcfl = ((realtype) (idals_mem->ncfl - idals_mem->ncfl0)) / ((realtype) nnid);
1594 lcfn = (rcfn > PT9);
1595 lcfl = (rcfl > PT9);
1596 if (!(lcfn || lcfl)) return(0);
1597 idals_mem->nwarn++;
1598 if (idals_mem->nwarn > 10) return(1);
1599 if (lcfn)
1600 IDAProcessError(IDA_mem, IDA_WARNING, "IDASLS", "idaLsPerf",
1601 MSG_LS_CFN_WARN, IDA_mem->ida_tn, rcfn);
1602 if (lcfl)
1603 IDAProcessError(IDA_mem, IDA_WARNING, "IDASLS", "idaLsPerf",
1604 MSG_LS_CFL_WARN, IDA_mem->ida_tn, rcfl);
1605 return(0);
1606 }
1607
1608
1609 /*---------------------------------------------------------------
1610 idaLsFree frees memory associates with the IDALs system
1611 solver interface.
1612 ---------------------------------------------------------------*/
idaLsFree(IDAMem IDA_mem)1613 int idaLsFree(IDAMem IDA_mem)
1614 {
1615 IDALsMem idals_mem;
1616
1617 /* Return immediately if IDA_mem or IDA_mem->ida_lmem are NULL */
1618 if (IDA_mem == NULL) return (IDALS_SUCCESS);
1619 if (IDA_mem->ida_lmem == NULL) return(IDALS_SUCCESS);
1620 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
1621
1622 /* Free N_Vector memory */
1623 if (idals_mem->ytemp) {
1624 N_VDestroy(idals_mem->ytemp);
1625 idals_mem->ytemp = NULL;
1626 }
1627 if (idals_mem->yptemp) {
1628 N_VDestroy(idals_mem->yptemp);
1629 idals_mem->yptemp = NULL;
1630 }
1631 if (idals_mem->x) {
1632 N_VDestroy(idals_mem->x);
1633 idals_mem->x = NULL;
1634 }
1635
1636 /* Nullify other N_Vector pointers */
1637 idals_mem->ycur = NULL;
1638 idals_mem->ypcur = NULL;
1639 idals_mem->rcur = NULL;
1640
1641 /* Nullify SUNMatrix pointer */
1642 idals_mem->J = NULL;
1643
1644 /* Free preconditioner memory (if applicable) */
1645 if (idals_mem->pfree) idals_mem->pfree(IDA_mem);
1646
1647 /* free IDALs interface structure */
1648 free(IDA_mem->ida_lmem);
1649
1650 return(IDALS_SUCCESS);
1651 }
1652
1653
1654 /*---------------------------------------------------------------
1655 idaLsInitializeCounters resets all counters from an
1656 IDALsMem structure.
1657 ---------------------------------------------------------------*/
idaLsInitializeCounters(IDALsMem idals_mem)1658 int idaLsInitializeCounters(IDALsMem idals_mem)
1659 {
1660 idals_mem->nje = 0;
1661 idals_mem->nreDQ = 0;
1662 idals_mem->npe = 0;
1663 idals_mem->nli = 0;
1664 idals_mem->nps = 0;
1665 idals_mem->ncfl = 0;
1666 idals_mem->njtsetup = 0;
1667 idals_mem->njtimes = 0;
1668 return(0);
1669 }
1670
1671
1672 /*---------------------------------------------------------------
1673 idaLs_AccessLMem
1674
1675 This routine unpacks the IDA_mem and idals_mem structures from
1676 the void* ida_mem pointer. If either is missing it returns
1677 IDALS_MEM_NULL or IDALS_LMEM_NULL.
1678 ---------------------------------------------------------------*/
idaLs_AccessLMem(void * ida_mem,const char * fname,IDAMem * IDA_mem,IDALsMem * idals_mem)1679 int idaLs_AccessLMem(void* ida_mem, const char* fname,
1680 IDAMem* IDA_mem, IDALsMem* idals_mem)
1681 {
1682 if (ida_mem==NULL) {
1683 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASLS",
1684 fname, MSG_LS_IDAMEM_NULL);
1685 return(IDALS_MEM_NULL);
1686 }
1687 *IDA_mem = (IDAMem) ida_mem;
1688 if ((*IDA_mem)->ida_lmem==NULL) {
1689 IDAProcessError(*IDA_mem, IDALS_LMEM_NULL, "IDASLS",
1690 fname, MSG_LS_LMEM_NULL);
1691 return(IDALS_LMEM_NULL);
1692 }
1693 *idals_mem = (IDALsMem) (*IDA_mem)->ida_lmem;
1694 return(IDALS_SUCCESS);
1695 }
1696
1697
1698
1699 /*================================================================
1700 PART II - backward problems
1701 ================================================================*/
1702
1703
1704 /*---------------------------------------------------------------
1705 IDASLS Exported functions -- Required
1706 ---------------------------------------------------------------*/
1707
1708 /* IDASetLinearSolverB specifies the iterative linear solver
1709 for backward integration */
IDASetLinearSolverB(void * ida_mem,int which,SUNLinearSolver LS,SUNMatrix A)1710 int IDASetLinearSolverB(void *ida_mem, int which,
1711 SUNLinearSolver LS, SUNMatrix A)
1712 {
1713 IDAMem IDA_mem;
1714 IDAadjMem IDAADJ_mem;
1715 IDABMem IDAB_mem;
1716 IDALsMemB idalsB_mem;
1717 void *ida_memB;
1718 int retval;
1719
1720 /* Check if ida_mem exists */
1721 if (ida_mem == NULL) {
1722 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASLS",
1723 "IDASetLinearSolverB", MSG_LS_IDAMEM_NULL);
1724 return(IDALS_MEM_NULL);
1725 }
1726 IDA_mem = (IDAMem) ida_mem;
1727
1728 /* Was ASA initialized? */
1729 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
1730 IDAProcessError(IDA_mem, IDALS_NO_ADJ, "IDASLS",
1731 "IDASetLinearSolverB", MSG_LS_NO_ADJ);
1732 return(IDALS_NO_ADJ);
1733 }
1734 IDAADJ_mem = IDA_mem->ida_adj_mem;
1735
1736 /* Check the value of which */
1737 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
1738 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASLS",
1739 "IDASetLinearSolverB", MSG_LS_BAD_WHICH);
1740 return(IDALS_ILL_INPUT);
1741 }
1742
1743 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
1744 IDAB_mem = IDAADJ_mem->IDAB_mem;
1745 while (IDAB_mem != NULL) {
1746 if( which == IDAB_mem->ida_index ) break;
1747 IDAB_mem = IDAB_mem->ida_next;
1748 }
1749
1750 /* Get memory for IDALsMemRecB */
1751 idalsB_mem = NULL;
1752 idalsB_mem = (IDALsMemB) malloc(sizeof(struct IDALsMemRecB));
1753 if (idalsB_mem == NULL) {
1754 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASLS",
1755 "IDASetLinearSolverB", MSG_LS_MEM_FAIL);
1756 return(IDALS_MEM_FAIL);
1757 }
1758
1759 /* initialize Jacobian and preconditioner functions */
1760 idalsB_mem->jacB = NULL;
1761 idalsB_mem->jacBS = NULL;
1762 idalsB_mem->jtsetupB = NULL;
1763 idalsB_mem->jtsetupBS = NULL;
1764 idalsB_mem->jtimesB = NULL;
1765 idalsB_mem->jtimesBS = NULL;
1766 idalsB_mem->psetB = NULL;
1767 idalsB_mem->psetBS = NULL;
1768 idalsB_mem->psolveB = NULL;
1769 idalsB_mem->psolveBS = NULL;
1770 idalsB_mem->P_dataB = NULL;
1771
1772 /* free any existing system solver attached to IDAB */
1773 if (IDAB_mem->ida_lfree) IDAB_mem->ida_lfree(IDAB_mem);
1774
1775 /* Attach lmemB data and lfreeB function. */
1776 IDAB_mem->ida_lmem = idalsB_mem;
1777 IDAB_mem->ida_lfree = idaLsFreeB;
1778
1779 /* set the linear solver for this backward problem */
1780 ida_memB = (void *)IDAB_mem->IDA_mem;
1781 retval = IDASetLinearSolver(ida_memB, LS, A);
1782 if (retval != IDALS_SUCCESS) {
1783 free(idalsB_mem);
1784 idalsB_mem = NULL;
1785 }
1786
1787 return(retval);
1788 }
1789
1790
1791 /*---------------------------------------------------------------
1792 IDASLS Exported functions -- Optional input/output
1793 ---------------------------------------------------------------*/
1794
IDASetJacFnB(void * ida_mem,int which,IDALsJacFnB jacB)1795 int IDASetJacFnB(void *ida_mem, int which, IDALsJacFnB jacB)
1796 {
1797 IDAMem IDA_mem;
1798 IDAadjMem IDAADJ_mem;
1799 IDABMem IDAB_mem;
1800 IDALsMemB idalsB_mem;
1801 void *ida_memB;
1802 int retval;
1803
1804 /* access relevant memory structures */
1805 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetJacFnB", &IDA_mem,
1806 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
1807 if (retval != IDALS_SUCCESS) return(retval);
1808
1809 /* set jacB function pointer */
1810 idalsB_mem->jacB = jacB;
1811
1812 /* call corresponding routine for IDAB_mem structure */
1813 ida_memB = (void*) IDAB_mem->IDA_mem;
1814 if (jacB != NULL) {
1815 retval = IDASetJacFn(ida_memB, idaLsJacBWrapper);
1816 } else {
1817 retval = IDASetJacFn(ida_memB, NULL);
1818 }
1819
1820 return(retval);
1821 }
1822
1823
IDASetJacFnBS(void * ida_mem,int which,IDALsJacFnBS jacBS)1824 int IDASetJacFnBS(void *ida_mem, int which, IDALsJacFnBS jacBS)
1825 {
1826 IDAMem IDA_mem;
1827 IDAadjMem IDAADJ_mem;
1828 IDABMem IDAB_mem;
1829 IDALsMemB idalsB_mem;
1830 void *ida_memB;
1831 int retval;
1832
1833 /* access relevant memory structures */
1834 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetJacFnBS", &IDA_mem,
1835 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
1836 if (retval != IDALS_SUCCESS) return(retval);
1837
1838 /* set jacBS function pointer */
1839 idalsB_mem->jacBS = jacBS;
1840
1841 /* call corresponding routine for IDAB_mem structure */
1842 ida_memB = (void*) IDAB_mem->IDA_mem;
1843 if (jacBS != NULL) {
1844 retval = IDASetJacFn(ida_memB, idaLsJacBSWrapper);
1845 } else {
1846 retval = IDASetJacFn(ida_memB, NULL);
1847 }
1848
1849 return(retval);
1850 }
1851
1852
IDASetEpsLinB(void * ida_mem,int which,realtype eplifacB)1853 int IDASetEpsLinB(void *ida_mem, int which, realtype eplifacB)
1854 {
1855 IDAadjMem IDAADJ_mem;
1856 IDAMem IDA_mem;
1857 IDABMem IDAB_mem;
1858 IDALsMemB idalsB_mem;
1859 void *ida_memB;
1860 int retval;
1861
1862 /* access relevant memory structures */
1863 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetEpsLinB", &IDA_mem,
1864 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
1865 if (retval != IDALS_SUCCESS) return(retval);
1866
1867 /* call corresponding routine for IDAB_mem structure */
1868 ida_memB = (void *) IDAB_mem->IDA_mem;
1869 return(IDASetEpsLin(ida_memB, eplifacB));
1870 }
1871
1872
IDASetLSNormFactorB(void * ida_mem,int which,realtype nrmfacB)1873 int IDASetLSNormFactorB(void *ida_mem, int which, realtype nrmfacB)
1874 {
1875 IDAadjMem IDAADJ_mem;
1876 IDAMem IDA_mem;
1877 IDABMem IDAB_mem;
1878 IDALsMemB idalsB_mem;
1879 void *ida_memB;
1880 int retval;
1881
1882 /* access relevant memory structures */
1883 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetLSNormFactorB", &IDA_mem,
1884 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
1885 if (retval != IDALS_SUCCESS) return(retval);
1886
1887 /* call corresponding routine for IDAB_mem structure */
1888 ida_memB = (void *) IDAB_mem->IDA_mem;
1889 return(IDASetLSNormFactor(ida_memB, nrmfacB));
1890 }
1891
1892
IDASetLinearSolutionScalingB(void * ida_mem,int which,booleantype onoffB)1893 int IDASetLinearSolutionScalingB(void *ida_mem, int which, booleantype onoffB)
1894 {
1895 IDAadjMem IDAADJ_mem;
1896 IDAMem IDA_mem;
1897 IDABMem IDAB_mem;
1898 IDALsMemB idalsB_mem;
1899 void *ida_memB;
1900 int retval;
1901
1902 /* access relevant memory structures */
1903 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetLinearSolutionScalingB",
1904 &IDA_mem, &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
1905 if (retval != IDALS_SUCCESS) return(retval);
1906
1907 /* call corresponding routine for IDAB_mem structure */
1908 ida_memB = (void *) IDAB_mem->IDA_mem;
1909 return(IDASetLinearSolutionScaling(ida_memB, onoffB));
1910 }
1911
1912
IDASetIncrementFactorB(void * ida_mem,int which,realtype dqincfacB)1913 int IDASetIncrementFactorB(void *ida_mem, int which, realtype dqincfacB)
1914 {
1915 IDAadjMem IDAADJ_mem;
1916 IDAMem IDA_mem;
1917 IDABMem IDAB_mem;
1918 IDALsMemB idalsB_mem;
1919 void *ida_memB;
1920 int retval;
1921
1922 /* access relevant memory structures */
1923 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetIncrementFactorB",
1924 &IDA_mem, &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
1925 if (retval != IDALS_SUCCESS) return(retval);
1926
1927 /* call corresponding routine for IDAB_mem structure */
1928 ida_memB = (void *) IDAB_mem->IDA_mem;
1929 return(IDASetIncrementFactor(ida_memB, dqincfacB));
1930 }
1931
1932
IDASetPreconditionerB(void * ida_mem,int which,IDALsPrecSetupFnB psetupB,IDALsPrecSolveFnB psolveB)1933 int IDASetPreconditionerB(void *ida_mem, int which,
1934 IDALsPrecSetupFnB psetupB,
1935 IDALsPrecSolveFnB psolveB)
1936 {
1937 IDAadjMem IDAADJ_mem;
1938 IDAMem IDA_mem;
1939 IDABMem IDAB_mem;
1940 void *ida_memB;
1941 IDALsMemB idalsB_mem;
1942 IDALsPrecSetupFn idals_psetup;
1943 IDALsPrecSolveFn idals_psolve;
1944 int retval;
1945
1946 /* access relevant memory structures */
1947 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetPreconditionerB",
1948 &IDA_mem, &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
1949 if (retval != IDALS_SUCCESS) return(retval);
1950
1951 /* Set preconditioners for the backward problem. */
1952 idalsB_mem->psetB = psetupB;
1953 idalsB_mem->psolveB = psolveB;
1954
1955 /* Call the corresponding "set" routine for the backward problem */
1956 ida_memB = (void *) IDAB_mem->IDA_mem;
1957 idals_psetup = (psetupB == NULL) ? NULL : idaLsPrecSetupB;
1958 idals_psolve = (psolveB == NULL) ? NULL : idaLsPrecSolveB;
1959 return(IDASetPreconditioner(ida_memB, idals_psetup, idals_psolve));
1960 }
1961
1962
IDASetPreconditionerBS(void * ida_mem,int which,IDALsPrecSetupFnBS psetupBS,IDALsPrecSolveFnBS psolveBS)1963 int IDASetPreconditionerBS(void *ida_mem, int which,
1964 IDALsPrecSetupFnBS psetupBS,
1965 IDALsPrecSolveFnBS psolveBS)
1966 {
1967 IDAadjMem IDAADJ_mem;
1968 IDAMem IDA_mem;
1969 IDABMem IDAB_mem;
1970 void *ida_memB;
1971 IDALsMemB idalsB_mem;
1972 IDALsPrecSetupFn idals_psetup;
1973 IDALsPrecSolveFn idals_psolve;
1974 int retval;
1975
1976 /* access relevant memory structures */
1977 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetPreconditionerBS",
1978 &IDA_mem, &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
1979 if (retval != IDALS_SUCCESS) return(retval);
1980
1981 /* Set preconditioners for the backward problem. */
1982 idalsB_mem->psetBS = psetupBS;
1983 idalsB_mem->psolveBS = psolveBS;
1984
1985 /* Call the corresponding "set" routine for the backward problem */
1986 ida_memB = (void *) IDAB_mem->IDA_mem;
1987 idals_psetup = (psetupBS == NULL) ? NULL : idaLsPrecSetupBS;
1988 idals_psolve = (psolveBS == NULL) ? NULL : idaLsPrecSolveBS;
1989 return(IDASetPreconditioner(ida_memB, idals_psetup, idals_psolve));
1990 }
1991
1992
IDASetJacTimesB(void * ida_mem,int which,IDALsJacTimesSetupFnB jtsetupB,IDALsJacTimesVecFnB jtimesB)1993 int IDASetJacTimesB(void *ida_mem, int which,
1994 IDALsJacTimesSetupFnB jtsetupB,
1995 IDALsJacTimesVecFnB jtimesB)
1996 {
1997 IDAadjMem IDAADJ_mem;
1998 IDAMem IDA_mem;
1999 IDABMem IDAB_mem;
2000 void *ida_memB;
2001 IDALsMemB idalsB_mem;
2002 IDALsJacTimesSetupFn idals_jtsetup;
2003 IDALsJacTimesVecFn idals_jtimes;
2004 int retval;
2005
2006 /* access relevant memory structures */
2007 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetJacTimesB", &IDA_mem,
2008 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2009 if (retval != IDALS_SUCCESS) return(retval);
2010
2011 /* Set jacobian routines for the backward problem. */
2012 idalsB_mem->jtsetupB = jtsetupB;
2013 idalsB_mem->jtimesB = jtimesB;
2014
2015 /* Call the corresponding "set" routine for the backward problem */
2016 ida_memB = (void *) IDAB_mem->IDA_mem;
2017 idals_jtsetup = (jtsetupB == NULL) ? NULL : idaLsJacTimesSetupB;
2018 idals_jtimes = (jtimesB == NULL) ? NULL : idaLsJacTimesVecB;
2019 return(IDASetJacTimes(ida_memB, idals_jtsetup, idals_jtimes));
2020 }
2021
2022
IDASetJacTimesBS(void * ida_mem,int which,IDALsJacTimesSetupFnBS jtsetupBS,IDALsJacTimesVecFnBS jtimesBS)2023 int IDASetJacTimesBS(void *ida_mem, int which,
2024 IDALsJacTimesSetupFnBS jtsetupBS,
2025 IDALsJacTimesVecFnBS jtimesBS)
2026 {
2027 IDAadjMem IDAADJ_mem;
2028 IDAMem IDA_mem;
2029 IDABMem IDAB_mem;
2030 void *ida_memB;
2031 IDALsMemB idalsB_mem;
2032 IDALsJacTimesSetupFn idals_jtsetup;
2033 IDALsJacTimesVecFn idals_jtimes;
2034 int retval;
2035
2036 /* access relevant memory structures */
2037 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetJacTimesBS", &IDA_mem,
2038 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2039 if (retval != IDALS_SUCCESS) return(retval);
2040
2041 /* Set jacobian routines for the backward problem. */
2042 idalsB_mem->jtsetupBS = jtsetupBS;
2043 idalsB_mem->jtimesBS = jtimesBS;
2044
2045 /* Call the corresponding "set" routine for the backward problem */
2046 ida_memB = (void *) IDAB_mem->IDA_mem;
2047 idals_jtsetup = (jtsetupBS == NULL) ? NULL : idaLsJacTimesSetupBS;
2048 idals_jtimes = (jtimesBS == NULL) ? NULL : idaLsJacTimesVecBS;
2049 return(IDASetJacTimes(ida_memB, idals_jtsetup, idals_jtimes));
2050 }
2051
2052
IDASetJacTimesResFnB(void * ida_mem,int which,IDAResFn jtimesResFn)2053 int IDASetJacTimesResFnB(void *ida_mem, int which, IDAResFn jtimesResFn)
2054 {
2055 IDAadjMem IDAADJ_mem;
2056 IDAMem IDA_mem;
2057 IDABMem IDAB_mem;
2058 IDALsMemB idalsB_mem;
2059 void *ida_memB;
2060 int retval;
2061
2062 /* access relevant memory structures */
2063 retval = idaLs_AccessLMemB(ida_mem, which, "IDASetJacTimesResFnB", &IDA_mem,
2064 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2065 if (retval != IDALS_SUCCESS) return(retval);
2066
2067 /* call corresponding routine for IDAB_mem structure */
2068 ida_memB = (void *) IDAB_mem->IDA_mem;
2069 return(IDASetJacTimesResFn(ida_memB, jtimesResFn));
2070 }
2071
2072
2073 /*-----------------------------------------------------------------
2074 IDASLS Private functions for backwards problems
2075 -----------------------------------------------------------------*/
2076
2077 /* idaLsJacBWrapper interfaces to the IDAJacFnB routine provided
2078 by the user. idaLsJacBWrapper is of type IDALsJacFn. */
idaLsJacBWrapper(realtype tt,realtype c_jB,N_Vector yyB,N_Vector ypB,N_Vector rrB,SUNMatrix JacB,void * ida_mem,N_Vector tmp1B,N_Vector tmp2B,N_Vector tmp3B)2079 static int idaLsJacBWrapper(realtype tt, realtype c_jB, N_Vector yyB,
2080 N_Vector ypB, N_Vector rrB, SUNMatrix JacB,
2081 void *ida_mem, N_Vector tmp1B,
2082 N_Vector tmp2B, N_Vector tmp3B)
2083 {
2084 IDAadjMem IDAADJ_mem;
2085 IDAMem IDA_mem;
2086 IDABMem IDAB_mem;
2087 IDALsMemB idalsB_mem;
2088 int retval;
2089
2090 /* access relevant memory structures */
2091 IDA_mem = NULL;
2092 IDAADJ_mem = NULL;
2093 idalsB_mem = NULL;
2094 IDAB_mem = NULL;
2095 retval = idaLs_AccessLMemBCur(ida_mem, "idaLsJacBWrapper", &IDA_mem,
2096 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2097
2098 /* Forward solution from interpolation */
2099 if (IDAADJ_mem->ia_noInterp == SUNFALSE) {
2100 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2101 IDAADJ_mem->ia_ypTmp, NULL, NULL);
2102 if (retval != IDA_SUCCESS) {
2103 IDAProcessError(IDAB_mem->IDA_mem, -1, "IDASLS",
2104 "idaLsJacBWrapper", MSG_LS_BAD_T);
2105 return(-1);
2106 }
2107 }
2108
2109 /* Call user's adjoint jacB routine */
2110 return(idalsB_mem->jacB(tt, c_jB, IDAADJ_mem->ia_yyTmp,
2111 IDAADJ_mem->ia_ypTmp, yyB, ypB,
2112 rrB, JacB, IDAB_mem->ida_user_data,
2113 tmp1B, tmp2B, tmp3B));
2114 }
2115
2116 /* idaLsJacBSWrapper interfaces to the IDAJacFnBS routine provided
2117 by the user. idaLsJacBSWrapper is of type IDALsJacFn. */
idaLsJacBSWrapper(realtype tt,realtype c_jB,N_Vector yyB,N_Vector ypB,N_Vector rrB,SUNMatrix JacB,void * ida_mem,N_Vector tmp1B,N_Vector tmp2B,N_Vector tmp3B)2118 static int idaLsJacBSWrapper(realtype tt, realtype c_jB, N_Vector yyB,
2119 N_Vector ypB, N_Vector rrB, SUNMatrix JacB,
2120 void *ida_mem, N_Vector tmp1B,
2121 N_Vector tmp2B, N_Vector tmp3B)
2122 {
2123 IDAadjMem IDAADJ_mem;
2124 IDAMem IDA_mem;
2125 IDABMem IDAB_mem;
2126 IDALsMemB idalsB_mem;
2127 int retval;
2128
2129 /* access relevant memory structures */
2130 IDA_mem = NULL;
2131 IDAADJ_mem = NULL;
2132 idalsB_mem = NULL;
2133 IDAB_mem = NULL;
2134 retval = idaLs_AccessLMemBCur(ida_mem, "idaLsJacBSWrapper", &IDA_mem,
2135 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2136
2137 /* Get forward solution from interpolation. */
2138 if(IDAADJ_mem->ia_noInterp == SUNFALSE) {
2139 if (IDAADJ_mem->ia_interpSensi)
2140 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2141 IDAADJ_mem->ia_ypTmp, IDAADJ_mem->ia_yySTmp,
2142 IDAADJ_mem->ia_ypSTmp);
2143 else
2144 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2145 IDAADJ_mem->ia_ypTmp, NULL, NULL);
2146
2147 if (retval != IDA_SUCCESS) {
2148 IDAProcessError(IDAB_mem->IDA_mem, -1, "IDASLS",
2149 "idaLsJacBSWrapper", MSG_LS_BAD_T);
2150 return(-1);
2151 }
2152 }
2153
2154 /* Call user's adjoint jacBS routine */
2155 return(idalsB_mem->jacBS(tt, c_jB, IDAADJ_mem->ia_yyTmp,
2156 IDAADJ_mem->ia_ypTmp, IDAADJ_mem->ia_yySTmp,
2157 IDAADJ_mem->ia_ypSTmp, yyB, ypB, rrB, JacB,
2158 IDAB_mem->ida_user_data, tmp1B, tmp2B, tmp3B));
2159 }
2160
2161
2162 /* idaLsPrecSetupB interfaces to the IDALsPrecSetupFnB
2163 routine provided by the user */
idaLsPrecSetupB(realtype tt,N_Vector yyB,N_Vector ypB,N_Vector rrB,realtype c_jB,void * ida_mem)2164 static int idaLsPrecSetupB(realtype tt, N_Vector yyB, N_Vector ypB,
2165 N_Vector rrB, realtype c_jB, void *ida_mem)
2166 {
2167 IDAMem IDA_mem;
2168 IDAadjMem IDAADJ_mem;
2169 IDALsMemB idalsB_mem;
2170 IDABMem IDAB_mem;
2171 int retval;
2172
2173 /* access relevant memory structures */
2174 IDA_mem = NULL;
2175 IDAADJ_mem = NULL;
2176 idalsB_mem = NULL;
2177 IDAB_mem = NULL;
2178 retval = idaLs_AccessLMemBCur(ida_mem, "idaLsPrecSetupB", &IDA_mem,
2179 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2180
2181 /* Get forward solution from interpolation. */
2182 if (IDAADJ_mem->ia_noInterp==SUNFALSE) {
2183 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2184 IDAADJ_mem->ia_ypTmp, NULL, NULL);
2185 if (retval != IDA_SUCCESS) {
2186 IDAProcessError(IDAB_mem->IDA_mem, -1, "IDASLS",
2187 "idaLsPrecSetupB", MSG_LS_BAD_T);
2188 return(-1);
2189 }
2190 }
2191
2192 /* Call user's adjoint precondB routine */
2193 return(idalsB_mem->psetB(tt, IDAADJ_mem->ia_yyTmp,
2194 IDAADJ_mem->ia_ypTmp, yyB, ypB, rrB,
2195 c_jB, IDAB_mem->ida_user_data));
2196 }
2197
2198
2199 /* idaLsPrecSetupBS interfaces to the IDALsPrecSetupFnBS routine
2200 provided by the user */
idaLsPrecSetupBS(realtype tt,N_Vector yyB,N_Vector ypB,N_Vector rrB,realtype c_jB,void * ida_mem)2201 static int idaLsPrecSetupBS(realtype tt, N_Vector yyB, N_Vector ypB,
2202 N_Vector rrB, realtype c_jB, void *ida_mem)
2203 {
2204 IDAMem IDA_mem;
2205 IDAadjMem IDAADJ_mem;
2206 IDALsMemB idalsB_mem;
2207 IDABMem IDAB_mem;
2208 int retval;
2209
2210 /* access relevant memory structures */
2211 IDA_mem = NULL;
2212 IDAADJ_mem = NULL;
2213 idalsB_mem = NULL;
2214 IDAB_mem = NULL;
2215 retval = idaLs_AccessLMemBCur(ida_mem, "idaLsPrecSetupBS", &IDA_mem,
2216 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2217
2218 /* Get forward solution from interpolation. */
2219 if(IDAADJ_mem->ia_noInterp == SUNFALSE) {
2220 if (IDAADJ_mem->ia_interpSensi)
2221 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2222 IDAADJ_mem->ia_ypTmp,
2223 IDAADJ_mem->ia_yySTmp,
2224 IDAADJ_mem->ia_ypSTmp);
2225 else
2226 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2227 IDAADJ_mem->ia_ypTmp, NULL, NULL);
2228 if (retval != IDA_SUCCESS) {
2229 IDAProcessError(IDAB_mem->IDA_mem, -1, "IDASLS",
2230 "idaLsPrecSetupBS", MSG_LS_BAD_T);
2231 return(-1);
2232 }
2233 }
2234
2235 /* Call user's adjoint precondBS routine */
2236 return(idalsB_mem->psetBS(tt, IDAADJ_mem->ia_yyTmp,
2237 IDAADJ_mem->ia_ypTmp,
2238 IDAADJ_mem->ia_yySTmp,
2239 IDAADJ_mem->ia_ypSTmp, yyB, ypB,
2240 rrB, c_jB, IDAB_mem->ida_user_data));
2241 }
2242
2243
2244 /* idaLsPrecSolveB interfaces to the IDALsPrecSolveFnB routine
2245 provided by the user */
idaLsPrecSolveB(realtype tt,N_Vector yyB,N_Vector ypB,N_Vector rrB,N_Vector rvecB,N_Vector zvecB,realtype c_jB,realtype deltaB,void * ida_mem)2246 static int idaLsPrecSolveB(realtype tt, N_Vector yyB, N_Vector ypB,
2247 N_Vector rrB, N_Vector rvecB,
2248 N_Vector zvecB, realtype c_jB,
2249 realtype deltaB, void *ida_mem)
2250 {
2251 IDAMem IDA_mem;
2252 IDAadjMem IDAADJ_mem;
2253 IDALsMemB idalsB_mem;
2254 IDABMem IDAB_mem;
2255 int retval;
2256
2257 /* access relevant memory structures */
2258 IDA_mem = NULL;
2259 IDAADJ_mem = NULL;
2260 idalsB_mem = NULL;
2261 IDAB_mem = NULL;
2262 retval = idaLs_AccessLMemBCur(ida_mem, "idaLsPrecSolveB", &IDA_mem,
2263 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2264
2265 /* Get forward solution from interpolation. */
2266 if (IDAADJ_mem->ia_noInterp==SUNFALSE) {
2267 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2268 IDAADJ_mem->ia_ypTmp, NULL, NULL);
2269 if (retval != IDA_SUCCESS) {
2270 IDAProcessError(IDAB_mem->IDA_mem, -1, "IDASLS",
2271 "idaLsPrecSolveB", MSG_LS_BAD_T);
2272 return(-1);
2273 }
2274 }
2275
2276 /* Call user's adjoint psolveB routine */
2277 return(idalsB_mem->psolveB(tt, IDAADJ_mem->ia_yyTmp,
2278 IDAADJ_mem->ia_ypTmp, yyB, ypB,
2279 rrB, rvecB, zvecB, c_jB, deltaB,
2280 IDAB_mem->ida_user_data));
2281 }
2282
2283
2284 /* idaLsPrecSolveBS interfaces to the IDALsPrecSolveFnBS routine
2285 provided by the user */
idaLsPrecSolveBS(realtype tt,N_Vector yyB,N_Vector ypB,N_Vector rrB,N_Vector rvecB,N_Vector zvecB,realtype c_jB,realtype deltaB,void * ida_mem)2286 static int idaLsPrecSolveBS(realtype tt, N_Vector yyB, N_Vector ypB,
2287 N_Vector rrB, N_Vector rvecB,
2288 N_Vector zvecB, realtype c_jB,
2289 realtype deltaB, void *ida_mem)
2290 {
2291 IDAMem IDA_mem;
2292 IDAadjMem IDAADJ_mem;
2293 IDALsMemB idalsB_mem;
2294 IDABMem IDAB_mem;
2295 int retval;
2296
2297 /* access relevant memory structures */
2298 IDA_mem = NULL;
2299 IDAADJ_mem = NULL;
2300 idalsB_mem = NULL;
2301 IDAB_mem = NULL;
2302 retval = idaLs_AccessLMemBCur(ida_mem, "idaLsPrecSolveBS", &IDA_mem,
2303 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2304
2305 /* Get forward solution from interpolation. */
2306 if(IDAADJ_mem->ia_noInterp == SUNFALSE) {
2307 if (IDAADJ_mem->ia_interpSensi)
2308 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2309 IDAADJ_mem->ia_ypTmp,
2310 IDAADJ_mem->ia_yySTmp,
2311 IDAADJ_mem->ia_ypSTmp);
2312 else
2313 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2314 IDAADJ_mem->ia_ypTmp, NULL, NULL);
2315 if (retval != IDA_SUCCESS) {
2316 IDAProcessError(IDAB_mem->IDA_mem, -1, "IDASLS",
2317 "idaLsPrecSolveBS", MSG_LS_BAD_T);
2318 return(-1);
2319 }
2320 }
2321
2322 /* Call user's adjoint psolveBS routine */
2323 return(idalsB_mem->psolveBS(tt, IDAADJ_mem->ia_yyTmp,
2324 IDAADJ_mem->ia_ypTmp,
2325 IDAADJ_mem->ia_yySTmp,
2326 IDAADJ_mem->ia_ypSTmp,
2327 yyB, ypB, rrB, rvecB, zvecB, c_jB,
2328 deltaB, IDAB_mem->ida_user_data));
2329 }
2330
2331
2332 /* idaLsJacTimesSetupB interfaces to the IDALsJacTimesSetupFnB
2333 routine provided by the user */
idaLsJacTimesSetupB(realtype tt,N_Vector yyB,N_Vector ypB,N_Vector rrB,realtype c_jB,void * ida_mem)2334 static int idaLsJacTimesSetupB(realtype tt, N_Vector yyB, N_Vector ypB,
2335 N_Vector rrB, realtype c_jB, void *ida_mem)
2336 {
2337 IDAMem IDA_mem;
2338 IDAadjMem IDAADJ_mem;
2339 IDALsMemB idalsB_mem;
2340 IDABMem IDAB_mem;
2341 int retval;
2342
2343 /* access relevant memory structures */
2344 IDA_mem = NULL;
2345 IDAADJ_mem = NULL;
2346 idalsB_mem = NULL;
2347 IDAB_mem = NULL;
2348 retval = idaLs_AccessLMemBCur(ida_mem, "idaLsJacTimesSetupB", &IDA_mem,
2349 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2350
2351 /* Get forward solution from interpolation. */
2352 if (IDAADJ_mem->ia_noInterp==SUNFALSE) {
2353 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2354 IDAADJ_mem->ia_ypTmp, NULL, NULL);
2355 if (retval != IDA_SUCCESS) {
2356 IDAProcessError(IDAB_mem->IDA_mem, -1, "IDASLS",
2357 "idaLsJacTimesSetupB", MSG_LS_BAD_T);
2358 return(-1);
2359 }
2360 }
2361 /* Call user's adjoint jtsetupB routine */
2362 return(idalsB_mem->jtsetupB(tt, IDAADJ_mem->ia_yyTmp,
2363 IDAADJ_mem->ia_ypTmp, yyB,
2364 ypB, rrB, c_jB, IDAB_mem->ida_user_data));
2365 }
2366
2367
2368 /* idaLsJacTimesSetupBS interfaces to the IDALsJacTimesSetupFnBS
2369 routine provided by the user */
idaLsJacTimesSetupBS(realtype tt,N_Vector yyB,N_Vector ypB,N_Vector rrB,realtype c_jB,void * ida_mem)2370 static int idaLsJacTimesSetupBS(realtype tt, N_Vector yyB, N_Vector ypB,
2371 N_Vector rrB, realtype c_jB, void *ida_mem)
2372 {
2373 IDAMem IDA_mem;
2374 IDAadjMem IDAADJ_mem;
2375 IDALsMemB idalsB_mem;
2376 IDABMem IDAB_mem;
2377 int retval;
2378
2379 /* access relevant memory structures */
2380 IDA_mem = NULL;
2381 IDAADJ_mem = NULL;
2382 idalsB_mem = NULL;
2383 IDAB_mem = NULL;
2384 retval = idaLs_AccessLMemBCur(ida_mem, "idaLsJacTimesSetupBS", &IDA_mem,
2385 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2386
2387 /* Get forward solution from interpolation. */
2388 if(IDAADJ_mem->ia_noInterp == SUNFALSE) {
2389 if (IDAADJ_mem->ia_interpSensi)
2390 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2391 IDAADJ_mem->ia_ypTmp,
2392 IDAADJ_mem->ia_yySTmp,
2393 IDAADJ_mem->ia_ypSTmp);
2394 else
2395 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2396 IDAADJ_mem->ia_ypTmp, NULL, NULL);
2397 if (retval != IDA_SUCCESS) {
2398 IDAProcessError(IDAB_mem->IDA_mem, -1, "IDASLS",
2399 "idaLsJacTimesSetupBS", MSG_LS_BAD_T);
2400 return(-1);
2401 }
2402 }
2403
2404 /* Call user's adjoint jtimesBS routine */
2405 return(idalsB_mem->jtsetupBS(tt, IDAADJ_mem->ia_yyTmp,
2406 IDAADJ_mem->ia_ypTmp,
2407 IDAADJ_mem->ia_yySTmp,
2408 IDAADJ_mem->ia_ypSTmp,
2409 yyB, ypB, rrB, c_jB,
2410 IDAB_mem->ida_user_data));
2411 }
2412
2413
2414 /* idaLsJacTimesVecB interfaces to the IDALsJacTimesVecFnB routine
2415 provided by the user */
idaLsJacTimesVecB(realtype tt,N_Vector yyB,N_Vector ypB,N_Vector rrB,N_Vector vB,N_Vector JvB,realtype c_jB,void * ida_mem,N_Vector tmp1B,N_Vector tmp2B)2416 static int idaLsJacTimesVecB(realtype tt, N_Vector yyB, N_Vector ypB,
2417 N_Vector rrB, N_Vector vB, N_Vector JvB,
2418 realtype c_jB, void *ida_mem,
2419 N_Vector tmp1B, N_Vector tmp2B)
2420 {
2421 IDAMem IDA_mem;
2422 IDAadjMem IDAADJ_mem;
2423 IDALsMemB idalsB_mem;
2424 IDABMem IDAB_mem;
2425 int retval;
2426
2427 /* access relevant memory structures */
2428 IDA_mem = NULL;
2429 IDAADJ_mem = NULL;
2430 idalsB_mem = NULL;
2431 IDAB_mem = NULL;
2432 retval = idaLs_AccessLMemBCur(ida_mem, "idaLsJacTimesVecB", &IDA_mem,
2433 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2434
2435 /* Get forward solution from interpolation. */
2436 if (IDAADJ_mem->ia_noInterp==SUNFALSE) {
2437 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2438 IDAADJ_mem->ia_ypTmp, NULL, NULL);
2439 if (retval != IDA_SUCCESS) {
2440 IDAProcessError(IDAB_mem->IDA_mem, -1, "IDASLS",
2441 "idaLsJacTimesVecB", MSG_LS_BAD_T);
2442 return(-1);
2443 }
2444 }
2445
2446 /* Call user's adjoint jtimesB routine */
2447 return(idalsB_mem->jtimesB(tt, IDAADJ_mem->ia_yyTmp,
2448 IDAADJ_mem->ia_ypTmp, yyB,
2449 ypB, rrB, vB, JvB, c_jB,
2450 IDAB_mem->ida_user_data,
2451 tmp1B, tmp2B));
2452 }
2453
2454
2455 /* idaLsJacTimesVecBS interfaces to the IDALsJacTimesVecFnBS routine
2456 provided by the user */
idaLsJacTimesVecBS(realtype tt,N_Vector yyB,N_Vector ypB,N_Vector rrB,N_Vector vB,N_Vector JvB,realtype c_jB,void * ida_mem,N_Vector tmp1B,N_Vector tmp2B)2457 static int idaLsJacTimesVecBS(realtype tt, N_Vector yyB, N_Vector ypB,
2458 N_Vector rrB, N_Vector vB, N_Vector JvB,
2459 realtype c_jB, void *ida_mem,
2460 N_Vector tmp1B, N_Vector tmp2B)
2461 {
2462 IDAMem IDA_mem;
2463 IDAadjMem IDAADJ_mem;
2464 IDALsMemB idalsB_mem;
2465 IDABMem IDAB_mem;
2466 int retval;
2467
2468 /* access relevant memory structures */
2469 IDA_mem = NULL;
2470 IDAADJ_mem = NULL;
2471 idalsB_mem = NULL;
2472 IDAB_mem = NULL;
2473 retval = idaLs_AccessLMemBCur(ida_mem, "idaLsJacTimesVecBS", &IDA_mem,
2474 &IDAADJ_mem, &IDAB_mem, &idalsB_mem);
2475
2476 /* Get forward solution from interpolation. */
2477 if(IDAADJ_mem->ia_noInterp == SUNFALSE) {
2478 if (IDAADJ_mem->ia_interpSensi)
2479 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2480 IDAADJ_mem->ia_ypTmp,
2481 IDAADJ_mem->ia_yySTmp,
2482 IDAADJ_mem->ia_ypSTmp);
2483 else
2484 retval = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
2485 IDAADJ_mem->ia_ypTmp, NULL, NULL);
2486 if (retval != IDA_SUCCESS) {
2487 IDAProcessError(IDAB_mem->IDA_mem, -1, "IDASLS",
2488 "idaLsJacTimesVecBS", MSG_LS_BAD_T);
2489 return(-1);
2490 }
2491 }
2492
2493 /* Call user's adjoint jtimesBS routine */
2494 return(idalsB_mem->jtimesBS(tt, IDAADJ_mem->ia_yyTmp,
2495 IDAADJ_mem->ia_ypTmp,
2496 IDAADJ_mem->ia_yySTmp,
2497 IDAADJ_mem->ia_ypSTmp,
2498 yyB, ypB, rrB, vB, JvB, c_jB,
2499 IDAB_mem->ida_user_data, tmp1B, tmp2B));
2500 }
2501
2502
2503 /* idaLsFreeB frees memory associated with the IDASLS wrapper */
idaLsFreeB(IDABMem IDAB_mem)2504 int idaLsFreeB(IDABMem IDAB_mem)
2505 {
2506 IDALsMemB idalsB_mem;
2507
2508 /* Return immediately if IDAB_mem or IDAB_mem->ida_lmem are NULL */
2509 if (IDAB_mem == NULL) return(IDALS_SUCCESS);
2510 if (IDAB_mem->ida_lmem == NULL) return(IDALS_SUCCESS);
2511 idalsB_mem = (IDALsMemB) IDAB_mem->ida_lmem;
2512
2513 /* free IDALsMemB interface structure */
2514 free(idalsB_mem);
2515
2516 return(IDALS_SUCCESS);
2517 }
2518
2519
2520 /* idaLs_AccessLMemB unpacks the IDA_mem, IDAADJ_mem, IDAB_mem and
2521 idalsB_mem structures from the void* ida_mem pointer.
2522 If any are missing it returns IDALS_MEM_NULL, IDALS_NO_ADJ,
2523 IDAS_ILL_INPUT, or IDALS_LMEMB_NULL. */
idaLs_AccessLMemB(void * ida_mem,int which,const char * fname,IDAMem * IDA_mem,IDAadjMem * IDAADJ_mem,IDABMem * IDAB_mem,IDALsMemB * idalsB_mem)2524 int idaLs_AccessLMemB(void *ida_mem, int which, const char *fname,
2525 IDAMem *IDA_mem, IDAadjMem *IDAADJ_mem,
2526 IDABMem *IDAB_mem, IDALsMemB *idalsB_mem)
2527 {
2528
2529 /* access IDAMem structure */
2530 if (ida_mem==NULL) {
2531 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASLS",
2532 fname, MSG_LS_IDAMEM_NULL);
2533 return(IDALS_MEM_NULL);
2534 }
2535 *IDA_mem = (IDAMem) ida_mem;
2536
2537 /* access IDAadjMem structure */
2538 if ((*IDA_mem)->ida_adjMallocDone == SUNFALSE) {
2539 IDAProcessError(*IDA_mem, IDALS_NO_ADJ, "IDASLS",
2540 fname, MSG_LS_NO_ADJ);
2541 return(IDALS_NO_ADJ);
2542 }
2543 *IDAADJ_mem = (*IDA_mem)->ida_adj_mem;
2544
2545 /* Check the value of which */
2546 if ( which >= (*IDAADJ_mem)->ia_nbckpbs ) {
2547 IDAProcessError(*IDA_mem, IDALS_ILL_INPUT, "IDASLS",
2548 fname, MSG_LS_BAD_WHICH);
2549 return(IDALS_ILL_INPUT);
2550 }
2551
2552 /* Find the IDABMem entry in the linked list corresponding to which */
2553 *IDAB_mem = (*IDAADJ_mem)->IDAB_mem;
2554 while ((*IDAB_mem) != NULL) {
2555 if ( which == (*IDAB_mem)->ida_index ) break;
2556 *IDAB_mem = (*IDAB_mem)->ida_next;
2557 }
2558
2559 /* access IDALsMemB structure */
2560 if ((*IDAB_mem)->ida_lmem == NULL) {
2561 IDAProcessError(*IDA_mem, IDALS_LMEMB_NULL, "IDASLS",
2562 fname, MSG_LS_LMEMB_NULL);
2563 return(IDALS_LMEMB_NULL);
2564 }
2565 *idalsB_mem = (IDALsMemB) ((*IDAB_mem)->ida_lmem);
2566
2567 return(IDALS_SUCCESS);
2568 }
2569
2570
2571 /* idaLs_AccessLMemBCur unpacks the ida_mem, ca_mem, idaB_mem and
2572 idalsB_mem structures from the void* ida_mem pointer.
2573 If any are missing it returns IDALS_MEM_NULL, IDALS_NO_ADJ,
2574 or IDALS_LMEMB_NULL. */
idaLs_AccessLMemBCur(void * ida_mem,const char * fname,IDAMem * IDA_mem,IDAadjMem * IDAADJ_mem,IDABMem * IDAB_mem,IDALsMemB * idalsB_mem)2575 int idaLs_AccessLMemBCur(void *ida_mem, const char *fname,
2576 IDAMem *IDA_mem, IDAadjMem *IDAADJ_mem,
2577 IDABMem *IDAB_mem, IDALsMemB *idalsB_mem)
2578 {
2579
2580 /* access IDAMem structure */
2581 if (ida_mem==NULL) {
2582 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASLS",
2583 fname, MSG_LS_IDAMEM_NULL);
2584 return(IDALS_MEM_NULL);
2585 }
2586 *IDA_mem = (IDAMem) ida_mem;
2587
2588 /* access IDAadjMem structure */
2589 if ((*IDA_mem)->ida_adjMallocDone == SUNFALSE) {
2590 IDAProcessError(*IDA_mem, IDALS_NO_ADJ, "IDASLS",
2591 fname, MSG_LS_NO_ADJ);
2592 return(IDALS_NO_ADJ);
2593 }
2594 *IDAADJ_mem = (*IDA_mem)->ida_adj_mem;
2595
2596 /* get current backward problem */
2597 if ((*IDAADJ_mem)->ia_bckpbCrt == NULL) {
2598 IDAProcessError(*IDA_mem, IDALS_LMEMB_NULL, "IDASLS",
2599 fname, MSG_LS_LMEMB_NULL);
2600 return(IDALS_LMEMB_NULL);
2601 }
2602 *IDAB_mem = (*IDAADJ_mem)->ia_bckpbCrt;
2603
2604 /* access IDALsMemB structure */
2605 if ((*IDAB_mem)->ida_lmem == NULL) {
2606 IDAProcessError(*IDA_mem, IDALS_LMEMB_NULL, "IDASLS",
2607 fname, MSG_LS_LMEMB_NULL);
2608 return(IDALS_LMEMB_NULL);
2609 }
2610 *idalsB_mem = (IDALsMemB) ((*IDAB_mem)->ida_lmem);
2611
2612 return(IDALS_SUCCESS);
2613 }
2614
2615
2616 /*---------------------------------------------------------------
2617 EOF
2618 ---------------------------------------------------------------*/
2619