1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 1.6 $
4 * $Date: 2010/12/01 22:43:33 $
5 * -----------------------------------------------------------------
6 * Programmer: Radu Serban @ LLNL
7 * -----------------------------------------------------------------
8 * Copyright (c) 2006, The Regents of the University of California.
9 * Produced at the Lawrence Livermore National Laboratory.
10 * All rights reserved.
11 * For details, see the LICENSE file.
12 * -----------------------------------------------------------------
13 * This is the implementation file for the KINDLS linear solvers
14 * -----------------------------------------------------------------
15 */
16
17 /*
18 * =================================================================
19 * IMPORTED HEADER FILES
20 * =================================================================
21 */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25
26 #include "kinsol_impl.h"
27 #include "kinsol_direct_impl.h"
28 #include <sundials/sundials_math.h>
29
30 /*
31 * =================================================================
32 * FUNCTION SPECIFIC CONSTANTS
33 * =================================================================
34 */
35
36 /* Constant for DQ Jacobian approximation */
37 #define MIN_INC_MULT RCONST(1000.0)
38
39 #define ZERO RCONST(0.0)
40 #define ONE RCONST(1.0)
41 #define TWO RCONST(2.0)
42
43 /*
44 * =================================================================
45 * READIBILITY REPLACEMENTS
46 * =================================================================
47 */
48
49 #define lrw1 (kin_mem->kin_lrw1)
50 #define liw1 (kin_mem->kin_liw1)
51 #define uround (kin_mem->kin_uround)
52 #define func (kin_mem->kin_func)
53 #define user_data (kin_mem->kin_user_data)
54 #define printfl (kin_mem->kin_printfl)
55 #define linit (kin_mem->kin_linit)
56 #define lsetup (kin_mem->kin_lsetup)
57 #define lsolve (kin_mem->kin_lsolve)
58 #define lfree (kin_mem->kin_lfree)
59 #define lmem (kin_mem->kin_lmem)
60 #define inexact_ls (kin_mem->kin_inexact_ls)
61 #define uu (kin_mem->kin_uu)
62 #define fval (kin_mem->kin_fval)
63 #define uscale (kin_mem->kin_uscale)
64 #define fscale (kin_mem->kin_fscale)
65 #define sqrt_relfunc (kin_mem->kin_sqrt_relfunc)
66 #define sJpnorm (kin_mem->kin_sJpnorm)
67 #define sfdotJp (kin_mem->kin_sfdotJp)
68 #define errfp (kin_mem->kin_errfp)
69 #define infofp (kin_mem->kin_infofp)
70 #define setupNonNull (kin_mem->kin_setupNonNull)
71 #define vtemp1 (kin_mem->kin_vtemp1)
72 #define vec_tmpl (kin_mem->kin_vtemp1)
73 #define vtemp2 (kin_mem->kin_vtemp2)
74
75 #define mtype (kindls_mem->d_type)
76 #define n (kindls_mem->d_n)
77 #define ml (kindls_mem->d_ml)
78 #define mu (kindls_mem->d_mu)
79 #define smu (kindls_mem->d_smu)
80 #define jacDQ (kindls_mem->d_jacDQ)
81 #define djac (kindls_mem->d_djac)
82 #define bjac (kindls_mem->d_bjac)
83 #define J (kindls_mem->d_J)
84 #define pivots (kindls_mem->d_pivots)
85 #define nje (kindls_mem->d_nje)
86 #define nfeDQ (kindls_mem->d_nfeDQ)
87 #define J_data (kindls_mem->d_J_data)
88 #define last_flag (kindls_mem->d_last_flag)
89
90 /*
91 * =================================================================
92 * EXPORTED FUNCTIONS
93 * =================================================================
94 */
95
96 /*
97 * -----------------------------------------------------------------
98 * KINDlsSetJacFn
99 * -----------------------------------------------------------------
100 */
101
KINDlsSetDenseJacFn(void * kinmem,KINDlsDenseJacFn jac)102 int KINDlsSetDenseJacFn(void *kinmem, KINDlsDenseJacFn jac)
103 {
104 KINMem kin_mem;
105 KINDlsMem kindls_mem;
106
107 /* Return immediately if kinmem is NULL */
108 if (kinmem == NULL)
109 {
110 KINProcessError(NULL, KINDLS_MEM_NULL, "KINDLS", "KINDlsSetDenseJacFn", MSGD_KINMEM_NULL);
111 return(KINDLS_MEM_NULL);
112 }
113 kin_mem = (KINMem) kinmem;
114
115 if (lmem == NULL)
116 {
117 KINProcessError(kin_mem, KINDLS_LMEM_NULL, "KINDLS", "KINDlsSetDenseJacFn", MSGD_LMEM_NULL);
118 return(KINDLS_LMEM_NULL);
119 }
120 kindls_mem = (KINDlsMem) lmem;
121
122 if (jac != NULL)
123 {
124 jacDQ = FALSE;
125 djac = jac;
126 }
127 else
128 {
129 jacDQ = TRUE;
130 }
131
132 return(KINDLS_SUCCESS);
133 }
134
KINDlsSetBandJacFn(void * kinmem,KINDlsBandJacFn jac)135 int KINDlsSetBandJacFn(void *kinmem, KINDlsBandJacFn jac)
136 {
137 KINMem kin_mem;
138 KINDlsMem kindls_mem;
139
140 /* Return immediately if kinmem is NULL */
141 if (kinmem == NULL)
142 {
143 KINProcessError(NULL, KINDLS_MEM_NULL, "KINDLS", "KINDlsSetBandJacFn", MSGD_KINMEM_NULL);
144 return(KINDLS_MEM_NULL);
145 }
146 kin_mem = (KINMem) kinmem;
147
148 if (lmem == NULL)
149 {
150 KINProcessError(kin_mem, KINDLS_LMEM_NULL, "KINDLS", "KINDlsSetBandJacFn", MSGD_LMEM_NULL);
151 return(KINDLS_LMEM_NULL);
152 }
153 kindls_mem = (KINDlsMem) lmem;
154
155 if (jac != NULL)
156 {
157 jacDQ = FALSE;
158 bjac = jac;
159 }
160 else
161 {
162 jacDQ = TRUE;
163 }
164
165 return(KINDLS_SUCCESS);
166 }
167
168 /*
169 * -----------------------------------------------------------------
170 * KINDlsGetWorkSpace
171 * -----------------------------------------------------------------
172 */
173
KINDlsGetWorkSpace(void * kinmem,long int * lenrwLS,long int * leniwLS)174 int KINDlsGetWorkSpace(void *kinmem, long int *lenrwLS, long int *leniwLS)
175 {
176 KINMem kin_mem;
177 KINDlsMem kindls_mem;
178
179 /* Return immediately if kinmem is NULL */
180 if (kinmem == NULL)
181 {
182 KINProcessError(NULL, KINDLS_MEM_NULL, "KINDLS", "KINBandGetWorkSpace", MSGD_KINMEM_NULL);
183 return(KINDLS_MEM_NULL);
184 }
185 kin_mem = (KINMem) kinmem;
186
187 if (lmem == NULL)
188 {
189 KINProcessError(kin_mem, KINDLS_LMEM_NULL, "KINDLS", "KINBandGetWorkSpace", MSGD_LMEM_NULL);
190 return(KINDLS_LMEM_NULL);
191 }
192 kindls_mem = (KINDlsMem) lmem;
193
194 if (mtype == SUNDIALS_DENSE)
195 {
196 *lenrwLS = n * n;
197 *leniwLS = n;
198 }
199 else if (mtype == SUNDIALS_BAND)
200 {
201 *lenrwLS = n * (smu + mu + 2 * ml + 2);
202 *leniwLS = n;
203 }
204
205 return(KINDLS_SUCCESS);
206 }
207
208 /*
209 * -----------------------------------------------------------------
210 * KINDlsGetNumJacEvals
211 * -----------------------------------------------------------------
212 */
213
KINDlsGetNumJacEvals(void * kinmem,long int * njevals)214 int KINDlsGetNumJacEvals(void *kinmem, long int *njevals)
215 {
216 KINMem kin_mem;
217 KINDlsMem kindls_mem;
218
219 /* Return immediately if kinmem is NULL */
220 if (kinmem == NULL)
221 {
222 KINProcessError(NULL, KINDLS_MEM_NULL, "KINDLS", "KINDlsGetNumJacEvals", MSGD_KINMEM_NULL);
223 return(KINDLS_MEM_NULL);
224 }
225 kin_mem = (KINMem) kinmem;
226
227 if (lmem == NULL)
228 {
229 KINProcessError(kin_mem, KINDLS_LMEM_NULL, "KINDLS", "KINDlsGetNumJacEvals", MSGD_LMEM_NULL);
230 return(KINDLS_LMEM_NULL);
231 }
232 kindls_mem = (KINDlsMem) lmem;
233
234 *njevals = nje;
235
236 return(KINDLS_SUCCESS);
237 }
238
239 /*
240 * -----------------------------------------------------------------
241 * KINDlsGetNumFuncEvals
242 * -----------------------------------------------------------------
243 */
244
KINDlsGetNumFuncEvals(void * kinmem,long int * nfevalsLS)245 int KINDlsGetNumFuncEvals(void *kinmem, long int *nfevalsLS)
246 {
247 KINMem kin_mem;
248 KINDlsMem kindls_mem;
249
250 /* Return immediately if kinmem is NULL */
251 if (kinmem == NULL)
252 {
253 KINProcessError(NULL, KINDLS_MEM_NULL, "KINDLS", "KINDlsGetNumFuncEvals", MSGD_KINMEM_NULL);
254 return(KINDLS_MEM_NULL);
255 }
256 kin_mem = (KINMem) kinmem;
257
258 if (lmem == NULL)
259 {
260 KINProcessError(kin_mem, KINDLS_LMEM_NULL, "KINDLS", "KINDlsGetNumGuncEvals", MSGD_LMEM_NULL);
261 return(KINDLS_LMEM_NULL);
262 }
263 kindls_mem = (KINDlsMem) lmem;
264
265 *nfevalsLS = nfeDQ;
266
267 return(KINDLS_SUCCESS);
268 }
269
270 /*
271 * -----------------------------------------------------------------
272 * KINDlsGetLastFlag
273 * -----------------------------------------------------------------
274 */
275
KINDlsGetLastFlag(void * kinmem,long int * flag)276 int KINDlsGetLastFlag(void *kinmem, long int *flag)
277 {
278 KINMem kin_mem;
279 KINDlsMem kindls_mem;
280
281 /* Return immediately if kinmem is NULL */
282 if (kinmem == NULL)
283 {
284 KINProcessError(NULL, KINDLS_MEM_NULL, "KINDLS", "KINDlsGetLastFlag", MSGD_KINMEM_NULL);
285 return(KINDLS_MEM_NULL);
286 }
287 kin_mem = (KINMem) kinmem;
288
289 if (lmem == NULL)
290 {
291 KINProcessError(kin_mem, KINDLS_LMEM_NULL, "KINDLS", "KINDlsGetLastFlag", MSGD_LMEM_NULL);
292 return(KINDLS_LMEM_NULL);
293 }
294 kindls_mem = (KINDlsMem) lmem;
295
296 *flag = last_flag;
297
298 return(KINDLS_SUCCESS);
299 }
300
301 /*
302 * -----------------------------------------------------------------
303 * KINDlsGetReturnFlagName
304 * -----------------------------------------------------------------
305 */
306
KINDlsGetReturnFlagName(long int flag)307 char *KINDlsGetReturnFlagName(long int flag)
308 {
309 char *name;
310
311 name = (char *)malloc(30 * sizeof(char));
312
313 switch (flag)
314 {
315 case KINDLS_SUCCESS:
316 sprintf(name, "KINDLS_SUCCESS");
317 break;
318 case KINDLS_MEM_NULL:
319 sprintf(name, "KINDLS_MEM_NULL");
320 break;
321 case KINDLS_LMEM_NULL:
322 sprintf(name, "KINDLS_LMEM_NULL");
323 break;
324 case KINDLS_ILL_INPUT:
325 sprintf(name, "KINDLS_ILL_INPUT");
326 break;
327 case KINDLS_MEM_FAIL:
328 sprintf(name, "KINDLS_MEM_FAIL");
329 break;
330 default:
331 sprintf(name, "NONE");
332 }
333
334 return(name);
335 }
336
337 /*
338 * =================================================================
339 * DQ JACOBIAN APPROXIMATIONS
340 * =================================================================
341 */
342
343
344
345 /*
346 * -----------------------------------------------------------------
347 * kinDlsDenseDQJac
348 * -----------------------------------------------------------------
349 * This routine generates a dense difference quotient approximation to
350 * the Jacobian of F(u). It assumes that a dense matrix of type
351 * DlsMat is stored column-wise, and that elements within each column
352 * are contiguous. The address of the jth column of J is obtained via
353 * the macro DENSE_COL and this pointer is associated with an N_Vector
354 * using the N_VGetArrayPointer/N_VSetArrayPointer functions.
355 * Finally, the actual computation of the jth column of the Jacobian is
356 * done with a call to N_VLinearSum.
357 *
358 * The increment used in the finite-difference approximation
359 * J_ij = ( F_i(u+sigma_j * e_j) - F_i(u) ) / sigma_j
360 * is
361 * sigma_j = max{|u_j|, |1/uscale_j|} * sqrt(uround)
362 *
363 * Note: uscale_j = 1/typ(u_j)
364 *
365 * NOTE: Any type of failure of the system function her leads to an
366 * unrecoverable failure of the Jacobian function and thus
367 * of the linear solver setup function, stopping KINSOL.
368 * -----------------------------------------------------------------
369 */
370
kinDlsDenseDQJac(long int N,N_Vector u,N_Vector fu,DlsMat Jac,void * data,N_Vector tmp1,N_Vector tmp2)371 int kinDlsDenseDQJac(long int N,
372 N_Vector u, N_Vector fu,
373 DlsMat Jac, void *data,
374 N_Vector tmp1, N_Vector tmp2)
375 {
376 realtype inc, inc_inv, ujsaved, ujscale, sign;
377 realtype *tmp2_data, *u_data, *uscale_data;
378 N_Vector ftemp, jthCol;
379 int retval = 0;
380 long int j;
381
382 KINMem kin_mem;
383 KINDlsMem kindls_mem;
384
385 /* data points to kin_mem */
386 kin_mem = (KINMem) data;
387 kindls_mem = (KINDlsMem) lmem;
388
389 /* Save pointer to the array in tmp2 */
390 tmp2_data = N_VGetArrayPointer(tmp2);
391
392 /* Rename work vectors for readibility */
393 ftemp = tmp1;
394 jthCol = tmp2;
395
396 /* Obtain pointers to the data for u and uscale */
397 u_data = N_VGetArrayPointer(u);
398 uscale_data = N_VGetArrayPointer(uscale);
399
400 /* This is the only for loop for 0..N-1 in KINSOL */
401
402 for (j = 0; j < N; j++)
403 {
404
405 /* Generate the jth col of Jac(u) */
406
407 N_VSetArrayPointer(DENSE_COL(Jac, j), jthCol);
408
409 ujsaved = u_data[j];
410 ujscale = ONE / uscale_data[j];
411 sign = (ujsaved >= ZERO) ? ONE : -ONE;
412 inc = sqrt_relfunc * MAX(ABS(ujsaved), ujscale) * sign;
413 u_data[j] += inc;
414
415 retval = func(u, ftemp, user_data);
416 nfeDQ++;
417 if (retval != 0)
418 {
419 break;
420 }
421
422 u_data[j] = ujsaved;
423
424 inc_inv = ONE / inc;
425 N_VLinearSum(inc_inv, ftemp, -inc_inv, fu, jthCol);
426
427 }
428
429 /* Restore original array pointer in tmp2 */
430 N_VSetArrayPointer(tmp2_data, tmp2);
431
432 return(retval);
433 }
434
435 /*
436 * -----------------------------------------------------------------
437 * kinDlsBandDQJac
438 * -----------------------------------------------------------------
439 * This routine generates a banded difference quotient approximation to
440 * the Jacobian of F(u). It assumes that a band matrix of type
441 * BandMat is stored column-wise, and that elements within each column
442 * are contiguous. This makes it possible to get the address of a column
443 * of J via the macro BAND_COL and to write a simple for loop to set
444 * each of the elements of a column in succession.
445 *
446 * NOTE: Any type of failure of the system function her leads to an
447 * unrecoverable failure of the Jacobian function and thus
448 * of the linear solver setup function, stopping KINSOL.
449 * -----------------------------------------------------------------
450 */
451
kinDlsBandDQJac(long int N,long int mupper,long int mlower,N_Vector u,N_Vector fu,DlsMat Jac,void * data,N_Vector tmp1,N_Vector tmp2)452 int kinDlsBandDQJac(long int N, long int mupper, long int mlower,
453 N_Vector u, N_Vector fu,
454 DlsMat Jac, void *data,
455 N_Vector tmp1, N_Vector tmp2)
456 {
457 realtype inc, inc_inv;
458 N_Vector futemp, utemp;
459 int retval;
460 long int group, i, j, width, ngroups, i1, i2;
461 realtype *col_j, *fu_data, *futemp_data, *u_data, *utemp_data, *uscale_data;
462
463 KINMem kin_mem;
464 KINDlsMem kindls_mem;
465
466 /* data points to kinmem */
467 kin_mem = (KINMem) data;
468 kindls_mem = (KINDlsMem) lmem;
469
470 /* Rename work vectors for use as temporary values of u and fu */
471 futemp = tmp1;
472 utemp = tmp2;
473
474 /* Obtain pointers to the data for ewt, fy, futemp, y, ytemp */
475 fu_data = N_VGetArrayPointer(fu);
476 futemp_data = N_VGetArrayPointer(futemp);
477 u_data = N_VGetArrayPointer(u);
478 uscale_data = N_VGetArrayPointer(uscale);
479 utemp_data = N_VGetArrayPointer(utemp);
480
481 /* Load utemp with u */
482 N_VScale(ONE, u, utemp);
483
484 /* Set bandwidth and number of column groups for band differencing */
485 width = mlower + mupper + 1;
486 ngroups = MIN(width, N);
487
488 for (group = 1; group <= ngroups; group++)
489 {
490
491 /* Increment all utemp components in group */
492 for (j = group - 1; j < N; j += width)
493 {
494 inc = sqrt_relfunc * MAX(ABS(u_data[j]), ABS(uscale_data[j]));
495 utemp_data[j] += inc;
496 }
497
498 /* Evaluate f with incremented u */
499 retval = func(utemp, futemp, user_data);
500 if (retval != 0)
501 {
502 return(-1);
503 }
504
505 /* Restore utemp components, then form and load difference quotients */
506 for (j = group - 1; j < N; j += width)
507 {
508 utemp_data[j] = u_data[j];
509 col_j = BAND_COL(Jac, j);
510 inc = sqrt_relfunc * MAX(ABS(u_data[j]), ABS(uscale_data[j]));
511 inc_inv = ONE / inc;
512 i1 = MAX(0, j - mupper);
513 i2 = MIN(j + mlower, N - 1);
514 for (i = i1; i <= i2; i++)
515 {
516 BAND_COL_ELEM(col_j, i, j) = inc_inv * (futemp_data[i] - fu_data[i]);
517 }
518 }
519 }
520
521 /* Increment counter nfeDQ */
522 nfeDQ += ngroups;
523
524 return(0);
525 }
526