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-2021, 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 IDASPILS
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 "ida_impl.h"
31 #include "ida_ls_impl.h"
32 #include "ida_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 functions 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 /*---------------------------------------------------------------
57 User-Callable Functions: initialization, reinit and free
58 ---------------------------------------------------------------*/
IDABBDPrecInit(void * ida_mem,sunindextype Nlocal,sunindextype mudq,sunindextype mldq,sunindextype mukeep,sunindextype mlkeep,realtype dq_rel_yy,IDABBDLocalFn Gres,IDABBDCommFn Gcomm)59 int IDABBDPrecInit(void *ida_mem, sunindextype Nlocal,
60 sunindextype mudq, sunindextype mldq,
61 sunindextype mukeep, sunindextype mlkeep,
62 realtype dq_rel_yy,
63 IDABBDLocalFn Gres, IDABBDCommFn Gcomm)
64 {
65 IDAMem IDA_mem;
66 IDALsMem idals_mem;
67 IBBDPrecData pdata;
68 sunindextype muk, mlk, storage_mu, lrw1, liw1;
69 long int lrw, liw;
70 int flag;
71
72 if (ida_mem == NULL) {
73 IDAProcessError(NULL, IDALS_MEM_NULL, "IDABBDPRE",
74 "IDABBDPrecInit", MSGBBD_MEM_NULL);
75 return(IDALS_MEM_NULL);
76 }
77 IDA_mem = (IDAMem) ida_mem;
78
79 /* Test if the LS linear solver interface has been created */
80 if (IDA_mem->ida_lmem == NULL) {
81 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDABBDPRE",
82 "IDABBDPrecInit", MSGBBD_LMEM_NULL);
83 return(IDALS_LMEM_NULL);
84 }
85 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
86
87 /* Test compatibility of NVECTOR package with the BBD preconditioner */
88 if(IDA_mem->ida_tempv1->ops->nvgetarraypointer == NULL) {
89 IDAProcessError(IDA_mem, IDALS_ILL_INPUT, "IDABBDPRE",
90 "IDABBDPrecInit", MSGBBD_BAD_NVECTOR);
91 return(IDALS_ILL_INPUT);
92 }
93
94 /* Allocate data memory. */
95 pdata = NULL;
96 pdata = (IBBDPrecData) malloc(sizeof *pdata);
97 if (pdata == NULL) {
98 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDABBDPRE",
99 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
100 return(IDALS_MEM_FAIL);
101 }
102
103 /* Set pointers to glocal and gcomm; load half-bandwidths. */
104 pdata->ida_mem = IDA_mem;
105 pdata->glocal = Gres;
106 pdata->gcomm = Gcomm;
107 pdata->mudq = SUNMIN(Nlocal-1, SUNMAX(0, mudq));
108 pdata->mldq = SUNMIN(Nlocal-1, SUNMAX(0, mldq));
109 muk = SUNMIN(Nlocal-1, SUNMAX(0, mukeep));
110 mlk = SUNMIN(Nlocal-1, SUNMAX(0, mlkeep));
111 pdata->mukeep = muk;
112 pdata->mlkeep = mlk;
113
114 /* Set extended upper half-bandwidth for PP (required for pivoting). */
115 storage_mu = SUNMIN(Nlocal-1, muk+mlk);
116
117 /* Allocate memory for preconditioner matrix. */
118 pdata->PP = NULL;
119 pdata->PP = SUNBandMatrixStorage(Nlocal, muk, mlk, storage_mu);
120 if (pdata->PP == NULL) {
121 free(pdata); pdata = NULL;
122 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDABBDPRE",
123 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
124 return(IDALS_MEM_FAIL);
125 }
126
127 /* Allocate memory for temporary N_Vectors */
128 pdata->zlocal = NULL;
129 pdata->zlocal = N_VNewEmpty_Serial(Nlocal);
130 if (pdata->zlocal == NULL) {
131 SUNMatDestroy(pdata->PP);
132 free(pdata); pdata = NULL;
133 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDABBDPRE",
134 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
135 return(IDALS_MEM_FAIL);
136 }
137 pdata->rlocal = NULL;
138 pdata->rlocal = N_VNewEmpty_Serial(Nlocal);
139 if (pdata->rlocal == NULL) {
140 N_VDestroy(pdata->zlocal);
141 SUNMatDestroy(pdata->PP);
142 free(pdata); pdata = NULL;
143 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDABBDPRE",
144 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
145 return(IDALS_MEM_FAIL);
146 }
147 pdata->tempv1 = NULL;
148 pdata->tempv1 = N_VClone(IDA_mem->ida_tempv1);
149 if (pdata->tempv1 == NULL){
150 N_VDestroy(pdata->rlocal);
151 N_VDestroy(pdata->zlocal);
152 SUNMatDestroy(pdata->PP);
153 free(pdata); pdata = NULL;
154 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDABBDPRE",
155 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
156 return(IDALS_MEM_FAIL);
157 }
158 pdata->tempv2 = NULL;
159 pdata->tempv2 = N_VClone(IDA_mem->ida_tempv1);
160 if (pdata->tempv2 == NULL){
161 N_VDestroy(pdata->rlocal);
162 N_VDestroy(pdata->zlocal);
163 N_VDestroy(pdata->tempv1);
164 SUNMatDestroy(pdata->PP);
165 free(pdata); pdata = NULL;
166 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDABBDPRE",
167 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
168 return(IDALS_MEM_FAIL);
169 }
170 pdata->tempv3 = NULL;
171 pdata->tempv3 = N_VClone(IDA_mem->ida_tempv1);
172 if (pdata->tempv3 == NULL){
173 N_VDestroy(pdata->rlocal);
174 N_VDestroy(pdata->zlocal);
175 N_VDestroy(pdata->tempv1);
176 N_VDestroy(pdata->tempv2);
177 SUNMatDestroy(pdata->PP);
178 free(pdata); pdata = NULL;
179 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDABBDPRE",
180 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
181 return(IDALS_MEM_FAIL);
182 }
183 pdata->tempv4 = NULL;
184 pdata->tempv4 = N_VClone(IDA_mem->ida_tempv1);
185 if (pdata->tempv4 == NULL){
186 N_VDestroy(pdata->rlocal);
187 N_VDestroy(pdata->zlocal);
188 N_VDestroy(pdata->tempv1);
189 N_VDestroy(pdata->tempv2);
190 N_VDestroy(pdata->tempv3);
191 SUNMatDestroy(pdata->PP);
192 free(pdata); pdata = NULL;
193 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDABBDPRE",
194 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
195 return(IDALS_MEM_FAIL);
196 }
197
198 /* Allocate memory for banded linear solver */
199 pdata->LS = NULL;
200 pdata->LS = SUNLinSol_Band(pdata->rlocal, pdata->PP);
201 if (pdata->LS == NULL) {
202 N_VDestroy(pdata->zlocal);
203 N_VDestroy(pdata->rlocal);
204 N_VDestroy(pdata->tempv1);
205 N_VDestroy(pdata->tempv2);
206 N_VDestroy(pdata->tempv3);
207 N_VDestroy(pdata->tempv4);
208 SUNMatDestroy(pdata->PP);
209 free(pdata); pdata = NULL;
210 IDAProcessError(IDA_mem, IDALS_MEM_FAIL, "IDABBDPRE",
211 "IDABBDPrecInit", MSGBBD_MEM_FAIL);
212 return(IDALS_MEM_FAIL);
213 }
214
215 /* initialize band linear solver object */
216 flag = SUNLinSolInitialize(pdata->LS);
217 if (flag != SUNLS_SUCCESS) {
218 N_VDestroy(pdata->zlocal);
219 N_VDestroy(pdata->rlocal);
220 N_VDestroy(pdata->tempv1);
221 N_VDestroy(pdata->tempv2);
222 N_VDestroy(pdata->tempv3);
223 N_VDestroy(pdata->tempv4);
224 SUNMatDestroy(pdata->PP);
225 SUNLinSolFree(pdata->LS);
226 free(pdata); pdata = NULL;
227 IDAProcessError(IDA_mem, IDALS_SUNLS_FAIL, "IDABBDPRE",
228 "IDABBDPrecInit", MSGBBD_SUNLS_FAIL);
229 return(IDALS_SUNLS_FAIL);
230 }
231
232 /* Set rel_yy based on input value dq_rel_yy (0 implies default). */
233 pdata->rel_yy = (dq_rel_yy > ZERO) ?
234 dq_rel_yy : SUNRsqrt(IDA_mem->ida_uround);
235
236 /* Store Nlocal to be used in IDABBDPrecSetup */
237 pdata->n_local = Nlocal;
238
239 /* Set work space sizes and initialize nge. */
240 pdata->rpwsize = 0;
241 pdata->ipwsize = 0;
242 if (IDA_mem->ida_tempv1->ops->nvspace) {
243 N_VSpace(IDA_mem->ida_tempv1, &lrw1, &liw1);
244 pdata->rpwsize += 4*lrw1;
245 pdata->ipwsize += 4*liw1;
246 }
247 if (pdata->rlocal->ops->nvspace) {
248 N_VSpace(pdata->rlocal, &lrw1, &liw1);
249 pdata->rpwsize += 2*lrw1;
250 pdata->ipwsize += 2*liw1;
251 }
252 if (pdata->PP->ops->space) {
253 flag = SUNMatSpace(pdata->PP, &lrw, &liw);
254 pdata->rpwsize += lrw;
255 pdata->ipwsize += liw;
256 }
257 if (pdata->LS->ops->space) {
258 flag = SUNLinSolSpace(pdata->LS, &lrw, &liw);
259 pdata->rpwsize += lrw;
260 pdata->ipwsize += liw;
261 }
262 pdata->nge = 0;
263
264 /* make sure pdata is free from any previous allocations */
265 if (idals_mem->pfree)
266 idals_mem->pfree(IDA_mem);
267
268 /* Point to the new pdata field in the LS memory */
269 idals_mem->pdata = pdata;
270
271 /* Attach the pfree function */
272 idals_mem->pfree = IDABBDPrecFree;
273
274 /* Attach preconditioner solve and setup functions */
275 flag = IDASetPreconditioner(ida_mem,
276 IDABBDPrecSetup,
277 IDABBDPrecSolve);
278
279 return(flag);
280 }
281
282
283 /*-------------------------------------------------------------*/
IDABBDPrecReInit(void * ida_mem,sunindextype mudq,sunindextype mldq,realtype dq_rel_yy)284 int IDABBDPrecReInit(void *ida_mem, sunindextype mudq,
285 sunindextype mldq, realtype dq_rel_yy)
286 {
287 IDAMem IDA_mem;
288 IDALsMem idals_mem;
289 IBBDPrecData pdata;
290 sunindextype Nlocal;
291
292 if (ida_mem == NULL) {
293 IDAProcessError(NULL, IDALS_MEM_NULL, "IDABBDPRE",
294 "IDABBDPrecReInit", MSGBBD_MEM_NULL);
295 return(IDALS_MEM_NULL);
296 }
297 IDA_mem = (IDAMem) ida_mem;
298
299 /* Test if the LS linear solver interface has been created */
300 if (IDA_mem->ida_lmem == NULL) {
301 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDABBDPRE",
302 "IDABBDPrecReInit", MSGBBD_LMEM_NULL);
303 return(IDALS_LMEM_NULL);
304 }
305 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
306
307 /* Test if the preconditioner data is non-NULL */
308 if (idals_mem->pdata == NULL) {
309 IDAProcessError(IDA_mem, IDALS_PMEM_NULL, "IDABBDPRE",
310 "IDABBDPrecReInit", MSGBBD_PMEM_NULL);
311 return(IDALS_PMEM_NULL);
312 }
313 pdata = (IBBDPrecData) idals_mem->pdata;
314
315 /* Load half-bandwidths. */
316 Nlocal = pdata->n_local;
317 pdata->mudq = SUNMIN(Nlocal-1, SUNMAX(0, mudq));
318 pdata->mldq = SUNMIN(Nlocal-1, SUNMAX(0, mldq));
319
320 /* Set rel_yy based on input value dq_rel_yy (0 implies default). */
321 pdata->rel_yy = (dq_rel_yy > ZERO) ?
322 dq_rel_yy : SUNRsqrt(IDA_mem->ida_uround);
323
324 /* Re-initialize nge */
325 pdata->nge = 0;
326
327 return(IDALS_SUCCESS);
328 }
329
330
331 /*-------------------------------------------------------------*/
IDABBDPrecGetWorkSpace(void * ida_mem,long int * lenrwBBDP,long int * leniwBBDP)332 int IDABBDPrecGetWorkSpace(void *ida_mem,
333 long int *lenrwBBDP,
334 long int *leniwBBDP)
335 {
336 IDAMem IDA_mem;
337 IDALsMem idals_mem;
338 IBBDPrecData pdata;
339
340 if (ida_mem == NULL) {
341 IDAProcessError(NULL, IDALS_MEM_NULL, "IDABBDPRE",
342 "IDABBDPrecGetWorkSpace", MSGBBD_MEM_NULL);
343 return(IDALS_MEM_NULL);
344 }
345 IDA_mem = (IDAMem) ida_mem;
346
347 if (IDA_mem->ida_lmem == NULL) {
348 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDABBDPRE",
349 "IDABBDPrecGetWorkSpace", MSGBBD_LMEM_NULL);
350 return(IDALS_LMEM_NULL);
351 }
352 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
353
354 if (idals_mem->pdata == NULL) {
355 IDAProcessError(IDA_mem, IDALS_PMEM_NULL, "IDABBDPRE",
356 "IDABBDPrecGetWorkSpace", MSGBBD_PMEM_NULL);
357 return(IDALS_PMEM_NULL);
358 }
359 pdata = (IBBDPrecData) idals_mem->pdata;
360
361 *lenrwBBDP = pdata->rpwsize;
362 *leniwBBDP = pdata->ipwsize;
363
364 return(IDALS_SUCCESS);
365 }
366
367
368 /*-------------------------------------------------------------*/
IDABBDPrecGetNumGfnEvals(void * ida_mem,long int * ngevalsBBDP)369 int IDABBDPrecGetNumGfnEvals(void *ida_mem,
370 long int *ngevalsBBDP)
371 {
372 IDAMem IDA_mem;
373 IDALsMem idals_mem;
374 IBBDPrecData pdata;
375
376 if (ida_mem == NULL) {
377 IDAProcessError(NULL, IDALS_MEM_NULL, "IDABBDPRE",
378 "IDABBDPrecGetNumGfnEvals", MSGBBD_MEM_NULL);
379 return(IDALS_MEM_NULL);
380 }
381 IDA_mem = (IDAMem) ida_mem;
382
383 if (IDA_mem->ida_lmem == NULL) {
384 IDAProcessError(IDA_mem, IDALS_LMEM_NULL, "IDABBDPRE",
385 "IDABBDPrecGetNumGfnEvals", MSGBBD_LMEM_NULL);
386 return(IDALS_LMEM_NULL);
387 }
388 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
389
390 if (idals_mem->pdata == NULL) {
391 IDAProcessError(IDA_mem, IDALS_PMEM_NULL, "IDABBDPRE",
392 "IDABBDPrecGetNumGfnEvals", MSGBBD_PMEM_NULL);
393 return(IDALS_PMEM_NULL);
394 }
395 pdata = (IBBDPrecData) idals_mem->pdata;
396
397 *ngevalsBBDP = pdata->nge;
398
399 return(IDALS_SUCCESS);
400 }
401
402
403 /*---------------------------------------------------------------
404 IDABBDPrecSetup:
405
406 IDABBDPrecSetup generates a band-block-diagonal preconditioner
407 matrix, where the local block (on this processor) is a band
408 matrix. Each local block is computed by a difference quotient
409 scheme via calls to the user-supplied routines glocal, gcomm.
410 After generating the block in the band matrix PP, this routine
411 does an LU factorization in place in PP.
412
413 The IDABBDPrecSetup parameters used here are as follows:
414
415 tt is the current value of the independent variable t.
416
417 yy is the current value of the dependent variable vector,
418 namely the predicted value of y(t).
419
420 yp is the current value of the derivative vector y',
421 namely the predicted value of y'(t).
422
423 c_j is the scalar in the system Jacobian, proportional to 1/hh.
424
425 bbd_data is the pointer to BBD memory set by IDABBDInit
426
427 The argument rr is not used.
428
429 Return value:
430 The value returned by this IDABBDPrecSetup function is a int
431 flag indicating whether it was successful. This value is
432 0 if successful,
433 > 0 for a recoverable error (step will be retried), or
434 < 0 for a nonrecoverable error (step fails).
435 ----------------------------------------------------------------*/
IDABBDPrecSetup(realtype tt,N_Vector yy,N_Vector yp,N_Vector rr,realtype c_j,void * bbd_data)436 static int IDABBDPrecSetup(realtype tt, N_Vector yy, N_Vector yp,
437 N_Vector rr, realtype c_j, void *bbd_data)
438 {
439 IBBDPrecData pdata;
440 IDAMem IDA_mem;
441 int retval;
442
443 pdata =(IBBDPrecData) bbd_data;
444
445 IDA_mem = (IDAMem) pdata->ida_mem;
446
447 /* Call IBBDDQJac for a new Jacobian calculation and store in PP. */
448 retval = SUNMatZero(pdata->PP);
449 retval = IBBDDQJac(pdata, tt, c_j, yy, yp, pdata->tempv1,
450 pdata->tempv2, pdata->tempv3, pdata->tempv4);
451 if (retval < 0) {
452 IDAProcessError(IDA_mem, -1, "IDABBDPRE", "IDABBDPrecSetup",
453 MSGBBD_FUNC_FAILED);
454 return(-1);
455 }
456 if (retval > 0) {
457 return(1);
458 }
459
460 /* Do LU factorization of matrix and return error flag */
461 retval = SUNLinSolSetup_Band(pdata->LS, pdata->PP);
462 return(retval);
463 }
464
465
466 /*---------------------------------------------------------------
467 IDABBDPrecSolve
468
469 The function IDABBDPrecSolve computes a solution to the linear
470 system P z = r, where P is the left preconditioner defined by
471 the routine IDABBDPrecSetup.
472
473 The IDABBDPrecSolve parameters used here are as follows:
474
475 rvec is the input right-hand side vector r.
476
477 zvec is the computed solution vector z.
478
479 bbd_data is the pointer to BBD data set by IDABBDInit.
480
481 The arguments tt, yy, yp, rr, c_j and delta are NOT used.
482
483 IDABBDPrecSolve returns the value returned from the linear
484 solver object.
485 ---------------------------------------------------------------*/
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)486 static int IDABBDPrecSolve(realtype tt, N_Vector yy, N_Vector yp,
487 N_Vector rr, N_Vector rvec, N_Vector zvec,
488 realtype c_j, realtype delta, void *bbd_data)
489 {
490 IBBDPrecData pdata;
491 int retval;
492
493 pdata = (IBBDPrecData) bbd_data;
494
495 /* Attach local data arrays for rvec and zvec to rlocal and zlocal */
496 N_VSetArrayPointer(N_VGetArrayPointer(rvec), pdata->rlocal);
497 N_VSetArrayPointer(N_VGetArrayPointer(zvec), pdata->zlocal);
498
499 /* Call banded solver object to do the work */
500 retval = SUNLinSolSolve(pdata->LS, pdata->PP, pdata->zlocal,
501 pdata->rlocal, ZERO);
502
503 /* Detach local data arrays from rlocal and zlocal */
504 N_VSetArrayPointer(NULL, pdata->rlocal);
505 N_VSetArrayPointer(NULL, pdata->zlocal);
506
507 return(retval);
508 }
509
510
511 /*-------------------------------------------------------------*/
IDABBDPrecFree(IDAMem IDA_mem)512 static int IDABBDPrecFree(IDAMem IDA_mem)
513 {
514 IDALsMem idals_mem;
515 IBBDPrecData pdata;
516
517 if (IDA_mem->ida_lmem == NULL) return(0);
518 idals_mem = (IDALsMem) IDA_mem->ida_lmem;
519
520 if (idals_mem->pdata == NULL) return(0);
521 pdata = (IBBDPrecData) idals_mem->pdata;
522
523 SUNLinSolFree(pdata->LS);
524 N_VDestroy(pdata->rlocal);
525 N_VDestroy(pdata->zlocal);
526 N_VDestroy(pdata->tempv1);
527 N_VDestroy(pdata->tempv2);
528 N_VDestroy(pdata->tempv3);
529 N_VDestroy(pdata->tempv4);
530 SUNMatDestroy(pdata->PP);
531
532 free(pdata);
533 pdata = NULL;
534
535 return(0);
536 }
537
538
539 /*---------------------------------------------------------------
540 IBBDDQJac
541
542 This routine generates a banded difference quotient approximation
543 to the local block of the Jacobian of G(t,y,y'). It assumes that
544 a band matrix of type SUNMatrix is stored column-wise, and that
545 elements within each column are contiguous.
546
547 All matrix elements are generated as difference quotients, by way
548 of calls to the user routine glocal. By virtue of the band
549 structure, the number of these calls is bandwidth + 1, where
550 bandwidth = mldq + mudq + 1. But the band matrix kept has
551 bandwidth = mlkeep + mukeep + 1. This routine also assumes that
552 the local elements of a vector are stored contiguously.
553
554 Return values are: 0 (success), > 0 (recoverable error),
555 or < 0 (nonrecoverable error).
556 ----------------------------------------------------------------*/
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)557 static int IBBDDQJac(IBBDPrecData pdata, realtype tt, realtype cj,
558 N_Vector yy, N_Vector yp, N_Vector gref,
559 N_Vector ytemp, N_Vector yptemp, N_Vector gtemp)
560 {
561 IDAMem IDA_mem;
562 realtype inc, inc_inv;
563 int retval;
564 sunindextype group, i, j, width, ngroups, i1, i2;
565 realtype *ydata, *ypdata, *ytempdata, *yptempdata, *grefdata, *gtempdata;
566 realtype *cnsdata = NULL, *ewtdata;
567 realtype *col_j, conj, yj, ypj, ewtj;
568
569 IDA_mem = (IDAMem) pdata->ida_mem;
570
571 /* Initialize ytemp and yptemp. */
572 N_VScale(ONE, yy, ytemp);
573 N_VScale(ONE, yp, yptemp);
574
575 /* Obtain pointers as required to the data array of vectors. */
576 ydata = N_VGetArrayPointer(yy);
577 ypdata = N_VGetArrayPointer(yp);
578 gtempdata = N_VGetArrayPointer(gtemp);
579 ewtdata = N_VGetArrayPointer(IDA_mem->ida_ewt);
580 if (IDA_mem->ida_constraintsSet)
581 cnsdata = N_VGetArrayPointer(IDA_mem->ida_constraints);
582 ytempdata = N_VGetArrayPointer(ytemp);
583 yptempdata= N_VGetArrayPointer(yptemp);
584 grefdata = N_VGetArrayPointer(gref);
585
586 /* Call gcomm and glocal to get base value of G(t,y,y'). */
587 if (pdata->gcomm != NULL) {
588 retval = pdata->gcomm(pdata->n_local, tt, yy, yp, IDA_mem->ida_user_data);
589 if (retval != 0) return(retval);
590 }
591
592 retval = pdata->glocal(pdata->n_local, tt, yy, yp, gref, IDA_mem->ida_user_data);
593 pdata->nge++;
594 if (retval != 0) return(retval);
595
596 /* Set bandwidth and number of column groups for band differencing. */
597 width = pdata->mldq + pdata->mudq + 1;
598 ngroups = SUNMIN(width, pdata->n_local);
599
600 /* Loop over groups. */
601 for(group = 1; group <= ngroups; group++) {
602
603 /* Loop over the components in this group. */
604 for(j = group-1; j < pdata->n_local; j += width) {
605 yj = ydata[j];
606 ypj = ypdata[j];
607 ewtj = ewtdata[j];
608
609 /* Set increment inc to yj based on rel_yy*abs(yj), with
610 adjustments using ypj and ewtj if this is small, and a further
611 adjustment to give it the same sign as hh*ypj. */
612 inc = pdata->rel_yy *
613 SUNMAX(SUNRabs(yj), SUNMAX( SUNRabs(IDA_mem->ida_hh*ypj), ONE/ewtj));
614 if (IDA_mem->ida_hh*ypj < ZERO) inc = -inc;
615 inc = (yj + inc) - yj;
616
617 /* Adjust sign(inc) again if yj has an inequality constraint. */
618 if (IDA_mem->ida_constraintsSet) {
619 conj = cnsdata[j];
620 if (SUNRabs(conj) == ONE) {if ((yj+inc)*conj < ZERO) inc = -inc;}
621 else if (SUNRabs(conj) == TWO) {if ((yj+inc)*conj <= ZERO) inc = -inc;}
622 }
623
624 /* Increment yj and ypj. */
625 ytempdata[j] += inc;
626 yptempdata[j] += cj*inc;
627
628 }
629
630 /* Evaluate G with incremented y and yp arguments. */
631 retval = pdata->glocal(pdata->n_local, tt, ytemp, yptemp,
632 gtemp, IDA_mem->ida_user_data);
633 pdata->nge++;
634 if (retval != 0) return(retval);
635
636 /* Loop over components of the group again; restore ytemp and yptemp. */
637 for(j = group-1; j < pdata->n_local; j += width) {
638 yj = ytempdata[j] = ydata[j];
639 ypj = yptempdata[j] = ypdata[j];
640 ewtj = ewtdata[j];
641
642 /* Set increment inc as before .*/
643 inc = pdata->rel_yy *
644 SUNMAX(SUNRabs(yj), SUNMAX( SUNRabs(IDA_mem->ida_hh*ypj), ONE/ewtj));
645 if (IDA_mem->ida_hh*ypj < ZERO) inc = -inc;
646 inc = (yj + inc) - yj;
647 if (IDA_mem->ida_constraintsSet) {
648 conj = cnsdata[j];
649 if (SUNRabs(conj) == ONE) {if ((yj+inc)*conj < ZERO) inc = -inc;}
650 else if (SUNRabs(conj) == TWO) {if ((yj+inc)*conj <= ZERO) inc = -inc;}
651 }
652
653 /* Form difference quotients and load into PP. */
654 inc_inv = ONE/inc;
655 col_j = SUNBandMatrix_Column(pdata->PP,j);
656 i1 = SUNMAX(0, j - pdata->mukeep);
657 i2 = SUNMIN(j + pdata->mlkeep, pdata->n_local-1);
658 for(i=i1; i <= i2; i++)
659 SM_COLUMN_ELEMENT_B(col_j,i,j) =
660 inc_inv * (gtempdata[i] - grefdata[i]);
661 }
662 }
663
664 return(0);
665 }
666