1 /*
2 * -----------------------------------------------------------------
3 * Programmer(s): Daniel R. Reynolds @ SMU
4 * Alan C. Hindmarsh and Radu Serban @ LLNL
5 * -----------------------------------------------------------------
6 * SUNDIALS Copyright Start
7 * Copyright (c) 2002-2019, Lawrence Livermore National Security
8 * and Southern Methodist University.
9 * All rights reserved.
10 *
11 * See the top-level LICENSE and NOTICE files for details.
12 *
13 * SPDX-License-Identifier: BSD-3-Clause
14 * SUNDIALS Copyright End
15 * -----------------------------------------------------------------
16 * This file contains implementations of routines for a
17 * band-block-diagonal preconditioner, i.e. a block-diagonal
18 * matrix with banded blocks, for use with IDA, the IDASLS
19 * linear solver interface.
20 *
21 * NOTE: With only one processor in use, a banded matrix results
22 * rather than a block-diagonal matrix with banded blocks.
23 * Diagonal blocking occurs at the processor level.
24 * -----------------------------------------------------------------
25 */
26
27 #include <stdio.h>
28 #include <stdlib.h>
29
30 #include "idas_impl.h"
31 #include "idas_ls_impl.h"
32 #include "idas_bbdpre_impl.h"
33 #include <sundials/sundials_math.h>
34 #include <nvector/nvector_serial.h>
35
36
37 #define ZERO RCONST(0.0)
38 #define ONE RCONST(1.0)
39 #define TWO RCONST(2.0)
40
41 /* Prototypes of IDABBDPrecSetup and IDABBDPrecSolve */
42 static int IDABBDPrecSetup(realtype tt, N_Vector yy, N_Vector yp,
43 N_Vector rr, realtype c_j, void *prec_data);
44 static int IDABBDPrecSolve(realtype tt, N_Vector yy, N_Vector yp,
45 N_Vector rr, N_Vector rvec, N_Vector zvec,
46 realtype c_j, realtype delta, void *prec_data);
47
48 /* Prototype for IDABBDPrecFree */
49 static int IDABBDPrecFree(IDAMem ida_mem);
50
51 /* Prototype for difference quotient Jacobian calculation routine */
52 static int IBBDDQJac(IBBDPrecData pdata, realtype tt, realtype cj,
53 N_Vector yy, N_Vector yp, N_Vector gref,
54 N_Vector ytemp, N_Vector yptemp, N_Vector gtemp);
55
56 /* Wrapper functions for adjoint code */
57 static int IDAAglocal(sunindextype NlocalB, realtype tt, N_Vector yyB,
58 N_Vector ypB, N_Vector gvalB, void *user_dataB);
59
60 static int IDAAgcomm(sunindextype NlocalB, realtype tt, N_Vector yyB,
61 N_Vector ypB, void *user_dataB);
62
63 /* Prototype for the pfree routine for backward problems. */
64 static int IDABBDPrecFreeB(IDABMem IDAB_mem);
65
66
67 /*================================================================
68 PART I - forward problems
69 ================================================================*/
70
71 /*---------------------------------------------------------------
72 User-Callable Functions: initialization, reinit and free
73 ---------------------------------------------------------------*/
IDABBDPrecInit(void * ida_mem,sunindextype Nlocal,sunindextype mudq,sunindextype mldq,sunindextype mukeep,sunindextype mlkeep,realtype dq_rel_yy,IDABBDLocalFn Gres,IDABBDCommFn Gcomm)74 int IDABBDPrecInit(void *ida_mem, sunindextype Nlocal,
75 sunindextype mudq, sunindextype mldq,
76 sunindextype mukeep, sunindextype mlkeep,
77 realtype dq_rel_yy,
78 IDABBDLocalFn Gres, IDABBDCommFn Gcomm)
79 {
80 IDAMem IDA_mem;
81 IDALsMem idals_mem;
82 IBBDPrecData pdata;
83 sunindextype muk, mlk, storage_mu, lrw1, liw1;
84 long int lrw, liw;
85 int flag;
86
87 if (ida_mem == NULL) {
88 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASBBDPRE",
89 "IDABBDPrecInit", MSGBBD_MEM_NULL);
90 return(IDALS_MEM_NULL);
91 }
92 IDA_mem = (IDAMem) ida_mem;
93
94 /* Test if the LS linear solver interface has been created */
95 if (IDA_mem->ida_lmem == NULL) {
96 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDASBBDPRE",
97 "IDABBDPrecInit", MSGBBD_LMEM_NULL);
98 return(IDALS_LMEM_NULL);
99 }
100 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
101
102 /* Test compatibility of NVECTOR package with the BBD preconditioner */
103 if(IDA_mem->ida_tempv1->ops->nvgetarraypointer == NULL) {
104 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASBBDPRE",
105 "IDABBDPrecInit", MSGBBD_BAD_NVECTOR);
106 return(IDALS_ILL_INPUT);
107 }
108
109 /* Allocate data memory. */
110 pdata = NULL;
111 pdata = (IBBDPrecData) malloc(sizeof *pdata);
112 if (pdata == NULL) {
113 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASBBDPRE",
114 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
115 return(IDALS_MEM_FAIL);
116 }
117
118 /* Set pointers to glocal and gcomm; load half-bandwidths. */
119 pdata->ida_mem = IDA_mem;
120 pdata->glocal = Gres;
121 pdata->gcomm = Gcomm;
122 pdata->mudq = SUNMIN(Nlocal-1, SUNMAX(0, mudq));
123 pdata->mldq = SUNMIN(Nlocal-1, SUNMAX(0, mldq));
124 muk = SUNMIN(Nlocal-1, SUNMAX(0, mukeep));
125 mlk = SUNMIN(Nlocal-1, SUNMAX(0, mlkeep));
126 pdata->mukeep = muk;
127 pdata->mlkeep = mlk;
128
129 /* Set extended upper half-bandwidth for PP (required for pivoting). */
130 storage_mu = SUNMIN(Nlocal-1, muk+mlk);
131
132 /* Allocate memory for preconditioner matrix. */
133 pdata->PP = NULL;
134 pdata->PP = SUNBandMatrixStorage(Nlocal, muk, mlk, storage_mu);
135 if (pdata->PP == NULL) {
136 free(pdata); pdata = NULL;
137 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASBBDPRE",
138 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
139 return(IDALS_MEM_FAIL);
140 }
141
142 /* Allocate memory for temporary N_Vectors */
143 pdata->zlocal = NULL;
144 pdata->zlocal = N_VNewEmpty_Serial(Nlocal);
145 if (pdata->zlocal == NULL) {
146 SUNMatDestroy(pdata->PP);
147 free(pdata); pdata = NULL;
148 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASBBDPRE",
149 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
150 return(IDALS_MEM_FAIL);
151 }
152 pdata->rlocal = NULL;
153 pdata->rlocal = N_VNewEmpty_Serial(Nlocal);
154 if (pdata->rlocal == NULL) {
155 N_VDestroy(pdata->zlocal);
156 SUNMatDestroy(pdata->PP);
157 free(pdata); pdata = NULL;
158 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASBBDPRE",
159 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
160 return(IDALS_MEM_FAIL);
161 }
162 pdata->tempv1 = NULL;
163 pdata->tempv1 = N_VClone(IDA_mem->ida_tempv1);
164 if (pdata->tempv1 == NULL){
165 N_VDestroy(pdata->rlocal);
166 N_VDestroy(pdata->zlocal);
167 SUNMatDestroy(pdata->PP);
168 free(pdata); pdata = NULL;
169 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASBBDPRE",
170 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
171 return(IDALS_MEM_FAIL);
172 }
173 pdata->tempv2 = NULL;
174 pdata->tempv2 = N_VClone(IDA_mem->ida_tempv1);
175 if (pdata->tempv2 == NULL){
176 N_VDestroy(pdata->rlocal);
177 N_VDestroy(pdata->zlocal);
178 N_VDestroy(pdata->tempv1);
179 SUNMatDestroy(pdata->PP);
180 free(pdata); pdata = NULL;
181 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASBBDPRE",
182 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
183 return(IDALS_MEM_FAIL);
184 }
185 pdata->tempv3 = NULL;
186 pdata->tempv3 = N_VClone(IDA_mem->ida_tempv1);
187 if (pdata->tempv3 == NULL){
188 N_VDestroy(pdata->rlocal);
189 N_VDestroy(pdata->zlocal);
190 N_VDestroy(pdata->tempv1);
191 N_VDestroy(pdata->tempv2);
192 SUNMatDestroy(pdata->PP);
193 free(pdata); pdata = NULL;
194 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASBBDPRE",
195 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
196 return(IDALS_MEM_FAIL);
197 }
198 pdata->tempv4 = NULL;
199 pdata->tempv4 = N_VClone(IDA_mem->ida_tempv1);
200 if (pdata->tempv4 == NULL){
201 N_VDestroy(pdata->rlocal);
202 N_VDestroy(pdata->zlocal);
203 N_VDestroy(pdata->tempv1);
204 N_VDestroy(pdata->tempv2);
205 N_VDestroy(pdata->tempv3);
206 SUNMatDestroy(pdata->PP);
207 free(pdata); pdata = NULL;
208 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASBBDPRE",
209 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
210 return(IDALS_MEM_FAIL);
211 }
212
213 /* Allocate memory for banded linear solver */
214 pdata->LS = NULL;
215 pdata->LS = SUNLinSol_Band(pdata->rlocal, pdata->PP);
216 if (pdata->LS == NULL) {
217 N_VDestroy(pdata->zlocal);
218 N_VDestroy(pdata->rlocal);
219 N_VDestroy(pdata->tempv1);
220 N_VDestroy(pdata->tempv2);
221 N_VDestroy(pdata->tempv3);
222 N_VDestroy(pdata->tempv4);
223 SUNMatDestroy(pdata->PP);
224 free(pdata); pdata = NULL;
225 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASBBDPRE",
226 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
227 return(IDALS_MEM_FAIL);
228 }
229
230 /* initialize band linear solver object */
231 flag = SUNLinSolInitialize(pdata->LS);
232 if (flag != SUNLS_SUCCESS) {
233 N_VDestroy(pdata->zlocal);
234 N_VDestroy(pdata->rlocal);
235 N_VDestroy(pdata->tempv1);
236 N_VDestroy(pdata->tempv2);
237 N_VDestroy(pdata->tempv3);
238 N_VDestroy(pdata->tempv4);
239 SUNMatDestroy(pdata->PP);
240 SUNLinSolFree(pdata->LS);
241 free(pdata); pdata = NULL;
242 IDAProcessError(IDA_mem, IDALS_SUNLS_FAIL, "IDASBBDPRE",
243 "IDABBDPrecInit", MSGBBD_SUNLS_FAIL);
244 return(IDALS_SUNLS_FAIL);
245 }
246
247 /* Set rel_yy based on input value dq_rel_yy (0 implies default). */
248 pdata->rel_yy = (dq_rel_yy > ZERO) ?
249 dq_rel_yy : SUNRsqrt(IDA_mem->ida_uround);
250
251 /* Store Nlocal to be used in IDABBDPrecSetup */
252 pdata->n_local = Nlocal;
253
254 /* Set work space sizes and initialize nge. */
255 pdata->rpwsize = 0;
256 pdata->ipwsize = 0;
257 if (IDA_mem->ida_tempv1->ops->nvspace) {
258 N_VSpace(IDA_mem->ida_tempv1, &lrw1, &liw1);
259 pdata->rpwsize += 4*lrw1;
260 pdata->ipwsize += 4*liw1;
261 }
262 if (pdata->rlocal->ops->nvspace) {
263 N_VSpace(pdata->rlocal, &lrw1, &liw1);
264 pdata->rpwsize += 2*lrw1;
265 pdata->ipwsize += 2*liw1;
266 }
267 if (pdata->PP->ops->space) {
268 flag = SUNMatSpace(pdata->PP, &lrw, &liw);
269 pdata->rpwsize += lrw;
270 pdata->ipwsize += liw;
271 }
272 if (pdata->LS->ops->space) {
273 flag = SUNLinSolSpace(pdata->LS, &lrw, &liw);
274 pdata->rpwsize += lrw;
275 pdata->ipwsize += liw;
276 }
277 pdata->nge = 0;
278
279 /* make sure pdata is free from any previous allocations */
280 if (idals_mem->pfree)
281 idals_mem->pfree(IDA_mem);
282
283 /* Point to the new pdata field in the LS memory */
284 idals_mem->pdata = pdata;
285
286 /* Attach the pfree function */
287 idals_mem->pfree = IDABBDPrecFree;
288
289 /* Attach preconditioner solve and setup functions */
290 flag = IDASetPreconditioner(ida_mem, IDABBDPrecSetup,
291 IDABBDPrecSolve);
292
293 return(flag);
294 }
295
296
297 /*-------------------------------------------------------------*/
IDABBDPrecReInit(void * ida_mem,sunindextype mudq,sunindextype mldq,realtype dq_rel_yy)298 int IDABBDPrecReInit(void *ida_mem, sunindextype mudq,
299 sunindextype mldq, realtype dq_rel_yy)
300 {
301 IDAMem IDA_mem;
302 IDALsMem idals_mem;
303 IBBDPrecData pdata;
304 sunindextype Nlocal;
305
306 if (ida_mem == NULL) {
307 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASBBDPRE",
308 "IDABBDPrecReInit", MSGBBD_MEM_NULL);
309 return(IDALS_MEM_NULL);
310 }
311 IDA_mem = (IDAMem) ida_mem;
312
313 /* Test if the LS linear solver interface has been created */
314 if (IDA_mem->ida_lmem == NULL) {
315 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDASBBDPRE",
316 "IDABBDPrecReInit", MSGBBD_LMEM_NULL);
317 return(IDALS_LMEM_NULL);
318 }
319 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
320
321 /* Test if the preconditioner data is non-NULL */
322 if (idals_mem->pdata == NULL) {
323 IDAProcessError(IDA_mem, IDALS_PMEM_NULL, "IDASBBDPRE",
324 "IDABBDPrecReInit", MSGBBD_PMEM_NULL);
325 return(IDALS_PMEM_NULL);
326 }
327 pdata = (IBBDPrecData) idals_mem->pdata;
328
329 /* Load half-bandwidths. */
330 Nlocal = pdata->n_local;
331 pdata->mudq = SUNMIN(Nlocal-1, SUNMAX(0, mudq));
332 pdata->mldq = SUNMIN(Nlocal-1, SUNMAX(0, mldq));
333
334 /* Set rel_yy based on input value dq_rel_yy (0 implies default). */
335 pdata->rel_yy = (dq_rel_yy > ZERO) ?
336 dq_rel_yy : SUNRsqrt(IDA_mem->ida_uround);
337
338 /* Re-initialize nge */
339 pdata->nge = 0;
340
341 return(IDALS_SUCCESS);
342 }
343
344
345 /*-------------------------------------------------------------*/
IDABBDPrecGetWorkSpace(void * ida_mem,long int * lenrwBBDP,long int * leniwBBDP)346 int IDABBDPrecGetWorkSpace(void *ida_mem,
347 long int *lenrwBBDP,
348 long int *leniwBBDP)
349 {
350 IDAMem IDA_mem;
351 IDALsMem idals_mem;
352 IBBDPrecData pdata;
353
354 if (ida_mem == NULL) {
355 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASBBDPRE",
356 "IDABBDPrecGetWorkSpace", MSGBBD_MEM_NULL);
357 return(IDALS_MEM_NULL);
358 }
359 IDA_mem = (IDAMem) ida_mem;
360
361 if (IDA_mem->ida_lmem == NULL) {
362 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDASBBDPRE",
363 "IDABBDPrecGetWorkSpace", MSGBBD_LMEM_NULL);
364 return(IDALS_LMEM_NULL);
365 }
366 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
367
368 if (idals_mem->pdata == NULL) {
369 IDAProcessError(IDA_mem, IDALS_PMEM_NULL, "IDASBBDPRE",
370 "IDABBDPrecGetWorkSpace", MSGBBD_PMEM_NULL);
371 return(IDALS_PMEM_NULL);
372 }
373 pdata = (IBBDPrecData) idals_mem->pdata;
374
375 *lenrwBBDP = pdata->rpwsize;
376 *leniwBBDP = pdata->ipwsize;
377
378 return(IDALS_SUCCESS);
379 }
380
381
382 /*-------------------------------------------------------------*/
IDABBDPrecGetNumGfnEvals(void * ida_mem,long int * ngevalsBBDP)383 int IDABBDPrecGetNumGfnEvals(void *ida_mem,
384 long int *ngevalsBBDP)
385 {
386 IDAMem IDA_mem;
387 IDALsMem idals_mem;
388 IBBDPrecData pdata;
389
390 if (ida_mem == NULL) {
391 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASBBDPRE",
392 "IDABBDPrecGetNumGfnEvals", MSGBBD_MEM_NULL);
393 return(IDALS_MEM_NULL);
394 }
395 IDA_mem = (IDAMem) ida_mem;
396
397 if (IDA_mem->ida_lmem == NULL) {
398 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDASBBDPRE",
399 "IDABBDPrecGetNumGfnEvals", MSGBBD_LMEM_NULL);
400 return(IDALS_LMEM_NULL);
401 }
402 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
403
404 if (idals_mem->pdata == NULL) {
405 IDAProcessError(IDA_mem, IDALS_PMEM_NULL, "IDASBBDPRE",
406 "IDABBDPrecGetNumGfnEvals", MSGBBD_PMEM_NULL);
407 return(IDALS_PMEM_NULL);
408 }
409 pdata = (IBBDPrecData) idals_mem->pdata;
410
411 *ngevalsBBDP = pdata->nge;
412
413 return(IDALS_SUCCESS);
414 }
415
416
417
418
419 /*---------------------------------------------------------------
420 IDABBDPrecSetup:
421
422 IDABBDPrecSetup generates a band-block-diagonal preconditioner
423 matrix, where the local block (on this processor) is a band
424 matrix. Each local block is computed by a difference quotient
425 scheme via calls to the user-supplied routines glocal, gcomm.
426 After generating the block in the band matrix PP, this routine
427 does an LU factorization in place in PP.
428
429 The IDABBDPrecSetup parameters used here are as follows:
430
431 tt is the current value of the independent variable t.
432
433 yy is the current value of the dependent variable vector,
434 namely the predicted value of y(t).
435
436 yp is the current value of the derivative vector y',
437 namely the predicted value of y'(t).
438
439 c_j is the scalar in the system Jacobian, proportional to 1/hh.
440
441 bbd_data is the pointer to BBD memory set by IDABBDInit
442
443 The argument rr is not used.
444
445 Return value:
446 The value returned by this IDABBDPrecSetup function is a int
447 flag indicating whether it was successful. This value is
448 0 if successful,
449 > 0 for a recoverable error (step will be retried), or
450 < 0 for a nonrecoverable error (step fails).
451 ----------------------------------------------------------------*/
IDABBDPrecSetup(realtype tt,N_Vector yy,N_Vector yp,N_Vector rr,realtype c_j,void * bbd_data)452 static int IDABBDPrecSetup(realtype tt, N_Vector yy, N_Vector yp,
453 N_Vector rr, realtype c_j, void *bbd_data)
454 {
455 sunindextype ier;
456 IBBDPrecData pdata;
457 IDAMem IDA_mem;
458 int retval;
459
460 pdata =(IBBDPrecData) bbd_data;
461
462 IDA_mem = (IDAMem) pdata->ida_mem;
463
464 /* Call IBBDDQJac for a new Jacobian calculation and store in PP. */
465 retval = SUNMatZero(pdata->PP);
466 retval = IBBDDQJac(pdata, tt, c_j, yy, yp, pdata->tempv1,
467 pdata->tempv2, pdata->tempv3, pdata->tempv4);
468 if (retval < 0) {
469 IDAProcessError(IDA_mem, -1, "IDASBBDPRE", "IDABBDPrecSetup",
470 MSGBBD_FUNC_FAILED);
471 return(-1);
472 }
473 if (retval > 0) {
474 return(+1);
475 }
476
477 /* Do LU factorization of matrix and return error flag */
478 ier = SUNLinSolSetup_Band(pdata->LS, pdata->PP);
479 return(ier);
480 }
481
482
483 /*---------------------------------------------------------------
484 IDABBDPrecSolve
485
486 The function IDABBDPrecSolve computes a solution to the linear
487 system P z = r, where P is the left preconditioner defined by
488 the routine IDABBDPrecSetup.
489
490 The IDABBDPrecSolve parameters used here are as follows:
491
492 rvec is the input right-hand side vector r.
493
494 zvec is the computed solution vector z.
495
496 bbd_data is the pointer to BBD data set by IDABBDInit.
497
498 The arguments tt, yy, yp, rr, c_j and delta are NOT used.
499
500 IDABBDPrecSolve returns the value returned from the linear
501 solver object.
502 ---------------------------------------------------------------*/
IDABBDPrecSolve(realtype tt,N_Vector yy,N_Vector yp,N_Vector rr,N_Vector rvec,N_Vector zvec,realtype c_j,realtype delta,void * bbd_data)503 static int IDABBDPrecSolve(realtype tt, N_Vector yy, N_Vector yp,
504 N_Vector rr, N_Vector rvec, N_Vector zvec,
505 realtype c_j, realtype delta, void *bbd_data)
506 {
507 IBBDPrecData pdata;
508 int retval;
509
510 pdata = (IBBDPrecData) bbd_data;
511
512 /* Attach local data arrays for rvec and zvec to rlocal and zlocal */
513 N_VSetArrayPointer(N_VGetArrayPointer(rvec), pdata->rlocal);
514 N_VSetArrayPointer(N_VGetArrayPointer(zvec), pdata->zlocal);
515
516 /* Call banded solver object to do the work */
517 retval = SUNLinSolSolve(pdata->LS, pdata->PP, pdata->zlocal,
518 pdata->rlocal, ZERO);
519
520 /* Detach local data arrays from rlocal and zlocal */
521 N_VSetArrayPointer(NULL, pdata->rlocal);
522 N_VSetArrayPointer(NULL, pdata->zlocal);
523
524 return(retval);
525 }
526
527
528 /*-------------------------------------------------------------*/
IDABBDPrecFree(IDAMem IDA_mem)529 static int IDABBDPrecFree(IDAMem IDA_mem)
530 {
531 IDALsMem idals_mem;
532 IBBDPrecData pdata;
533
534 if (IDA_mem->ida_lmem == NULL) return(0);
535 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
536
537 if (idals_mem->pdata == NULL) return(0);
538 pdata = (IBBDPrecData) idals_mem->pdata;
539
540 SUNLinSolFree(pdata->LS);
541 N_VDestroy(pdata->rlocal);
542 N_VDestroy(pdata->zlocal);
543 N_VDestroy(pdata->tempv1);
544 N_VDestroy(pdata->tempv2);
545 N_VDestroy(pdata->tempv3);
546 N_VDestroy(pdata->tempv4);
547 SUNMatDestroy(pdata->PP);
548
549 free(pdata);
550 pdata = NULL;
551
552 return(0);
553 }
554
555
556 /*---------------------------------------------------------------
557 IBBDDQJac
558
559 This routine generates a banded difference quotient approximation
560 to the local block of the Jacobian of G(t,y,y'). It assumes that
561 a band matrix of type SUNMatrix is stored column-wise, and that
562 elements within each column are contiguous.
563
564 All matrix elements are generated as difference quotients, by way
565 of calls to the user routine glocal. By virtue of the band
566 structure, the number of these calls is bandwidth + 1, where
567 bandwidth = mldq + mudq + 1. But the band matrix kept has
568 bandwidth = mlkeep + mukeep + 1. This routine also assumes that
569 the local elements of a vector are stored contiguously.
570
571 Return values are: 0 (success), > 0 (recoverable error),
572 or < 0 (nonrecoverable error).
573 ----------------------------------------------------------------*/
IBBDDQJac(IBBDPrecData pdata,realtype tt,realtype cj,N_Vector yy,N_Vector yp,N_Vector gref,N_Vector ytemp,N_Vector yptemp,N_Vector gtemp)574 static int IBBDDQJac(IBBDPrecData pdata, realtype tt, realtype cj,
575 N_Vector yy, N_Vector yp, N_Vector gref,
576 N_Vector ytemp, N_Vector yptemp, N_Vector gtemp)
577 {
578 IDAMem IDA_mem;
579 realtype inc, inc_inv;
580 int retval;
581 sunindextype group, i, j, width, ngroups, i1, i2;
582 realtype *ydata, *ypdata, *ytempdata, *yptempdata, *grefdata, *gtempdata;
583 realtype *cnsdata = NULL, *ewtdata;
584 realtype *col_j, conj, yj, ypj, ewtj;
585
586 IDA_mem = (IDAMem) pdata->ida_mem;
587
588 /* Initialize ytemp and yptemp. */
589 N_VScale(ONE, yy, ytemp);
590 N_VScale(ONE, yp, yptemp);
591
592 /* Obtain pointers as required to the data array of vectors. */
593 ydata = N_VGetArrayPointer(yy);
594 ypdata = N_VGetArrayPointer(yp);
595 gtempdata = N_VGetArrayPointer(gtemp);
596 ewtdata = N_VGetArrayPointer(IDA_mem->ida_ewt);
597 if (IDA_mem->ida_constraints != NULL)
598 cnsdata = N_VGetArrayPointer(IDA_mem->ida_constraints);
599 ytempdata = N_VGetArrayPointer(ytemp);
600 yptempdata= N_VGetArrayPointer(yptemp);
601 grefdata = N_VGetArrayPointer(gref);
602
603 /* Call gcomm and glocal to get base value of G(t,y,y'). */
604 if (pdata->gcomm != NULL) {
605 retval = pdata->gcomm(pdata->n_local, tt, yy, yp, IDA_mem->ida_user_data);
606 if (retval != 0) return(retval);
607 }
608
609 retval = pdata->glocal(pdata->n_local, tt, yy, yp, gref, IDA_mem->ida_user_data);
610 pdata->nge++;
611 if (retval != 0) return(retval);
612
613 /* Set bandwidth and number of column groups for band differencing. */
614 width = pdata->mldq + pdata->mudq + 1;
615 ngroups = SUNMIN(width, pdata->n_local);
616
617 /* Loop over groups. */
618 for(group = 1; group <= ngroups; group++) {
619
620 /* Loop over the components in this group. */
621 for(j = group-1; j < pdata->n_local; j += width) {
622 yj = ydata[j];
623 ypj = ypdata[j];
624 ewtj = ewtdata[j];
625
626 /* Set increment inc to yj based on rel_yy*abs(yj), with
627 adjustments using ypj and ewtj if this is small, and a further
628 adjustment to give it the same sign as hh*ypj. */
629 inc = pdata->rel_yy *
630 SUNMAX(SUNRabs(yj), SUNMAX( SUNRabs(IDA_mem->ida_hh*ypj), ONE/ewtj));
631 if (IDA_mem->ida_hh*ypj < ZERO) inc = -inc;
632 inc = (yj + inc) - yj;
633
634 /* Adjust sign(inc) again if yj has an inequality constraint. */
635 if (IDA_mem->ida_constraints != NULL) {
636 conj = cnsdata[j];
637 if (SUNRabs(conj) == ONE) {if ((yj+inc)*conj < ZERO) inc = -inc;}
638 else if (SUNRabs(conj) == TWO) {if ((yj+inc)*conj <= ZERO) inc = -inc;}
639 }
640
641 /* Increment yj and ypj. */
642 ytempdata[j] += inc;
643 yptempdata[j] += cj*inc;
644
645 }
646
647 /* Evaluate G with incremented y and yp arguments. */
648 retval = pdata->glocal(pdata->n_local, tt, ytemp, yptemp,
649 gtemp, IDA_mem->ida_user_data);
650 pdata->nge++;
651 if (retval != 0) return(retval);
652
653 /* Loop over components of the group again; restore ytemp and yptemp. */
654 for(j = group-1; j < pdata->n_local; j += width) {
655 yj = ytempdata[j] = ydata[j];
656 ypj = yptempdata[j] = ypdata[j];
657 ewtj = ewtdata[j];
658
659 /* Set increment inc as before .*/
660 inc = pdata->rel_yy *
661 SUNMAX(SUNRabs(yj), SUNMAX( SUNRabs(IDA_mem->ida_hh*ypj), ONE/ewtj));
662 if (IDA_mem->ida_hh*ypj < ZERO) inc = -inc;
663 inc = (yj + inc) - yj;
664 if (IDA_mem->ida_constraints != NULL) {
665 conj = cnsdata[j];
666 if (SUNRabs(conj) == ONE) {if ((yj+inc)*conj < ZERO) inc = -inc;}
667 else if (SUNRabs(conj) == TWO) {if ((yj+inc)*conj <= ZERO) inc = -inc;}
668 }
669
670 /* Form difference quotients and load into PP. */
671 inc_inv = ONE/inc;
672 col_j = SUNBandMatrix_Column(pdata->PP,j);
673 i1 = SUNMAX(0, j-pdata->mukeep);
674 i2 = SUNMIN(j + pdata->mlkeep, pdata->n_local-1);
675 for(i = i1; i <= i2; i++)
676 SM_COLUMN_ELEMENT_B(col_j,i,j) =
677 inc_inv * (gtempdata[i] - grefdata[i]);
678 }
679 }
680
681 return(0);
682 }
683
684
685 /*================================================================
686 PART II - backward problems
687 ================================================================*/
688
689 /*---------------------------------------------------------------
690 User-Callable Functions: initialization, reinit and free
691 ---------------------------------------------------------------*/
IDABBDPrecInitB(void * ida_mem,int which,sunindextype NlocalB,sunindextype mudqB,sunindextype mldqB,sunindextype mukeepB,sunindextype mlkeepB,realtype dq_rel_yyB,IDABBDLocalFnB glocalB,IDABBDCommFnB gcommB)692 int IDABBDPrecInitB(void *ida_mem, int which, sunindextype NlocalB,
693 sunindextype mudqB, sunindextype mldqB,
694 sunindextype mukeepB, sunindextype mlkeepB,
695 realtype dq_rel_yyB, IDABBDLocalFnB glocalB,
696 IDABBDCommFnB gcommB)
697 {
698 IDAMem IDA_mem;
699 IDAadjMem IDAADJ_mem;
700 IDABMem IDAB_mem;
701 IDABBDPrecDataB idabbdB_mem;
702 void *ida_memB;
703 int flag;
704
705 /* Check if ida_mem is allright. */
706 if (ida_mem == NULL) {
707 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASBBDPRE",
708 "IDABBDPrecInitB", MSG_LS_IDAMEM_NULL);
709 return(IDALS_MEM_NULL);
710 }
711 IDA_mem = (IDAMem) ida_mem;
712
713 /* Is ASA initialized? */
714 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
715 IDAProcessError(IDA_mem, IDALS_NO_ADJ, "IDASBBDPRE",
716 "IDABBDPrecInitB", MSG_LS_NO_ADJ);
717 return(IDALS_NO_ADJ);
718 }
719 IDAADJ_mem = IDA_mem->ida_adj_mem;
720
721 /* Check the value of which */
722 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
723 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASBBDPRE",
724 "IDABBDPrecInitB", MSG_LS_BAD_WHICH);
725 return(IDALS_ILL_INPUT);
726 }
727
728 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
729 IDAB_mem = IDAADJ_mem->IDAB_mem;
730 while (IDAB_mem != NULL) {
731 if( which == IDAB_mem->ida_index ) break;
732 /* advance */
733 IDAB_mem = IDAB_mem->ida_next;
734 }
735 /* ida_mem corresponding to 'which' problem. */
736 ida_memB = (void *) IDAB_mem->IDA_mem;
737
738 /* Initialize the BBD preconditioner for this backward problem. */
739 flag = IDABBDPrecInit(ida_memB, NlocalB, mudqB, mldqB, mukeepB,
740 mlkeepB, dq_rel_yyB, IDAAglocal, IDAAgcomm);
741 if (flag != IDA_SUCCESS) return(flag);
742
743 /* Allocate memory for IDABBDPrecDataB to store the user-provided
744 functions which will be called from the wrappers */
745 idabbdB_mem = NULL;
746 idabbdB_mem = (IDABBDPrecDataB) malloc(sizeof(* idabbdB_mem));
747 if (idabbdB_mem == NULL) {
748 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDASBBDPRE",
749 "IDABBDPrecInitB", MSGBBD_MEM_FAIL);
750 return(IDALS_MEM_FAIL);
751 }
752
753 /* set pointers to user-provided functions */
754 idabbdB_mem->glocalB = glocalB;
755 idabbdB_mem->gcommB = gcommB;
756
757 /* Attach pmem and pfree */
758 IDAB_mem->ida_pmem = idabbdB_mem;
759 IDAB_mem->ida_pfree = IDABBDPrecFreeB;
760
761 return(IDALS_SUCCESS);
762 }
763
764
765 /*-------------------------------------------------------------*/
IDABBDPrecReInitB(void * ida_mem,int which,sunindextype mudqB,sunindextype mldqB,realtype dq_rel_yyB)766 int IDABBDPrecReInitB(void *ida_mem, int which, sunindextype mudqB,
767 sunindextype mldqB, realtype dq_rel_yyB)
768 {
769 IDAMem IDA_mem;
770 IDAadjMem IDAADJ_mem;
771 IDABMem IDAB_mem;
772 void *ida_memB;
773 int flag;
774
775 /* Check if ida_mem is allright. */
776 if (ida_mem == NULL) {
777 IDAProcessError(NULL, IDALS_MEM_NULL, "IDASBBDPRE",
778 "IDABBDPrecReInitB", MSG_LS_IDAMEM_NULL);
779 return(IDALS_MEM_NULL);
780 }
781 IDA_mem = (IDAMem) ida_mem;
782
783 /* Is ASA initialized? */
784 if (IDA_mem->ida_adjMallocDone == SUNFALSE) {
785 IDAProcessError(IDA_mem, IDALS_NO_ADJ, "IDASBBDPRE",
786 "IDABBDPrecReInitB", MSG_LS_NO_ADJ);
787 return(IDALS_NO_ADJ);
788 }
789 IDAADJ_mem = IDA_mem->ida_adj_mem;
790
791 /* Check the value of which */
792 if ( which >= IDAADJ_mem->ia_nbckpbs ) {
793 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDASBBDPRE",
794 "IDABBDPrecReInitB", MSG_LS_BAD_WHICH);
795 return(IDALS_ILL_INPUT);
796 }
797
798 /* Find the IDABMem entry in the linked list corresponding to 'which'. */
799 IDAB_mem = IDAADJ_mem->IDAB_mem;
800 while (IDAB_mem != NULL) {
801 if( which == IDAB_mem->ida_index ) break;
802 /* advance */
803 IDAB_mem = IDAB_mem->ida_next;
804 }
805 /* ida_mem corresponding to 'which' backward problem. */
806 ida_memB = (void *) IDAB_mem->IDA_mem;
807
808 /* ReInitialize the BBD preconditioner for this backward problem. */
809 flag = IDABBDPrecReInit(ida_memB, mudqB, mldqB, dq_rel_yyB);
810 return(flag);
811 }
812
813
814 /*-------------------------------------------------------------*/
IDABBDPrecFreeB(IDABMem IDAB_mem)815 static int IDABBDPrecFreeB(IDABMem IDAB_mem)
816 {
817 free(IDAB_mem->ida_pmem);
818 IDAB_mem->ida_pmem = NULL;
819 return(0);
820 }
821
822
823 /*----------------------------------------------------------------
824 Wrapper functions
825 ----------------------------------------------------------------*/
826
827 /*----------------------------------------------------------------
828 IDAAglocal
829
830 This routine interfaces to the IDALocalFnB routine
831 provided by the user.
832 ----------------------------------------------------------------*/
IDAAglocal(sunindextype NlocalB,realtype tt,N_Vector yyB,N_Vector ypB,N_Vector gvalB,void * ida_mem)833 static int IDAAglocal(sunindextype NlocalB, realtype tt, N_Vector yyB,
834 N_Vector ypB, N_Vector gvalB, void *ida_mem)
835 {
836 IDAMem IDA_mem;
837 IDAadjMem IDAADJ_mem;
838 IDABMem IDAB_mem;
839 IDABBDPrecDataB idabbdB_mem;
840 int flag;
841
842 IDA_mem = (IDAMem) ida_mem;
843 IDAADJ_mem = IDA_mem->ida_adj_mem;
844
845 /* Get current backward problem. */
846 IDAB_mem = IDAADJ_mem->ia_bckpbCrt;
847
848 /* Get the preconditioner's memory. */
849 idabbdB_mem = (IDABBDPrecDataB) IDAB_mem->ida_pmem;
850
851 /* Get forward solution from interpolation. */
852 if (IDAADJ_mem->ia_noInterp == SUNFALSE) {
853 flag = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
854 IDAADJ_mem->ia_ypTmp, NULL, NULL);
855 if (flag != IDA_SUCCESS) {
856 IDAProcessError(IDA_mem, -1, "IDASBBDPRE", "IDAAglocal",
857 MSGBBD_BAD_T);
858 return(-1);
859 }
860 }
861 /* Call user's adjoint LocalFnB function. */
862 return idabbdB_mem->glocalB(NlocalB, tt, IDAADJ_mem->ia_yyTmp,
863 IDAADJ_mem->ia_ypTmp, yyB, ypB,
864 gvalB, IDAB_mem->ida_user_data);
865 }
866
867
868 /*----------------------------------------------------------------
869 IDAAgcomm
870
871 This routine interfaces to the IDACommFnB routine
872 provided by the user.
873 ----------------------------------------------------------------*/
IDAAgcomm(sunindextype NlocalB,realtype tt,N_Vector yyB,N_Vector ypB,void * ida_mem)874 static int IDAAgcomm(sunindextype NlocalB, realtype tt,
875 N_Vector yyB, N_Vector ypB, void *ida_mem)
876 {
877 IDAMem IDA_mem;
878 IDAadjMem IDAADJ_mem;
879 IDABMem IDAB_mem;
880 IDABBDPrecDataB idabbdB_mem;
881 int flag;
882
883 IDA_mem = (IDAMem) ida_mem;
884 IDAADJ_mem = IDA_mem->ida_adj_mem;
885
886 /* Get current backward problem. */
887 IDAB_mem = IDAADJ_mem->ia_bckpbCrt;
888
889 /* Get the preconditioner's memory. */
890 idabbdB_mem = (IDABBDPrecDataB) IDAB_mem->ida_pmem;
891 if (idabbdB_mem->gcommB == NULL) return(0);
892
893 /* Get forward solution from interpolation. */
894 if (IDAADJ_mem->ia_noInterp == SUNFALSE) {
895 flag = IDAADJ_mem->ia_getY(IDA_mem, tt, IDAADJ_mem->ia_yyTmp,
896 IDAADJ_mem->ia_ypTmp, NULL, NULL);
897 if (flag != IDA_SUCCESS) {
898 IDAProcessError(IDA_mem, -1, "IDASBBDPRE", "IDAAgcomm",
899 MSGBBD_BAD_T);
900 return(-1);
901 }
902 }
903
904 /* Call user's adjoint CommFnB routine */
905 return idabbdB_mem->gcommB(NlocalB, tt, IDAADJ_mem->ia_yyTmp,
906 IDAADJ_mem->ia_ypTmp, yyB, ypB,
907 IDAB_mem->ida_user_data);
908 }
909