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 Bi-CGSTAB (SPBCG) 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_spbcgs.h>
29 #include "kinsol_spils_impl.h"
30
31 #include <sundials/sundials_spbcgs.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 /* KINSpbcg linit, lsetup, lsolve, and lfree routines */
49
50 static int KINSpbcgInit(KINMem kin_mem);
51 static int KINSpbcgSetup(KINMem kin_mem);
52 static int KINSpbcgSolve(KINMem kin_mem, N_Vector xx,
53 N_Vector bb, realtype *sJpnorm, realtype *sFdotJp);
54 static void KINSpbcgFree(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 : KINSpbcg
106 * -----------------------------------------------------------------
107 * This routine allocates and initializes the memory record and
108 * sets function fields specific to the SPBCG linear solver module.
109 * KINSpbcg sets the kin_linit, kin_lsetup, kin_lsolve, and
110 * kin_lfree fields in *kinmem to be KINSpbcgInit, KINSpbcgSetup,
111 * KINSpbcgSolve, and KINSpbcgFree, 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 SpbcgMalloc to allocate memory for the module
115 * SPBCG. It sets setupNonNull in (*kin_mem) and sets various
116 * fields in the KINSpilsMemRec structure.
117 * Finally, KINSpbcg allocates memory for local vectors, and calls
118 * SpbcgMalloc to allocate memory for the Spbcg solver.
119 * -----------------------------------------------------------------
120 */
121
KINSpbcg(void * kinmem,int maxl)122 int KINSpbcg(void *kinmem, int maxl)
123 {
124 KINMem kin_mem;
125 KINSpilsMem kinspils_mem;
126 SpbcgMem spbcg_mem;
127 int maxl1;
128
129 if (kinmem == NULL){
130 KINProcessError(NULL, KINSPILS_MEM_NULL, "KINSPILS", "KINSpbcg", 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", "KINSpbcg", 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 = KINSpbcgInit;
152 lsetup = KINSpbcgSetup;
153 lsolve = KINSpbcgSolve;
154 lfree = KINSpbcgFree;
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", "KINSpbcg", MSGS_MEM_FAIL);
161 return(KINSPILS_MEM_FAIL);
162 }
163
164 /* Set ILS type */
165 kinspils_mem->s_type = SPILS_SPBCG;
166
167 /* set SPBCG 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 SPBCG parameters */
186
187 kinspils_mem->s_pretype = PREC_NONE;
188 kinspils_mem->s_last_flag = KINSPILS_SUCCESS;
189
190 /* Call SpbcgMalloc to allocate workspace for SPBCG */
191
192 /* vec_tmpl passed as template vector */
193 spbcg_mem = NULL;
194 spbcg_mem = SpbcgMalloc(maxl1, vec_tmpl);
195 if (spbcg_mem == NULL) {
196 KINProcessError(NULL, KINSPILS_MEM_FAIL, "KINSPILS", "KINSpbcg", MSGS_MEM_FAIL);
197 free(kinspils_mem); kinspils_mem = NULL;
198 return(KINSPILS_MEM_FAIL);
199 }
200
201 /* This is an iterative linear solver */
202
203 inexact_ls = TRUE;
204
205 /* Attach SPBCG memory to spils memory structure */
206 spils_mem = (void *) spbcg_mem;
207
208 /* attach linear solver memory to KINSOL memory */
209 lmem = kinspils_mem;
210
211 return(KINSPILS_SUCCESS);
212 }
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 : KINSpbcgInit
229 * -----------------------------------------------------------------
230 * This routine initializes variables associated with the SPBCG
231 * iterative linear solver. Mmemory allocation was done previously
232 * in KINSpbcg.
233 * -----------------------------------------------------------------
234 */
235
KINSpbcgInit(KINMem kin_mem)236 static int KINSpbcgInit(KINMem kin_mem)
237 {
238 KINSpilsMem kinspils_mem;
239 SpbcgMem spbcg_mem;
240
241 kinspils_mem = (KINSpilsMem) lmem;
242 spbcg_mem = (SpbcgMem) spils_mem;
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", "KINSpbcgInit",
272 MSG_NOL_FAIL);
273 return(KIN_ILL_INPUT);
274 }
275
276 /* Set maxl in the SPBCG memory in case it was changed by the user */
277 spbcg_mem->l_max = maxl;
278
279 last_flag = KINSPILS_SUCCESS;
280 return(0);
281 }
282
283 /*
284 * -----------------------------------------------------------------
285 * Function : KINSpbcgSetup
286 * -----------------------------------------------------------------
287 * This routine does the setup operations for the SPBCG linear
288 * solver, that is, it is an interface to the user-supplied
289 * routine pset.
290 * -----------------------------------------------------------------
291 */
292
KINSpbcgSetup(KINMem kin_mem)293 static int KINSpbcgSetup(KINMem kin_mem)
294 {
295 KINSpilsMem kinspils_mem;
296 int ret;
297
298 kinspils_mem = (KINSpilsMem) lmem;
299
300 /* call pset routine */
301
302 ret = pset(uu, uscale, fval, fscale, P_data, vtemp1, vtemp2);
303
304 last_flag = ret;
305
306 npe++;
307 nnilset = nni;
308
309 /* return the same value ret that pset returned */
310
311 return(ret);
312 }
313
314 /*
315 * -----------------------------------------------------------------
316 * Function : KINSpbcgSolve
317 * -----------------------------------------------------------------
318 * This routine handles the call to the generic SPBCG solver routine
319 * called SpbcgSolve for the solution of the linear system Ax = b.
320 *
321 * Appropriate variables are passed to SpbcgSolve and the counters
322 * nli, nps, and ncfl are incremented, and the return value is set
323 * according to the success of SpbcgSolve. The success flag is
324 * returned if SpbcgSolve converged, or if the residual was reduced.
325 * Of the other error conditions, only preconditioner solver
326 * failure is specifically returned. Otherwise a generic flag is
327 * returned to denote failure of this routine.
328 * -----------------------------------------------------------------
329 */
330
KINSpbcgSolve(KINMem kin_mem,N_Vector xx,N_Vector bb,realtype * sJpnorm,realtype * sFdotJp)331 static int KINSpbcgSolve(KINMem kin_mem, N_Vector xx, N_Vector bb,
332 realtype *sJpnorm, realtype *sFdotJp)
333 {
334 KINSpilsMem kinspils_mem;
335 SpbcgMem spbcg_mem;
336 int ret, nli_inc, nps_inc;
337 realtype res_norm;
338
339 kinspils_mem = (KINSpilsMem) lmem;
340 spbcg_mem = (SpbcgMem) spils_mem;
341
342 /* Set initial guess to xx = 0. bb is set, by the routine
343 calling KINSpbcgSolve, to the RHS vector for the system
344 to be solved. */
345
346 N_VConst(ZERO, xx);
347
348 new_uu = TRUE; /* set flag required for user Jacobian routine */
349
350 /* call SpbcgSolve */
351
352 ret = SpbcgSolve(spbcg_mem, kin_mem, xx, bb, pretype, eps,
353 kin_mem, fscale, fscale, KINSpilsAtimes,
354 KINSpilsPSolve, &res_norm, &nli_inc, &nps_inc);
355
356 /* increment counters nli, nps, and ncfl
357 (nni is updated in the KINSol main iteration loop) */
358
359 nli = nli + (long int) nli_inc;
360 nps = nps + (long int) nps_inc;
361
362 if (printfl > 2)
363 KINPrintInfo(kin_mem, PRNT_NLI, "KINSPBCG", "KINSpbcgSolve", INFO_NLI, nli_inc);
364
365 if (ret != 0) ncfl++;
366 last_flag = ret;
367
368 if ( (ret != 0) && (ret != SPBCG_RES_REDUCED) ) {
369
370 /* Handle all failure returns from SpbcgSolve */
371
372 switch(ret) {
373 case SPBCG_PSOLVE_FAIL_REC:
374 case SPBCG_ATIMES_FAIL_REC:
375 return(1);
376 break;
377 case SPBCG_CONV_FAIL:
378 case SPBCG_MEM_NULL:
379 case SPBCG_ATIMES_FAIL_UNREC:
380 case SPBCG_PSOLVE_FAIL_UNREC:
381 return(-1);
382 break;
383 }
384 }
385
386 /* SpbcgSolve returned either SPBCG_SUCCESS or SPBCG_RES_REDUCED.
387
388 Compute the terms sJpnorm and sFdotJp for use in the linesearch
389 routine and in KINForcingTerm. Both of these terms are subsequently
390 corrected if the step is reduced by constraints or the linesearch.
391
392 sJpnorm is the norm of the scaled product (scaled by fscale) of the
393 current Jacobian matrix J and the step vector p (= solution vector xx).
394
395 sFdotJp is the dot product of the scaled f vector and the scaled
396 vector J*p, where the scaling uses fscale. */
397
398 ret = KINSpilsAtimes(kin_mem, xx, bb);
399 if (ret > 0) {
400 last_flag = SPBCG_ATIMES_FAIL_REC;
401 return(1);
402 }
403 else if (ret < 0) {
404 last_flag = SPBCG_ATIMES_FAIL_UNREC;
405 return(-1);
406 }
407
408 *sJpnorm = N_VWL2Norm(bb, fscale);
409 N_VProd(bb, fscale, bb);
410 N_VProd(bb, fscale, bb);
411 *sFdotJp = N_VDotProd(fval, bb);
412
413 if (printfl > 2) KINPrintInfo(kin_mem, PRNT_EPS, "KINSPBCG",
414 "KINSpbcgSolve", INFO_EPS, res_norm, eps);
415
416 return(0);
417 }
418
419 /*
420 * -----------------------------------------------------------------
421 * Function : KINSpbcgFree
422 * -----------------------------------------------------------------
423 * Frees memory specific to the SPBCG linear solver module.
424 * -----------------------------------------------------------------
425 */
426
KINSpbcgFree(KINMem kin_mem)427 static void KINSpbcgFree(KINMem kin_mem)
428 {
429 KINSpilsMem kinspils_mem;
430 SpbcgMem spbcg_mem;
431
432 kinspils_mem = (KINSpilsMem) lmem;
433
434 spbcg_mem = (SpbcgMem) spils_mem;
435 SpbcgFree(spbcg_mem);
436
437 if (kinspils_mem->s_pfree != NULL) (kinspils_mem->s_pfree)(kin_mem);
438
439 free(kinspils_mem); kinspils_mem = NULL;
440 }
441