1 /* -----------------------------------------------------------------
2 * Programmer(s): David J. Gardner @ LLNL
3 * Allan Taylor, Alan Hindmarsh, Radu Serban, and
4 * Aaron Collier @ LLNL
5 * -----------------------------------------------------------------
6 * SUNDIALS Copyright Start
7 * Copyright (c) 2002-2020, 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 KINSol and the
19 * KINLS linear solver interface.
20 *
21 * Note: With only one process, a banded matrix results
22 * rather than a b-b-d matrix with banded blocks. Diagonal
23 * blocking occurs at the process level.
24 * -----------------------------------------------------------------*/
25
26 #include <stdio.h>
27 #include <stdlib.h>
28
29 #include "kinsol_impl.h"
30 #include "kinsol_ls_impl.h"
31 #include "kinsol_bbdpre_impl.h"
32
33 #include <sundials/sundials_math.h>
34 #include <nvector/nvector_serial.h>
35
36 #define ZERO RCONST(0.0)
37 #define ONE RCONST(1.0)
38
39 /* Prototypes of functions KINBBDPrecSetup and KINBBDPrecSolve */
40 static int KINBBDPrecSetup(N_Vector uu, N_Vector uscale,
41 N_Vector fval, N_Vector fscale,
42 void *pdata);
43
44 static int KINBBDPrecSolve(N_Vector uu, N_Vector uscale,
45 N_Vector fval, N_Vector fscale,
46 N_Vector vv, void *pdata);
47
48 /* Prototype for KINBBDPrecFree */
49 static int KINBBDPrecFree(KINMem kin_mem);
50
51 /* Prototype for difference quotient jacobian calculation routine */
52 static int KBBDDQJac(KBBDPrecData pdata,
53 N_Vector uu, N_Vector uscale,
54 N_Vector gu, N_Vector gtemp, N_Vector utemp);
55
56 /*------------------------------------------------------------------
57 user-callable functions
58 ------------------------------------------------------------------*/
59
60 /*------------------------------------------------------------------
61 KINBBDPrecInit
62 ------------------------------------------------------------------*/
KINBBDPrecInit(void * kinmem,sunindextype Nlocal,sunindextype mudq,sunindextype mldq,sunindextype mukeep,sunindextype mlkeep,realtype dq_rel_uu,KINBBDLocalFn gloc,KINBBDCommFn gcomm)63 int KINBBDPrecInit(void *kinmem, sunindextype Nlocal,
64 sunindextype mudq, sunindextype mldq,
65 sunindextype mukeep, sunindextype mlkeep,
66 realtype dq_rel_uu,
67 KINBBDLocalFn gloc, KINBBDCommFn gcomm)
68 {
69 KINMem kin_mem;
70 KINLsMem kinls_mem;
71 KBBDPrecData pdata;
72 sunindextype muk, mlk, storage_mu, lrw1, liw1;
73 long int lrw, liw;
74 int flag;
75
76 if (kinmem == NULL) {
77 KINProcessError(NULL, KINLS_MEM_NULL, "KINBBDPRE",
78 "KINBBDPrecInit", MSGBBD_MEM_NULL);
79 return(KINLS_MEM_NULL);
80 }
81 kin_mem = (KINMem) kinmem;
82
83 /* Test if the LS linear solver interface has been created */
84 if (kin_mem->kin_lmem == NULL) {
85 KINProcessError(kin_mem, KINLS_LMEM_NULL, "KINBBDPRE",
86 "KINBBDPrecInit", MSGBBD_LMEM_NULL);
87 return(KINLS_LMEM_NULL);
88 }
89 kinls_mem = (KINLsMem) kin_mem->kin_lmem;
90
91 /* Test compatibility of NVECTOR package with the BBD preconditioner */
92 /* Note: Do NOT need to check for N_VScale since has already been checked for in KINSOL */
93 if (kin_mem->kin_vtemp1->ops->nvgetarraypointer == NULL) {
94 KINProcessError(kin_mem, KINLS_ILL_INPUT, "KINBBDPRE",
95 "KINBBDPrecInit", MSGBBD_BAD_NVECTOR);
96 return(KINLS_ILL_INPUT);
97 }
98
99 /* Allocate data memory */
100 pdata = NULL;
101 pdata = (KBBDPrecData) malloc(sizeof *pdata);
102 if (pdata == NULL) {
103 KINProcessError(kin_mem, KINLS_MEM_FAIL,
104 "KINBBDPRE", "KINBBDPrecInit", MSGBBD_MEM_FAIL);
105 return(KINLS_MEM_FAIL);
106 }
107
108 /* Set pointers to gloc and gcomm; load half-bandwidths */
109 pdata->kin_mem = kinmem;
110 pdata->gloc = gloc;
111 pdata->gcomm = gcomm;
112 pdata->mudq = SUNMIN(Nlocal-1, SUNMAX(0, mudq));
113 pdata->mldq = SUNMIN(Nlocal-1, SUNMAX(0, mldq));
114 muk = SUNMIN(Nlocal-1, SUNMAX(0, mukeep));
115 mlk = SUNMIN(Nlocal-1, SUNMAX(0, mlkeep));
116 pdata->mukeep = muk;
117 pdata->mlkeep = mlk;
118
119 /* Set extended upper half-bandwidth for PP (required for pivoting) */
120 storage_mu = SUNMIN(Nlocal-1, muk+mlk);
121
122 /* Allocate memory for preconditioner matrix */
123 pdata->PP = NULL;
124 pdata->PP = SUNBandMatrixStorage(Nlocal, muk, mlk, storage_mu);
125 if (pdata->PP == NULL) {
126 free(pdata); pdata = NULL;
127 KINProcessError(kin_mem, KINLS_MEM_FAIL, "KINBBDPRE",
128 "KINBBDPrecInit", MSGBBD_MEM_FAIL);
129 return(KINLS_MEM_FAIL);
130 }
131
132 /* Allocate memory for temporary N_Vectors */
133 pdata->zlocal = NULL;
134 pdata->zlocal = N_VNew_Serial(Nlocal);
135 if (pdata->zlocal == NULL) {
136 SUNMatDestroy(pdata->PP);
137 free(pdata); pdata = NULL;
138 KINProcessError(kin_mem, KINLS_MEM_FAIL, "KINBBDPRE",
139 "KINBBDPrecInit", MSGBBD_MEM_FAIL);
140 return(KINLS_MEM_FAIL);
141 }
142
143 pdata->rlocal = NULL;
144 pdata->rlocal = N_VNewEmpty_Serial(Nlocal); /* empty vector */
145 if (pdata->rlocal == NULL) {
146 N_VDestroy(pdata->zlocal);
147 SUNMatDestroy(pdata->PP);
148 free(pdata); pdata = NULL;
149 KINProcessError(kin_mem, KINLS_MEM_FAIL, "KINBBDPRE",
150 "KINBBDPrecInit", MSGBBD_MEM_FAIL);
151 return(KINLS_MEM_FAIL);
152 }
153
154 pdata->tempv1 = NULL;
155 pdata->tempv1 = N_VClone(kin_mem->kin_vtemp1);
156 if (pdata->tempv1 == NULL) {
157 N_VDestroy(pdata->zlocal);
158 N_VDestroy(pdata->rlocal);
159 SUNMatDestroy(pdata->PP);
160 free(pdata); pdata = NULL;
161 KINProcessError(kin_mem, KINLS_MEM_FAIL, "KINBBDPRE",
162 "KINBBDPrecInit", MSGBBD_MEM_FAIL);
163 return(KINLS_MEM_FAIL);
164 }
165
166 pdata->tempv2 = NULL;
167 pdata->tempv2 = N_VClone(kin_mem->kin_vtemp1);
168 if (pdata->tempv2 == NULL) {
169 N_VDestroy(pdata->zlocal);
170 N_VDestroy(pdata->rlocal);
171 N_VDestroy(pdata->tempv1);
172 SUNMatDestroy(pdata->PP);
173 free(pdata); pdata = NULL;
174 KINProcessError(kin_mem, KINLS_MEM_FAIL, "KINBBDPRE",
175 "KINBBDPrecInit", MSGBBD_MEM_FAIL);
176 return(KINLS_MEM_FAIL);
177 }
178
179 pdata->tempv3 = NULL;
180 pdata->tempv3 = N_VClone(kin_mem->kin_vtemp1);
181 if (pdata->tempv3 == NULL) {
182 N_VDestroy(pdata->zlocal);
183 N_VDestroy(pdata->rlocal);
184 N_VDestroy(pdata->tempv1);
185 N_VDestroy(pdata->tempv2);
186 SUNMatDestroy(pdata->PP);
187 free(pdata); pdata = NULL;
188 KINProcessError(kin_mem, KINLS_MEM_FAIL, "KINBBDPRE",
189 "KINBBDPrecInit", MSGBBD_MEM_FAIL);
190 return(KINLS_MEM_FAIL);
191 }
192
193 /* Allocate memory for banded linear solver */
194 pdata->LS = NULL;
195 pdata->LS = SUNLinSol_Band(pdata->zlocal, pdata->PP);
196 if (pdata->LS == NULL) {
197 N_VDestroy(pdata->zlocal);
198 N_VDestroy(pdata->rlocal);
199 N_VDestroy(pdata->tempv1);
200 N_VDestroy(pdata->tempv2);
201 N_VDestroy(pdata->tempv3);
202 SUNMatDestroy(pdata->PP);
203 free(pdata); pdata = NULL;
204 KINProcessError(kin_mem, KINLS_MEM_FAIL, "KINBBDPRE",
205 "KINBBDPrecInit", MSGBBD_MEM_FAIL);
206 return(KINLS_MEM_FAIL);
207 }
208
209 /* initialize band linear solver object */
210 flag = SUNLinSolInitialize(pdata->LS);
211 if (flag != SUNLS_SUCCESS) {
212 N_VDestroy(pdata->zlocal);
213 N_VDestroy(pdata->rlocal);
214 N_VDestroy(pdata->tempv1);
215 N_VDestroy(pdata->tempv2);
216 N_VDestroy(pdata->tempv3);
217 SUNMatDestroy(pdata->PP);
218 SUNLinSolFree(pdata->LS);
219 free(pdata); pdata = NULL;
220 KINProcessError(kin_mem, KINLS_SUNLS_FAIL, "KINBBDPRE",
221 "KINBBDPrecInit", MSGBBD_SUNLS_FAIL);
222 return(KINLS_SUNLS_FAIL);
223 }
224
225 /* Set rel_uu based on input value dq_rel_uu (0 implies default) */
226 pdata->rel_uu = (dq_rel_uu > ZERO) ? dq_rel_uu : SUNRsqrt(kin_mem->kin_uround);
227
228 /* Store Nlocal to be used in KINBBDPrecSetup */
229 pdata->n_local = Nlocal;
230
231 /* Set work space sizes and initialize nge */
232 pdata->rpwsize = 0;
233 pdata->ipwsize = 0;
234 if (kin_mem->kin_vtemp1->ops->nvspace) {
235 N_VSpace(kin_mem->kin_vtemp1, &lrw1, &liw1);
236 pdata->rpwsize += 3*lrw1;
237 pdata->ipwsize += 3*liw1;
238 }
239 if (pdata->zlocal->ops->nvspace) {
240 N_VSpace(pdata->zlocal, &lrw1, &liw1);
241 pdata->rpwsize += lrw1;
242 pdata->ipwsize += liw1;
243 }
244 if (pdata->rlocal->ops->nvspace) {
245 N_VSpace(pdata->rlocal, &lrw1, &liw1);
246 pdata->rpwsize += lrw1;
247 pdata->ipwsize += liw1;
248 }
249 if (pdata->PP->ops->space) {
250 flag = SUNMatSpace(pdata->PP, &lrw, &liw);
251 pdata->rpwsize += lrw;
252 pdata->ipwsize += liw;
253 }
254 if (pdata->LS->ops->space) {
255 flag = SUNLinSolSpace(pdata->LS, &lrw, &liw);
256 pdata->rpwsize += lrw;
257 pdata->ipwsize += liw;
258 }
259 pdata->nge = 0;
260
261 /* make sure pdata is free from any previous allocations */
262 if (kinls_mem->pfree != NULL)
263 kinls_mem->pfree(kin_mem);
264
265 /* Point to the new pdata field in the LS memory */
266 kinls_mem->pdata = pdata;
267
268 /* Attach the pfree function */
269 kinls_mem->pfree = KINBBDPrecFree;
270
271 /* Attach preconditioner solve and setup functions */
272 flag = KINSetPreconditioner(kinmem, KINBBDPrecSetup,
273 KINBBDPrecSolve);
274
275 return(flag);
276 }
277
278
279 /*------------------------------------------------------------------
280 KINBBDPrecGetWorkSpace
281 ------------------------------------------------------------------*/
KINBBDPrecGetWorkSpace(void * kinmem,long int * lenrwBBDP,long int * leniwBBDP)282 int KINBBDPrecGetWorkSpace(void *kinmem,
283 long int *lenrwBBDP,
284 long int *leniwBBDP)
285 {
286 KINMem kin_mem;
287 KINLsMem kinls_mem;
288 KBBDPrecData pdata;
289
290 if (kinmem == NULL) {
291 KINProcessError(NULL, KINLS_MEM_NULL, "KINBBDPRE",
292 "KINBBDPrecGetWorkSpace", MSGBBD_MEM_NULL);
293 return(KINLS_MEM_NULL);
294 }
295 kin_mem = (KINMem) kinmem;
296
297 if (kin_mem->kin_lmem == NULL) {
298 KINProcessError(kin_mem, KINLS_LMEM_NULL, "KINBBDPRE",
299 "KINBBDPrecGetWorkSpace", MSGBBD_LMEM_NULL);
300 return(KINLS_LMEM_NULL);
301 }
302 kinls_mem = (KINLsMem) kin_mem->kin_lmem;
303
304 if (kinls_mem->pdata == NULL) {
305 KINProcessError(kin_mem, KINLS_PMEM_NULL, "KINBBDPRE",
306 "KINBBDPrecGetWorkSpace", MSGBBD_PMEM_NULL);
307 return(KINLS_PMEM_NULL);
308 }
309 pdata = (KBBDPrecData) kinls_mem->pdata;
310
311 *lenrwBBDP = pdata->rpwsize;
312 *leniwBBDP = pdata->ipwsize;
313
314 return(KINLS_SUCCESS);
315 }
316
317 /*------------------------------------------------------------------
318 KINBBDPrecGetNumGfnEvals
319 -------------------------------------------------------------------*/
KINBBDPrecGetNumGfnEvals(void * kinmem,long int * ngevalsBBDP)320 int KINBBDPrecGetNumGfnEvals(void *kinmem,
321 long int *ngevalsBBDP)
322 {
323 KINMem kin_mem;
324 KINLsMem kinls_mem;
325 KBBDPrecData pdata;
326
327 if (kinmem == NULL) {
328 KINProcessError(NULL, KINLS_MEM_NULL, "KINBBDPRE",
329 "KINBBDPrecGetNumGfnEvals", MSGBBD_MEM_NULL);
330 return(KINLS_MEM_NULL);
331 }
332 kin_mem = (KINMem) kinmem;
333
334 if (kin_mem->kin_lmem == NULL) {
335 KINProcessError(kin_mem, KINLS_LMEM_NULL, "KINBBDPRE",
336 "KINBBDPrecGetNumGfnEvals", MSGBBD_LMEM_NULL);
337 return(KINLS_LMEM_NULL);
338 }
339 kinls_mem = (KINLsMem) kin_mem->kin_lmem;
340
341 if (kinls_mem->pdata == NULL) {
342 KINProcessError(kin_mem, KINLS_PMEM_NULL, "KINBBDPRE",
343 "KINBBDPrecGetNumGfnEvals", MSGBBD_PMEM_NULL);
344 return(KINLS_PMEM_NULL);
345 }
346 pdata = (KBBDPrecData) kinls_mem->pdata;
347
348 *ngevalsBBDP = pdata->nge;
349
350 return(KINLS_SUCCESS);
351 }
352
353
354 /*------------------------------------------------------------------
355 KINBBDPrecSetup
356
357 KINBBDPrecSetup generates and factors a banded block of the
358 preconditioner matrix on each processor, via calls to the
359 user-supplied gloc and gcomm functions. It uses difference
360 quotient approximations to the Jacobian elements.
361
362 KINBBDPrecSetup calculates a new Jacobian, stored in banded
363 matrix PP and does an LU factorization of P in place in PP.
364
365 The parameters of KINBBDPrecSetup are as follows:
366
367 uu is the current value of the dependent variable vector,
368 namely the solutin to func(uu)=0
369
370 uscale is the dependent variable scaling vector (i.e. uu)
371
372 fval is the vector f(u)
373
374 fscale is the function scaling vector
375
376 bbd_data is the pointer to BBD data set by KINBBDInit.
377
378 Note: The value to be returned by the KINBBDPrecSetup function
379 is a flag indicating whether it was successful. This value is:
380 0 if successful,
381 > 0 for a recoverable error - step will be retried.
382 ------------------------------------------------------------------*/
KINBBDPrecSetup(N_Vector uu,N_Vector uscale,N_Vector fval,N_Vector fscale,void * bbd_data)383 static int KINBBDPrecSetup(N_Vector uu, N_Vector uscale,
384 N_Vector fval, N_Vector fscale,
385 void *bbd_data)
386 {
387 KBBDPrecData pdata;
388 KINMem kin_mem;
389 int retval;
390
391 pdata = (KBBDPrecData) bbd_data;
392
393 kin_mem = (KINMem) pdata->kin_mem;
394
395 /* Call KBBDDQJac for a new Jacobian calculation and store in PP */
396 retval = SUNMatZero(pdata->PP);
397 if (retval != 0) {
398 KINProcessError(kin_mem, -1, "KINBBDPRE", "KINBBDPrecSetup",
399 MSGBBD_SUNMAT_FAIL);
400 return(-1);
401 }
402
403 retval = KBBDDQJac(pdata, uu, uscale,
404 pdata->tempv1, pdata->tempv2, pdata->tempv3);
405 if (retval != 0) {
406 KINProcessError(kin_mem, -1, "KINBBDPRE", "KINBBDPrecSetup",
407 MSGBBD_FUNC_FAILED);
408 return(-1);
409 }
410
411 /* Do LU factorization of P and return error flag */
412 retval = SUNLinSolSetup_Band(pdata->LS, pdata->PP);
413 return(retval);
414 }
415
416 /*------------------------------------------------------------------
417 INBBDPrecSolve
418
419 KINBBDPrecSolve solves a linear system P z = r, with the
420 banded blocked preconditioner matrix P generated and factored
421 by KINBBDPrecSetup. Here, r comes in as vv and z is
422 returned in vv as well.
423
424 The parameters for KINBBDPrecSolve are as follows:
425
426 uu an N_Vector giving the current iterate for the system
427
428 uscale an N_Vector giving the diagonal entries of the
429 uu scaling matrix
430
431 fval an N_Vector giving the current function value
432
433 fscale an N_Vector giving the diagonal entries of the
434 function scaling matrix
435
436 vv vector initially set to the right-hand side vector r, but
437 which upon return contains a solution of the linear system
438 P*z = r
439
440 bbd_data is the pointer to BBD data set by KINBBDInit.
441
442 Note: The value returned by the KINBBDPrecSolve function is a
443 flag returned from the lienar solver object.
444 ------------------------------------------------------------------*/
445
KINBBDPrecSolve(N_Vector uu,N_Vector uscale,N_Vector fval,N_Vector fscale,N_Vector vv,void * bbd_data)446 static int KINBBDPrecSolve(N_Vector uu, N_Vector uscale,
447 N_Vector fval, N_Vector fscale,
448 N_Vector vv, void *bbd_data)
449 {
450 KBBDPrecData pdata;
451 realtype *vd;
452 realtype *zd;
453 int i, retval;
454
455 pdata = (KBBDPrecData) bbd_data;
456
457 /* Get data pointers */
458 vd = N_VGetArrayPointer(vv);
459 zd = N_VGetArrayPointer(pdata->zlocal);
460
461 /* Attach local data array for vv to rlocal */
462 N_VSetArrayPointer(vd, pdata->rlocal);
463
464 /* Call banded solver object to do the work */
465 retval = SUNLinSolSolve(pdata->LS, pdata->PP, pdata->zlocal,
466 pdata->rlocal, ZERO);
467
468 /* Copy result into vv */
469 for (i=0; i<pdata->n_local; i++)
470 vd[i] = zd[i];
471
472 return(retval);
473 }
474
475
476 /*------------------------------------------------------------------
477 KINBBDPrecFree
478 ------------------------------------------------------------------*/
KINBBDPrecFree(KINMem kin_mem)479 static int KINBBDPrecFree(KINMem kin_mem)
480 {
481 KINLsMem kinls_mem;
482 KBBDPrecData pdata;
483
484 if (kin_mem->kin_lmem == NULL) return(0);
485 kinls_mem = (KINLsMem) kin_mem->kin_lmem;
486
487 if (kinls_mem->pdata == NULL) return(0);
488 pdata = (KBBDPrecData) kinls_mem->pdata;
489
490 SUNLinSolFree(pdata->LS);
491 N_VDestroy(pdata->zlocal);
492 N_VDestroy(pdata->rlocal);
493 N_VDestroy(pdata->tempv1);
494 N_VDestroy(pdata->tempv2);
495 N_VDestroy(pdata->tempv3);
496 SUNMatDestroy(pdata->PP);
497
498 free(pdata);
499 pdata = NULL;
500
501 return(0);
502 }
503
504
505 /*------------------------------------------------------------------
506 KBBDDQJac
507
508 This routine generates a banded difference quotient
509 approximation to the Jacobian of f(u). It assumes that a band
510 matrix of type SUNMatrix is stored column-wise, and that elements
511 within each column are contiguous. All matrix elements are
512 generated as difference quotients, by way of calls to the user
513 routine gloc. By virtue of the band structure, the number of
514 these calls is bandwidth + 1, where bandwidth = ml + mu + 1.
515 This routine also assumes that the local elements of a vector
516 are stored contiguously.
517 ------------------------------------------------------------------*/
KBBDDQJac(KBBDPrecData pdata,N_Vector uu,N_Vector uscale,N_Vector gu,N_Vector gtemp,N_Vector utemp)518 static int KBBDDQJac(KBBDPrecData pdata,
519 N_Vector uu, N_Vector uscale,
520 N_Vector gu, N_Vector gtemp, N_Vector utemp)
521 {
522 KINMem kin_mem;
523 realtype inc, inc_inv;
524 int retval;
525 sunindextype group, i, j, width, ngroups, i1, i2;
526 realtype *udata, *uscdata, *gudata, *gtempdata, *utempdata, *col_j;
527
528 kin_mem = (KINMem) pdata->kin_mem;
529
530 /* load utemp with uu = predicted solution vector */
531 N_VScale(ONE, uu, utemp);
532
533 /* set pointers to the data for all vectors */
534 udata = N_VGetArrayPointer(uu);
535 uscdata = N_VGetArrayPointer(uscale);
536 gudata = N_VGetArrayPointer(gu);
537 gtempdata = N_VGetArrayPointer(gtemp);
538 utempdata = N_VGetArrayPointer(utemp);
539
540 /* Call gcomm and gloc to get base value of g(uu) */
541 if (pdata->gcomm != NULL) {
542 retval = pdata->gcomm(pdata->n_local, uu, kin_mem->kin_user_data);
543 if (retval != 0) return(retval);
544 }
545
546 retval = pdata->gloc(pdata->n_local, uu, gu, kin_mem->kin_user_data);
547 pdata->nge++;
548 if (retval != 0) return(retval);
549
550 /* Set bandwidth and number of column groups for band differencing */
551 width = pdata->mldq + pdata->mudq + 1;
552 ngroups = SUNMIN(width, pdata->n_local);
553
554 /* Loop over groups */
555 for(group = 1; group <= ngroups; group++) {
556
557 /* increment all u_j in group */
558 for(j = group - 1; j < pdata->n_local; j += width) {
559 inc = pdata->rel_uu * SUNMAX(SUNRabs(udata[j]), (ONE / uscdata[j]));
560 utempdata[j] += inc;
561 }
562
563 /* Evaluate g with incremented u */
564 retval = pdata->gloc(pdata->n_local, utemp, gtemp, kin_mem->kin_user_data);
565 pdata->nge++;
566 if (retval != 0) return(retval);
567
568 /* restore utemp, then form and load difference quotients */
569 for (j = group - 1; j < pdata->n_local; j += width) {
570 utempdata[j] = udata[j];
571 col_j = SUNBandMatrix_Column(pdata->PP,j);
572 inc = pdata->rel_uu * SUNMAX(SUNRabs(udata[j]) , (ONE / uscdata[j]));
573 inc_inv = ONE / inc;
574 i1 = SUNMAX(0, (j - pdata->mukeep));
575 i2 = SUNMIN((j + pdata->mlkeep), (pdata->n_local - 1));
576 for (i = i1; i <= i2; i++)
577 SM_COLUMN_ELEMENT_B(col_j, i, j) = inc_inv * (gtempdata[i] - gudata[i]);
578 }
579 }
580
581 return(0);
582 }
583