1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 4075 $
4 * $Date: 2014-04-24 10:46:58 -0700 (Thu, 24 Apr 2014) $
5 * -----------------------------------------------------------------
6 * Programmer(s): Radu Serban and Aaron Collier @ 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 KINSOL scaled,
19 * preconditioned GMRES linear solver, KINSpgmr.
20 * -----------------------------------------------------------------
21 */
22
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <stdarg.h>
26
27 #include "kinsol_impl.h"
28 #include <kinsol/kinsol_spgmr.h>
29 #include "kinsol_spils_impl.h"
30
31 #include <sundials/sundials_spgmr.h>
32 #include <sundials/sundials_math.h>
33
34 /*
35 * -----------------------------------------------------------------
36 * private constants
37 * -----------------------------------------------------------------
38 */
39
40 #define ZERO RCONST(0.0)
41
42 /*
43 * -----------------------------------------------------------------
44 * function prototypes
45 * -----------------------------------------------------------------
46 */
47
48 /* KINSpgmr linit, lsetup, lsolve, and lfree routines */
49
50 static int KINSpgmrInit(KINMem kin_mem);
51 static int KINSpgmrSetup(KINMem kin_mem);
52 static int KINSpgmrSolve(KINMem kin_mem, N_Vector xx,
53 N_Vector bb, realtype *sJpnorm, realtype *sFdotJp);
54 static void KINSpgmrFree(KINMem kin_mem);
55
56 /*
57 * -----------------------------------------------------------------
58 * readability replacements
59 * -----------------------------------------------------------------
60 */
61
62 #define nni (kin_mem->kin_nni)
63 #define nnilset (kin_mem->kin_nnilset)
64 #define func (kin_mem->kin_func)
65 #define user_data (kin_mem->kin_user_data)
66 #define printfl (kin_mem->kin_printfl)
67 #define linit (kin_mem->kin_linit)
68 #define lsetup (kin_mem->kin_lsetup)
69 #define lsolve (kin_mem->kin_lsolve)
70 #define lfree (kin_mem->kin_lfree)
71 #define lmem (kin_mem->kin_lmem)
72 #define inexact_ls (kin_mem->kin_inexact_ls)
73 #define uu (kin_mem->kin_uu)
74 #define fval (kin_mem->kin_fval)
75 #define uscale (kin_mem->kin_uscale)
76 #define fscale (kin_mem->kin_fscale)
77 #define sqrt_relfunc (kin_mem->kin_sqrt_relfunc)
78 #define jacCurrent (kin_mem->kin_jacCurrent)
79 #define eps (kin_mem->kin_eps)
80 #define errfp (kin_mem->kin_errfp)
81 #define infofp (kin_mem->kin_infofp)
82 #define setupNonNull (kin_mem->kin_setupNonNull)
83 #define vtemp1 (kin_mem->kin_vtemp1)
84 #define vec_tmpl (kin_mem->kin_vtemp1)
85 #define vtemp2 (kin_mem->kin_vtemp2)
86 #define strategy (kin_mem->kin_globalstrategy)
87
88 #define pretype (kinspils_mem->s_pretype)
89 #define gstype (kinspils_mem->s_gstype)
90 #define nli (kinspils_mem->s_nli)
91 #define npe (kinspils_mem->s_npe)
92 #define nps (kinspils_mem->s_nps)
93 #define ncfl (kinspils_mem->s_ncfl)
94 #define njtimes (kinspils_mem->s_njtimes)
95 #define nfes (kinspils_mem->s_nfes)
96 #define new_uu (kinspils_mem->s_new_uu)
97 #define spils_mem (kinspils_mem->s_spils_mem)
98
99 #define jtimesDQ (kinspils_mem->s_jtimesDQ)
100 #define jtimes (kinspils_mem->s_jtimes)
101 #define J_data (kinspils_mem->s_J_data)
102
103 #define last_flag (kinspils_mem->s_last_flag)
104
105 /*
106 * -----------------------------------------------------------------
107 * Function : KINSpgmr
108 * -----------------------------------------------------------------
109 * This routine allocates and initializes the memory record and
110 * sets function fields specific to the SPGMR linear solver module.
111 * KINSpgmr sets the kin_linit, kin_lsetup, kin_lsolve, and
112 * kin_lfree fields in *kinmem to be KINSpgmrInit, KINSpgmrSetup,
113 * KINSpgmrSolve, and KINSpgmrFree, respectively. It allocates
114 * memory for a structure of type KINSpilsMemRec and sets the
115 * kin_lmem field in *kinmem to the address of this structure. It
116 * also calls SpgmrMalloc to allocate memory for the module
117 * SPGMR. In summary, KINSpgmr sets various fields in the
118 * KINSpilsMemRec structure.
119 * -----------------------------------------------------------------
120 */
121
KINSpgmr(void * kinmem,int maxl)122 int KINSpgmr(void *kinmem, int maxl)
123 {
124 KINMem kin_mem;
125 KINSpilsMem kinspils_mem;
126 SpgmrMem spgmr_mem;
127 int maxl1;
128
129 if (kinmem == NULL){
130 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpgmr", MSGS_KINMEM_NULL);
131 return(KINSPILS_MEM_NULL);
132 }
133 kin_mem = (KINMem) kinmem;
134
135 /* check for required vector operations */
136
137 /* Note: do NOT need to check for N_VLinearSum, N_VProd, N_VScale, N_VDiv,
138 or N_VWL2Norm because they are required by KINSOL */
139
140 if ((vec_tmpl->ops->nvconst == NULL) ||
141 (vec_tmpl->ops->nvdotprod == NULL) ||
142 (vec_tmpl->ops->nvl1norm == NULL)) {
143 KINProcessError(NULL, KINSPILS_ILL_INPUT, "KINSPILS", "KINSpgmr", MSGS_BAD_NVECTOR);
144 return(KINSPILS_ILL_INPUT);
145 }
146
147 if (lfree != NULL) lfree(kin_mem);
148
149 /* set four main function fields in kin_mem */
150
151 linit = KINSpgmrInit;
152 lsetup = KINSpgmrSetup;
153 lsolve = KINSpgmrSolve;
154 lfree = KINSpgmrFree;
155
156 /* get memory for KINSpilsMemRec */
157 kinspils_mem = NULL;
158 kinspils_mem = (KINSpilsMem) malloc(sizeof(struct KINSpilsMemRec));
159 if (kinspils_mem == NULL){
160 KINProcessError(NULL, KINSPILS_MEM_FAIL, "KINSPILS", "KINSpgmr", MSGS_MEM_FAIL);
161 return(KINSPILS_MEM_FAIL);
162 }
163
164 /* Set ILS type */
165 kinspils_mem->s_type = SPILS_SPGMR;
166
167 /* set SPGMR parameters that were passed in call sequence */
168
169 maxl1 = (maxl <= 0) ? KINSPILS_MAXL : maxl;
170 kinspils_mem->s_maxl = maxl1;
171
172 /* Set defaults for Jacobian-related fileds */
173
174 jtimesDQ = TRUE;
175 jtimes = NULL;
176 J_data = NULL;
177
178 /* Set defaults for preconditioner-related fields */
179
180 kinspils_mem->s_pset = NULL;
181 kinspils_mem->s_psolve = NULL;
182 kinspils_mem->s_pfree = NULL;
183 kinspils_mem->s_P_data = kin_mem->kin_user_data;
184
185 /* Set default values for the rest of the SPGMR parameters */
186
187 kinspils_mem->s_pretype = PREC_NONE;
188 kinspils_mem->s_gstype = MODIFIED_GS;
189 kinspils_mem->s_maxlrst = 0;
190 kinspils_mem->s_last_flag = KINSPILS_SUCCESS;
191
192 /* Call SpgmrMalloc to allocate workspace for SPGMR */
193
194 /* vec_tmpl passed as template vector */
195 spgmr_mem = NULL;
196 spgmr_mem = SpgmrMalloc(maxl1, vec_tmpl);
197 if (spgmr_mem == NULL) {
198 KINProcessError(NULL, KINSPILS_MEM_FAIL, "KINSPILS", "KINSpgmr", MSGS_MEM_FAIL);
199 free(kinspils_mem); kinspils_mem = NULL;
200 return(KINSPILS_MEM_FAIL);
201 }
202
203 /* This is an iterative linear solver */
204
205 inexact_ls = TRUE;
206
207 /* Attach SPGMR memory to spils memory structure */
208 spils_mem = (void *) spgmr_mem;
209
210 /* attach linear solver memory to KINSOL memory */
211 lmem = kinspils_mem;
212
213 return(KINSPILS_SUCCESS);
214 }
215
216 /*
217 * -----------------------------------------------------------------
218 * additional readability replacements
219 * -----------------------------------------------------------------
220 */
221
222 #define maxl (kinspils_mem->s_maxl)
223 #define maxlrst (kinspils_mem->s_maxlrst)
224 #define pset (kinspils_mem->s_pset)
225 #define psolve (kinspils_mem->s_psolve)
226 #define P_data (kinspils_mem->s_P_data)
227
228 /*
229 * -----------------------------------------------------------------
230 * Function : KINSpgmrInit
231 * -----------------------------------------------------------------
232 * This routine initializes variables associated with the GMRES
233 * linear solver. Memory allocation was done previously in
234 * KINSpgmr.
235 * -----------------------------------------------------------------
236 */
237
KINSpgmrInit(KINMem kin_mem)238 static int KINSpgmrInit(KINMem kin_mem)
239 {
240 KINSpilsMem kinspils_mem;
241
242 kinspils_mem = (KINSpilsMem) lmem;
243
244 /* initialize counters */
245
246 npe = nli = nps = ncfl = 0;
247 njtimes = nfes = 0;
248
249 /* set preconditioner type */
250
251 if (psolve != NULL) {
252 pretype = PREC_RIGHT;
253 } else {
254 pretype = PREC_NONE;
255 }
256
257 /* set setupNonNull to TRUE iff there is preconditioning with setup */
258
259 setupNonNull = (psolve != NULL) && (pset != NULL);
260
261 /* Set Jacobian-related fields, based on jtimesDQ */
262
263 if (jtimesDQ) {
264 jtimes = KINSpilsDQJtimes;
265 J_data = kin_mem;
266 } else {
267 J_data = user_data;
268 }
269
270 if ( (strategy == KIN_PICARD) && jtimesDQ ) {
271 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSpgmrInit",
272 MSG_NOL_FAIL);
273 return(KIN_ILL_INPUT);
274 }
275
276 last_flag = KINSPILS_SUCCESS;
277 return(0);
278 }
279
280 /*
281 * -----------------------------------------------------------------
282 * Function : KINSpgmrSetup
283 * -----------------------------------------------------------------
284 * This routine does the setup operations for the SPGMR linear
285 * solver, that is, it is an interface to the user-supplied
286 * routine pset.
287 * -----------------------------------------------------------------
288 */
289
KINSpgmrSetup(KINMem kin_mem)290 static int KINSpgmrSetup(KINMem kin_mem)
291 {
292 KINSpilsMem kinspils_mem;
293 int ret;
294
295 kinspils_mem = (KINSpilsMem) lmem;
296
297 /* call pset routine */
298
299 ret = pset(uu, uscale, fval, fscale, P_data, vtemp1, vtemp2);
300
301 last_flag = ret;
302
303 npe++;
304 nnilset = nni;
305
306 /* return the same value ret that pset returned */
307
308 return(ret);
309 }
310
311 /*
312 * -----------------------------------------------------------------
313 * Function : KINSpgmrSolve
314 * -----------------------------------------------------------------
315 * This routine handles the call to the generic SPGMR solver
316 * SpgmrSolve for the solution of the linear system Ax = b.
317 *
318 * Appropriate variables are passed to SpgmrSolve and the counters
319 * nli, nps, and ncfl are incremented, and the return value is set
320 * according to the success of SpgmrSolve. The success flag is
321 * returned if SpgmrSolve converged, or if the residual was reduced.
322 * Of the other error conditions, only preconditioner solver
323 * failure is specifically returned. Otherwise a generic flag is
324 * returned to denote failure of this routine.
325 * -----------------------------------------------------------------
326 */
327
KINSpgmrSolve(KINMem kin_mem,N_Vector xx,N_Vector bb,realtype * sJpnorm,realtype * sFdotJp)328 static int KINSpgmrSolve(KINMem kin_mem, N_Vector xx, N_Vector bb,
329 realtype *sJpnorm, realtype *sFdotJp)
330 {
331 KINSpilsMem kinspils_mem;
332 SpgmrMem spgmr_mem;
333 int ret, nli_inc, nps_inc;
334 realtype res_norm;
335
336 kinspils_mem = (KINSpilsMem) lmem;
337
338 spgmr_mem = (SpgmrMem) spils_mem;
339
340 /* Set initial guess to xx = 0. bb is set, by the routine
341 calling KINSpgmrSolve, to the RHS vector for the system
342 to be solved. */
343
344 N_VConst(ZERO, xx);
345
346 new_uu = TRUE; /* set flag required for user Jacobian routine */
347
348 /* call SpgmrSolve */
349
350 ret = SpgmrSolve(spgmr_mem, kin_mem, xx, bb, pretype, gstype, eps,
351 maxlrst, kin_mem, fscale, fscale, KINSpilsAtimes,
352 KINSpilsPSolve, &res_norm, &nli_inc, &nps_inc);
353
354 /* increment counters nli, nps, and ncfl
355 (nni is updated in the KINSol main iteration loop) */
356
357 nli = nli + (long int) nli_inc;
358 nps = nps + (long int) nps_inc;
359
360 if (printfl > 2)
361 KINPrintInfo(kin_mem, PRNT_NLI, "KINSPGMR", "KINSpgmrSolve", INFO_NLI, nli_inc);
362
363 if (ret != 0) ncfl++;
364 last_flag = ret;
365
366 if ( (ret != 0) && (ret != SPGMR_RES_REDUCED) ) {
367
368 /* Handle all failure returns from SpgmrSolve */
369
370 switch(ret) {
371 case SPGMR_PSOLVE_FAIL_REC:
372 case SPGMR_ATIMES_FAIL_REC:
373 return(1);
374 break;
375 case SPGMR_CONV_FAIL:
376 case SPGMR_QRFACT_FAIL:
377 case SPGMR_MEM_NULL:
378 case SPGMR_GS_FAIL:
379 case SPGMR_QRSOL_FAIL:
380 case SPGMR_ATIMES_FAIL_UNREC:
381 case SPGMR_PSOLVE_FAIL_UNREC:
382 return(-1);
383 break;
384 }
385 }
386
387 /* SpgmrSolve returned either SPGMR_SUCCESS or SPGMR_RES_REDUCED.
388
389 Compute the terms sJpnorm and sFdotJp for use in the linesearch
390 routine and in KINForcingTerm. Both of these terms are subsequently
391 corrected if the step is reduced by constraints or the linesearch.
392
393 sJpnorm is the norm of the scaled product (scaled by fscale) of the
394 current Jacobian matrix J and the step vector p (= solution vector xx).
395
396 sFdotJp is the dot product of the scaled f vector and the scaled
397 vector J*p, where the scaling uses fscale. */
398
399 ret = KINSpilsAtimes(kin_mem, xx, bb);
400 if (ret > 0) {
401 last_flag = SPGMR_ATIMES_FAIL_REC;
402 return(1);
403 }
404 else if (ret < 0) {
405 last_flag = SPGMR_ATIMES_FAIL_UNREC;
406 return(-1);
407 }
408
409 *sJpnorm = N_VWL2Norm(bb, fscale);
410 N_VProd(bb, fscale, bb);
411 N_VProd(bb, fscale, bb);
412 *sFdotJp = N_VDotProd(fval, bb);
413
414 if (printfl > 2) KINPrintInfo(kin_mem, PRNT_EPS, "KINSPGMR",
415 "KINSpgmrSolve", INFO_EPS, res_norm, eps);
416
417 return(0);
418 }
419
420 /*
421 * -----------------------------------------------------------------
422 * Function : KINSpgmrFree
423 * -----------------------------------------------------------------
424 * This routine frees memory specific to the SPGMR linear solver.
425 * -----------------------------------------------------------------
426 */
427
KINSpgmrFree(KINMem kin_mem)428 static void KINSpgmrFree(KINMem kin_mem)
429 {
430 KINSpilsMem kinspils_mem;
431 SpgmrMem spgmr_mem;
432
433 kinspils_mem = (KINSpilsMem) lmem;
434
435 spgmr_mem = (SpgmrMem) spils_mem;
436 SpgmrFree(spgmr_mem);
437
438 if (kinspils_mem->s_pfree != NULL) (kinspils_mem->s_pfree)(kin_mem);
439
440 free(kinspils_mem); kinspils_mem = NULL;
441 }
442