1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 4397 $
4 * $Date: 2015-02-28 14:03:10 -0800 (Sat, 28 Feb 2015) $
5 * -----------------------------------------------------------------
6 * Programmer(s): Aaron Collier and 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 KINSOL interface to the
19 * scaled, preconditioned TFQMR (SPTFQMR) iterative linear solver.
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_sptfqmr.h>
29 #include "kinsol_spils_impl.h"
30
31 #include <sundials/sundials_sptfqmr.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 /* KINSptfqmr linit, lsetup, lsolve, and lfree routines */
49
50 static int KINSptfqmrInit(KINMem kin_mem);
51 static int KINSptfqmrSetup(KINMem kin_mem);
52 static int KINSptfqmrSolve(KINMem kin_mem, N_Vector xx,
53 N_Vector bb, realtype *sJpnorm, realtype *sFdotJp);
54 static void KINSptfqmrFree(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 eps (kin_mem->kin_eps)
79 #define errfp (kin_mem->kin_errfp)
80 #define infofp (kin_mem->kin_infofp)
81 #define setupNonNull (kin_mem->kin_setupNonNull)
82 #define vtemp1 (kin_mem->kin_vtemp1)
83 #define vec_tmpl (kin_mem->kin_vtemp1)
84 #define vtemp2 (kin_mem->kin_vtemp2)
85 #define strategy (kin_mem->kin_globalstrategy)
86
87 #define pretype (kinspils_mem->s_pretype)
88 #define nli (kinspils_mem->s_nli)
89 #define npe (kinspils_mem->s_npe)
90 #define nps (kinspils_mem->s_nps)
91 #define ncfl (kinspils_mem->s_ncfl)
92 #define njtimes (kinspils_mem->s_njtimes)
93 #define nfes (kinspils_mem->s_nfes)
94 #define new_uu (kinspils_mem->s_new_uu)
95 #define spils_mem (kinspils_mem->s_spils_mem)
96
97 #define jtimesDQ (kinspils_mem->s_jtimesDQ)
98 #define jtimes (kinspils_mem->s_jtimes)
99 #define J_data (kinspils_mem->s_J_data)
100
101 #define last_flag (kinspils_mem->s_last_flag)
102
103 /*
104 * -----------------------------------------------------------------
105 * Function : KINSptfqmr
106 * -----------------------------------------------------------------
107 * This routine allocates and initializes the memory record and
108 * sets function fields specific to the SPTFQMR linear solver module.
109 * KINSptfqmr sets the kin_linit, kin_lsetup, kin_lsolve, and
110 * kin_lfree fields in *kinmem to be KINSptfqmrInit, KINSptfqmrSetup,
111 * KINSptfqmrSolve, and KINSptfqmrFree, respectively. It allocates
112 * memory for a structure of type KINSpilsMemRec and sets the
113 * kin_lmem field in *kinmem to the address of this structure. It
114 * also calls SptfqmrMalloc to allocate memory for the module
115 * SPTFQMR. It sets setupNonNull in (*kin_mem),
116 * and sets various fields in the KINSpilsMemRec structure.
117 * Finally, KINSptfqmr allocates memory for local vectors, and calls
118 * SptfqmrMalloc to allocate memory for the Sptfqmr solver.
119 * -----------------------------------------------------------------
120 */
121
KINSptfqmr(void * kinmem,int maxl)122 int KINSptfqmr(void *kinmem, int maxl)
123 {
124 KINMem kin_mem;
125 KINSpilsMem kinspils_mem;
126 SptfqmrMem sptfqmr_mem;
127 int maxl1;
128
129 if (kinmem == NULL){
130 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSptfqmr", 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", "KINSptfqmr", 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 = KINSptfqmrInit;
152 lsetup = KINSptfqmrSetup;
153 lsolve = KINSptfqmrSolve;
154 lfree = KINSptfqmrFree;
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", "KINSptfqmr", MSGS_MEM_FAIL);
161 return(KINSPILS_MEM_FAIL);
162 }
163
164 /* Set ILS type */
165 kinspils_mem->s_type = SPILS_SPTFQMR;
166
167 /* set SPTFQMR 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 SPTFQMR parameters */
186
187 kinspils_mem->s_pretype = PREC_NONE;
188 kinspils_mem->s_last_flag = KINSPILS_SUCCESS;
189
190 /* Call SptfqmrMalloc to allocate workspace for SPTFQMR */
191
192 /* vec_tmpl passed as template vector */
193
194 sptfqmr_mem = NULL;
195 sptfqmr_mem = SptfqmrMalloc(maxl1, vec_tmpl);
196 if (sptfqmr_mem == NULL) {
197 KINProcessError(NULL, KINSPILS_MEM_FAIL, "KINSPILS", "KINSptfqmr", MSGS_MEM_FAIL);
198 free(kinspils_mem); kinspils_mem = NULL;
199 return(KINSPILS_MEM_FAIL);
200 }
201
202 /* this is an iterative linear solver */
203
204 inexact_ls = TRUE;
205
206 /* Attach SPTFQMR memory to spils memory structure */
207 spils_mem = (void *) sptfqmr_mem;
208
209 /* attach linear solver memory to KINSOL memory */
210 lmem = kinspils_mem;
211
212 return(KINSPILS_SUCCESS);
213 }
214
215 /*
216 * -----------------------------------------------------------------
217 * additional readability replacements
218 * -----------------------------------------------------------------
219 */
220
221 #define maxl (kinspils_mem->s_maxl)
222 #define pset (kinspils_mem->s_pset)
223 #define psolve (kinspils_mem->s_psolve)
224 #define P_data (kinspils_mem->s_P_data)
225
226 /*
227 * -----------------------------------------------------------------
228 * Function : KINSptfqmrInit
229 * -----------------------------------------------------------------
230 * This routine initializes variables associated with the SPTFQMR
231 * iterative linear solver. Memory allocation was done previously
232 * in KINSptfqmr.
233 * -----------------------------------------------------------------
234 */
235
KINSptfqmrInit(KINMem kin_mem)236 static int KINSptfqmrInit(KINMem kin_mem)
237 {
238 KINSpilsMem kinspils_mem;
239 kinspils_mem = (KINSpilsMem) lmem;
240
241 /* initialize counters */
242
243 npe = nli = nps = ncfl = 0;
244 njtimes = nfes = 0;
245
246 /* set preconditioner type */
247
248 if (psolve != NULL) {
249 pretype = PREC_RIGHT;
250 } else {
251 pretype = PREC_NONE;
252 }
253
254 /* set setupNonNull to TRUE iff there is preconditioning with setup */
255
256 setupNonNull = ((psolve != NULL) && (pset != NULL));
257
258 /* Set Jacobian-related fields, based on jtimesDQ */
259
260 if (jtimesDQ) {
261 jtimes = KINSpilsDQJtimes;
262 J_data = kin_mem;
263 } else {
264 J_data = user_data;
265 }
266
267 if ( (strategy == KIN_PICARD) && jtimesDQ ) {
268 KINProcessError(kin_mem, KIN_ILL_INPUT, "KINSOL", "KINSptfqmrInit",
269 MSG_NOL_FAIL);
270 return(KIN_ILL_INPUT);
271 }
272
273 last_flag = KINSPILS_SUCCESS;
274 return(0);
275 }
276
277 /*
278 * -----------------------------------------------------------------
279 * Function : KINSptfqmrSetup
280 * -----------------------------------------------------------------
281 * This routine does the setup operations for the SPTFQMR linear
282 * solver, that is, it is an interface to the user-supplied
283 * routine pset.
284 * -----------------------------------------------------------------
285 */
286
KINSptfqmrSetup(KINMem kin_mem)287 static int KINSptfqmrSetup(KINMem kin_mem)
288 {
289 KINSpilsMem kinspils_mem;
290 int ret;
291
292 kinspils_mem = (KINSpilsMem) lmem;
293
294 /* call pset routine */
295
296 ret = pset(uu, uscale, fval, fscale, P_data, vtemp1, vtemp2);
297
298 last_flag = ret;
299
300 npe++;
301 nnilset = nni;
302
303 /* return the same value ret that pset returned */
304
305 return(ret);
306 }
307
308 /*
309 * -----------------------------------------------------------------
310 * Function : KINSptfqmrSolve
311 * -----------------------------------------------------------------
312 * This routine handles the call to the generic SPTFQMR solver routine
313 * called SptfqmrSolve for the solution of the linear system Ax = b.
314 *
315 * Appropriate variables are passed to SptfqmrSolve and the counters
316 * nli, nps, and ncfl are incremented, and the return value is set
317 * according to the success of SptfqmrSolve. The success flag is
318 * returned if SptfqmrSolve converged, or if the residual was reduced.
319 * Of the other error conditions, only preconditioner solver
320 * failure is specifically returned. Otherwise a generic flag is
321 * returned to denote failure of this routine.
322 * -----------------------------------------------------------------
323 */
324
KINSptfqmrSolve(KINMem kin_mem,N_Vector xx,N_Vector bb,realtype * sJpnorm,realtype * sFdotJp)325 static int KINSptfqmrSolve(KINMem kin_mem, N_Vector xx, N_Vector bb,
326 realtype *sJpnorm, realtype *sFdotJp)
327 {
328 KINSpilsMem kinspils_mem;
329 SptfqmrMem sptfqmr_mem;
330 int ret, nli_inc, nps_inc;
331 realtype res_norm;
332
333 kinspils_mem = (KINSpilsMem) lmem;
334 sptfqmr_mem = (SptfqmrMem) spils_mem;
335
336 /* Set initial guess to xx = 0. bb is set, by the routine
337 calling KINSptfqmrSolve, to the RHS vector for the system
338 to be solved. */
339
340 N_VConst(ZERO, xx);
341
342 new_uu = TRUE; /* set flag required for user Jacobian routine */
343
344 /* call SptfqmrSolve */
345
346 ret = SptfqmrSolve(sptfqmr_mem, kin_mem, xx, bb, pretype, eps,
347 kin_mem, fscale, fscale, KINSpilsAtimes,
348 KINSpilsPSolve, &res_norm, &nli_inc, &nps_inc);
349
350 /* increment counters nli, nps, and ncfl
351 (nni is updated in the KINSol main iteration loop) */
352
353 nli = nli + (long int) nli_inc;
354 nps = nps + (long int) nps_inc;
355
356 if (printfl > 2)
357 KINPrintInfo(kin_mem, PRNT_NLI, "KINSPTFQMR", "KINSptfqmrSolve", INFO_NLI, nli_inc);
358
359 if (ret != 0) ncfl++;
360 last_flag = ret;
361
362 if ( (ret != 0) && (ret != SPTFQMR_RES_REDUCED) ) {
363
364 /* Handle all failure returns from SptfqmrSolve */
365
366 switch(ret) {
367 case SPTFQMR_PSOLVE_FAIL_REC:
368 case SPTFQMR_ATIMES_FAIL_REC:
369 return(1);
370 break;
371 case SPTFQMR_CONV_FAIL:
372 case SPTFQMR_MEM_NULL:
373 case SPTFQMR_ATIMES_FAIL_UNREC:
374 case SPTFQMR_PSOLVE_FAIL_UNREC:
375 return(-1);
376 break;
377 }
378 }
379
380 /* SptfqmrSolve returned either SPTFQMR_SUCCESS or SPTFQMR_RES_REDUCED.
381
382 Compute the terms sJpnorm and sFdotJp for use in the linesearch
383 routine and in KINForcingTerm. Both of these terms are subsequently
384 corrected if the step is reduced by constraints or the linesearch.
385
386 sJpnorm is the norm of the scaled product (scaled by fscale) of the
387 current Jacobian matrix J and the step vector p (= solution vector xx).
388
389 sFdotJp is the dot product of the scaled f vector and the scaled
390 vector J*p, where the scaling uses fscale. */
391
392 ret = KINSpilsAtimes(kin_mem, xx, bb);
393 if (ret > 0) {
394 last_flag = SPTFQMR_ATIMES_FAIL_REC;
395 return(1);
396 }
397 else if (ret < 0) {
398 last_flag = SPTFQMR_ATIMES_FAIL_UNREC;
399 return(-1);
400 }
401
402 *sJpnorm = N_VWL2Norm(bb, fscale);
403 N_VProd(bb, fscale, bb);
404 N_VProd(bb, fscale, bb);
405 *sFdotJp = N_VDotProd(fval, bb);
406
407 if (printfl > 2) KINPrintInfo(kin_mem, PRNT_EPS, "KINSPTFQMR",
408 "KINSptfqmrSolve", INFO_EPS, res_norm, eps);
409
410 return(0);
411 }
412
413 /*
414 * -----------------------------------------------------------------
415 * Function : KINSptfqmrFree
416 * -----------------------------------------------------------------
417 * Frees memory specific to the SPTFQMR linear solver module.
418 * -----------------------------------------------------------------
419 */
420
KINSptfqmrFree(KINMem kin_mem)421 static void KINSptfqmrFree(KINMem kin_mem)
422 {
423 KINSpilsMem kinspils_mem;
424 SptfqmrMem sptfqmr_mem;
425
426 kinspils_mem = (KINSpilsMem) lmem;
427
428 sptfqmr_mem = (SptfqmrMem) spils_mem;
429 SptfqmrFree(sptfqmr_mem);
430
431 if (kinspils_mem->s_pfree != NULL) (kinspils_mem->s_pfree)(kin_mem);
432
433 free(kinspils_mem); kinspils_mem = NULL;
434 }
435