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