1 /*
2 * -----------------------------------------------------------------
3 * $Revision: 4272 $
4 * $Date: 2014-12-02 11:19:41 -0800 (Tue, 02 Dec 2014) $
5 * -----------------------------------------------------------------
6 * Programmer(s): 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 CVSPGMR linear solver.
19 * -----------------------------------------------------------------
20 */
21
22 #include <stdio.h>
23 #include <stdlib.h>
24
25 #include <cvodes/cvodes_spgmr.h>
26 #include "cvodes_spils_impl.h"
27 #include "cvodes_impl.h"
28
29 #include <sundials/sundials_spgmr.h>
30 #include <sundials/sundials_math.h>
31
32 /* Constants */
33
34 #define ZERO RCONST(0.0)
35 #define ONE RCONST(1.0)
36
37 /* CVSPGMR linit, lsetup, lsolve, and lfree routines */
38
39 static int CVSpgmrInit(CVodeMem cv_mem);
40
41 static int CVSpgmrSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
42 N_Vector fpred, booleantype *jcurPtr, N_Vector vtemp1,
43 N_Vector vtemp2, N_Vector vtemp3);
44
45 static int CVSpgmrSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
46 N_Vector ynow, N_Vector fnow);
47
48 static void CVSpgmrFree(CVodeMem cv_mem);
49
50 /* CVSPGMR lfreeB function */
51
52 static void CVSpgmrFreeB(CVodeBMem cvB_mem);
53
54 /*
55 * ================================================================
56 *
57 * PART I - forward problems
58 *
59 * ================================================================
60 */
61
62
63 /* Readability Replacements */
64
65 #define tq (cv_mem->cv_tq)
66 #define nst (cv_mem->cv_nst)
67 #define tn (cv_mem->cv_tn)
68 #define h (cv_mem->cv_h)
69 #define gamma (cv_mem->cv_gamma)
70 #define gammap (cv_mem->cv_gammap)
71 #define f (cv_mem->cv_f)
72 #define user_data (cv_mem->cv_user_data)
73 #define ewt (cv_mem->cv_ewt)
74 #define mnewt (cv_mem->cv_mnewt)
75 #define ropt (cv_mem->cv_ropt)
76 #define linit (cv_mem->cv_linit)
77 #define lsetup (cv_mem->cv_lsetup)
78 #define lsolve (cv_mem->cv_lsolve)
79 #define lfree (cv_mem->cv_lfree)
80 #define lmem (cv_mem->cv_lmem)
81 #define vec_tmpl (cv_mem->cv_tempv)
82 #define setupNonNull (cv_mem->cv_setupNonNull)
83
84 #define sqrtN (cvspils_mem->s_sqrtN)
85 #define ytemp (cvspils_mem->s_ytemp)
86 #define x (cvspils_mem->s_x)
87 #define ycur (cvspils_mem->s_ycur)
88 #define fcur (cvspils_mem->s_fcur)
89 #define delta (cvspils_mem->s_delta)
90 #define deltar (cvspils_mem->s_deltar)
91 #define npe (cvspils_mem->s_npe)
92 #define nli (cvspils_mem->s_nli)
93 #define nps (cvspils_mem->s_nps)
94 #define ncfl (cvspils_mem->s_ncfl)
95 #define nstlpre (cvspils_mem->s_nstlpre)
96 #define njtimes (cvspils_mem->s_njtimes)
97 #define nfes (cvspils_mem->s_nfes)
98 #define spils_mem (cvspils_mem->s_spils_mem)
99
100 #define jtimesDQ (cvspils_mem->s_jtimesDQ)
101 #define jtimes (cvspils_mem->s_jtimes)
102 #define j_data (cvspils_mem->s_j_data)
103
104 #define last_flag (cvspils_mem->s_last_flag)
105
106 /*
107 * -----------------------------------------------------------------
108 * CVSpgmr
109 * -----------------------------------------------------------------
110 * This routine initializes the memory record and sets various function
111 * fields specific to the Spgmr linear solver module. CVSpgmr first
112 * calls the existing lfree routine if this is not NULL. It then sets
113 * the cv_linit, cv_lsetup, cv_lsolve, cv_lfree fields in (*cvode_mem)
114 * to be CVSpgmrInit, CVSpgmrSetup, CVSpgmrSolve, and CVSpgmrFree,
115 * respectively. It allocates memory for a structure of type
116 * CVSpilsMemRec and sets the cv_lmem field in (*cvode_mem) to the
117 * address of this structure. It sets setupNonNull in (*cvode_mem),
118 * and sets various fields in the CVSpilsMemRec structure.
119 * Finally, CVSpgmr allocates memory for ytemp and x, and calls
120 * SpgmrMalloc to allocate memory for the Spgmr solver.
121 * -----------------------------------------------------------------
122 */
123
CVSpgmr(void * cvode_mem,int pretype,int maxl)124 int CVSpgmr(void *cvode_mem, int pretype, int maxl)
125 {
126 CVodeMem cv_mem;
127 CVSpilsMem cvspils_mem;
128 SpgmrMem spgmr_mem;
129 int mxl;
130
131 /* Return immediately if cvode_mem is NULL */
132 if (cvode_mem == NULL) {
133 cvProcessError(NULL, CVSPILS_MEM_NULL, "CVSPGMR", "CVSpgmr", MSGS_CVMEM_NULL);
134 return(CVSPILS_MEM_NULL);
135 }
136 cv_mem = (CVodeMem) cvode_mem;
137
138 /* Check if N_VDotProd is present */
139 if(vec_tmpl->ops->nvdotprod == NULL) {
140 cvProcessError(cv_mem, CVSPILS_ILL_INPUT, "CVSPGMR", "CVSpgmr", MSGS_BAD_NVECTOR);
141 return(CVSPILS_ILL_INPUT);
142 }
143
144 if (lfree != NULL) lfree(cv_mem);
145
146 /* Set four main function fields in cv_mem */
147 linit = CVSpgmrInit;
148 lsetup = CVSpgmrSetup;
149 lsolve = CVSpgmrSolve;
150 lfree = CVSpgmrFree;
151
152 /* Get memory for CVSpilsMemRec */
153 cvspils_mem = NULL;
154 cvspils_mem = (CVSpilsMem) malloc(sizeof(struct CVSpilsMemRec));
155 if (cvspils_mem == NULL) {
156 cvProcessError(cv_mem, CVSPILS_MEM_FAIL, "CVSPGMR", "CVSpgmr", MSGS_MEM_FAIL);
157 return(CVSPILS_MEM_FAIL);
158 }
159
160 /* Set ILS type */
161 cvspils_mem->s_type = SPILS_SPGMR;
162
163 /* Set Spgmr parameters that have been passed in call sequence */
164 cvspils_mem->s_pretype = pretype;
165 mxl = cvspils_mem->s_maxl = (maxl <= 0) ? CVSPILS_MAXL : maxl;
166
167 /* Set defaults for Jacobian-related fileds */
168 jtimesDQ = TRUE;
169 jtimes = NULL;
170 j_data = NULL;
171
172 /* Set defaults for preconditioner-related fields */
173 cvspils_mem->s_pset = NULL;
174 cvspils_mem->s_psolve = NULL;
175 cvspils_mem->s_pfree = NULL;
176 cvspils_mem->s_P_data = cv_mem->cv_user_data;
177
178 /* Set default values for the rest of the Spgmr parameters */
179 cvspils_mem->s_gstype = MODIFIED_GS;
180 cvspils_mem->s_eplifac = CVSPILS_EPLIN;
181
182 cvspils_mem->s_last_flag = CVSPILS_SUCCESS;
183
184 setupNonNull = FALSE;
185
186 /* Check for legal pretype */
187 if ((pretype != PREC_NONE) && (pretype != PREC_LEFT) &&
188 (pretype != PREC_RIGHT) && (pretype != PREC_BOTH)) {
189 cvProcessError(cv_mem, CVSPILS_ILL_INPUT, "CVSPGMR", "CVSpgmr", MSGS_BAD_PRETYPE);
190 free(cvspils_mem); cvspils_mem = NULL;
191 return(CVSPILS_ILL_INPUT);
192 }
193
194 /* Allocate memory for ytemp and x */
195 ytemp = N_VClone(vec_tmpl);
196 if (ytemp == NULL) {
197 cvProcessError(cv_mem, CVSPILS_MEM_FAIL, "CVSPGMR", "CVSpgmr", MSGS_MEM_FAIL);
198 free(cvspils_mem); cvspils_mem = NULL;
199 return(CVSPILS_MEM_FAIL);
200 }
201 x = N_VClone(vec_tmpl);
202 if (x == NULL) {
203 cvProcessError(cv_mem, CVSPILS_MEM_FAIL, "CVSPGMR", "CVSpgmr", MSGS_MEM_FAIL);
204 N_VDestroy(ytemp);
205 free(cvspils_mem); cvspils_mem = NULL;
206 return(CVSPILS_MEM_FAIL);
207 }
208
209 /* Compute sqrtN from a dot product */
210 N_VConst(ONE, ytemp);
211 sqrtN = SUNRsqrt( N_VDotProd(ytemp, ytemp) );
212
213 /* Call SpgmrMalloc to allocate workspace for Spgmr */
214 spgmr_mem = NULL;
215 spgmr_mem = SpgmrMalloc(mxl, vec_tmpl);
216 if (spgmr_mem == NULL) {
217 cvProcessError(cv_mem, CVSPILS_MEM_FAIL, "CVSPGMR", "CVSpgmr", MSGS_MEM_FAIL);
218 N_VDestroy(ytemp);
219 N_VDestroy(x);
220 free(cvspils_mem); cvspils_mem = NULL;
221 return(CVSPILS_MEM_FAIL);
222 }
223
224 /* Attach SPGMR memory to spils memory structure */
225 spils_mem = (void *) spgmr_mem;
226
227 /* Attach linear solver memory to integrator memory */
228 lmem = cvspils_mem;
229
230 return(CVSPILS_SUCCESS);
231 }
232
233
234 /* Additional readability Replacements */
235
236 #define pretype (cvspils_mem->s_pretype)
237 #define gstype (cvspils_mem->s_gstype)
238 #define eplifac (cvspils_mem->s_eplifac)
239 #define maxl (cvspils_mem->s_maxl)
240 #define psolve (cvspils_mem->s_psolve)
241 #define pset (cvspils_mem->s_pset)
242 #define P_data (cvspils_mem->s_P_data)
243
244 /*
245 * -----------------------------------------------------------------
246 * CVSpgmrInit
247 * -----------------------------------------------------------------
248 * This routine does remaining initializations specific to the Spgmr
249 * linear solver.
250 * -----------------------------------------------------------------
251 */
252
CVSpgmrInit(CVodeMem cv_mem)253 static int CVSpgmrInit(CVodeMem cv_mem)
254 {
255 CVSpilsMem cvspils_mem;
256 cvspils_mem = (CVSpilsMem) lmem;
257
258 /* Initialize counters */
259 npe = nli = nps = ncfl = nstlpre = 0;
260 njtimes = nfes = 0;
261
262 /* Check for legal combination pretype - psolve */
263 if ((pretype != PREC_NONE) && (psolve == NULL)) {
264 cvProcessError(cv_mem, -1, "CVSPGMR", "CVSpgmrInit", MSGS_PSOLVE_REQ);
265 last_flag = CVSPILS_ILL_INPUT;
266 return(-1);
267 }
268
269 /* Set setupNonNull = TRUE iff there is preconditioning (pretype != PREC_NONE)
270 and there is a preconditioning setup phase (pset != NULL) */
271 setupNonNull = (pretype != PREC_NONE) && (pset != NULL);
272
273 /* Set Jacobian-related fields, based on jtimesDQ */
274 if (jtimesDQ) {
275 jtimes = CVSpilsDQJtimes;
276 j_data = cv_mem;
277 } else {
278 j_data = user_data;
279 }
280
281 last_flag = CVSPILS_SUCCESS;
282 return(0);
283 }
284
285 /*
286 * -----------------------------------------------------------------
287 * CVSpgmrSetup
288 * -----------------------------------------------------------------
289 * This routine does the setup operations for the Spgmr linear solver.
290 * It makes a decision as to whether or not to signal for re-evaluation
291 * of Jacobian data in the pset routine, based on various state
292 * variables, then it calls pset. If we signal for re-evaluation,
293 * then we reset jcur = *jcurPtr to TRUE, regardless of the pset output.
294 * In any case, if jcur == TRUE, we increment npe and save nst in nstlpre.
295 * -----------------------------------------------------------------
296 */
297
CVSpgmrSetup(CVodeMem cv_mem,int convfail,N_Vector ypred,N_Vector fpred,booleantype * jcurPtr,N_Vector vtemp1,N_Vector vtemp2,N_Vector vtemp3)298 static int CVSpgmrSetup(CVodeMem cv_mem, int convfail, N_Vector ypred,
299 N_Vector fpred, booleantype *jcurPtr, N_Vector vtemp1,
300 N_Vector vtemp2, N_Vector vtemp3)
301 {
302 booleantype jbad, jok;
303 realtype dgamma;
304 int retval;
305 CVSpilsMem cvspils_mem;
306
307 cvspils_mem = (CVSpilsMem) lmem;
308
309 /* Use nst, gamma/gammap, and convfail to set J eval. flag jok */
310 dgamma = SUNRabs((gamma/gammap) - ONE);
311 jbad = (nst == 0) || (nst > nstlpre + CVSPILS_MSBPRE) ||
312 ((convfail == CV_FAIL_BAD_J) && (dgamma < CVSPILS_DGMAX)) ||
313 (convfail == CV_FAIL_OTHER);
314 *jcurPtr = jbad;
315 jok = !jbad;
316
317 /* Call pset routine and possibly reset jcur */
318 retval = pset(tn, ypred, fpred, jok, jcurPtr, gamma, P_data,
319 vtemp1, vtemp2, vtemp3);
320 if (retval < 0) {
321 cvProcessError(cv_mem, SPGMR_PSET_FAIL_UNREC, "CVSPGMR", "CVSpgmrSetup", MSGS_PSET_FAILED);
322 last_flag = SPGMR_PSET_FAIL_UNREC;
323 }
324 if (retval > 0) {
325 last_flag = SPGMR_PSET_FAIL_REC;
326 }
327
328 if (jbad) *jcurPtr = TRUE;
329
330 /* If jcur = TRUE, increment npe and save nst value */
331 if (*jcurPtr) {
332 npe++;
333 nstlpre = nst;
334 }
335
336 last_flag = SPGMR_SUCCESS;
337
338 /* Return the same value that pset returned */
339 return(retval);
340 }
341
342 /*
343 * -----------------------------------------------------------------
344 * CVSpgmrSolve
345 * -----------------------------------------------------------------
346 * This routine handles the call to the generic solver SpgmrSolve
347 * for the solution of the linear system Ax = b with the SPGMR method,
348 * without restarts. The solution x is returned in the vector b.
349 *
350 * If the WRMS norm of b is small, we return x = b (if this is the first
351 * Newton iteration) or x = 0 (if a later Newton iteration).
352 *
353 * Otherwise, we set the tolerance parameter and initial guess (x = 0),
354 * call SpgmrSolve, and copy the solution x into b. The x-scaling and
355 * b-scaling arrays are both equal to weight, and no restarts are allowed.
356 *
357 * The counters nli, nps, and ncfl are incremented, and the return value
358 * is set according to the success of SpgmrSolve. The success flag is
359 * returned if SpgmrSolve converged, or if this is the first Newton
360 * iteration and the residual norm was reduced below its initial value.
361 * -----------------------------------------------------------------
362 */
363
CVSpgmrSolve(CVodeMem cv_mem,N_Vector b,N_Vector weight,N_Vector ynow,N_Vector fnow)364 static int CVSpgmrSolve(CVodeMem cv_mem, N_Vector b, N_Vector weight,
365 N_Vector ynow, N_Vector fnow)
366 {
367 realtype bnorm, res_norm;
368 CVSpilsMem cvspils_mem;
369 SpgmrMem spgmr_mem;
370 int nli_inc, nps_inc, retval;
371
372 cvspils_mem = (CVSpilsMem) lmem;
373
374 spgmr_mem = (SpgmrMem) spils_mem;
375
376 /* Test norm(b); if small, return x = 0 or x = b */
377 deltar = eplifac*tq[4];
378
379 bnorm = N_VWrmsNorm(b, weight);
380 if (bnorm <= deltar) {
381 if (mnewt > 0) N_VConst(ZERO, b);
382 return(0);
383 }
384
385 /* Set vectors ycur and fcur for use by the Atimes and Psolve routines */
386 ycur = ynow;
387 fcur = fnow;
388
389 /* Set inputs delta and initial guess x = 0 to SpgmrSolve */
390 delta = deltar * sqrtN;
391 N_VConst(ZERO, x);
392
393 /* Call SpgmrSolve and copy x to b */
394 retval = SpgmrSolve(spgmr_mem, cv_mem, x, b, pretype, gstype, delta, 0,
395 cv_mem, weight, weight, CVSpilsAtimes, CVSpilsPSolve,
396 &res_norm, &nli_inc, &nps_inc);
397
398 N_VScale(ONE, x, b);
399
400 /* Increment counters nli, nps, and ncfl */
401 nli += nli_inc;
402 nps += nps_inc;
403 if (retval != SPGMR_SUCCESS) ncfl++;
404
405 /* Interpret return value from SpgmrSolve */
406
407 last_flag = retval;
408
409 switch(retval) {
410
411 case SPGMR_SUCCESS:
412 return(0);
413 break;
414 case SPGMR_RES_REDUCED:
415 if (mnewt == 0) return(0);
416 else return(1);
417 break;
418 case SPGMR_CONV_FAIL:
419 return(1);
420 break;
421 case SPGMR_QRFACT_FAIL:
422 return(1);
423 break;
424 case SPGMR_PSOLVE_FAIL_REC:
425 return(1);
426 break;
427 case SPGMR_ATIMES_FAIL_REC:
428 return(1);
429 break;
430 case SPGMR_MEM_NULL:
431 return(-1);
432 break;
433 case SPGMR_ATIMES_FAIL_UNREC:
434 cvProcessError(cv_mem, SPGMR_ATIMES_FAIL_UNREC, "CVSPGMR", "CVSpgmrSolve", MSGS_JTIMES_FAILED);
435 return(-1);
436 break;
437 case SPGMR_PSOLVE_FAIL_UNREC:
438 cvProcessError(cv_mem, SPGMR_PSOLVE_FAIL_UNREC, "CVSPGMR", "CVSpgmrSolve", MSGS_PSOLVE_FAILED);
439 return(-1);
440 break;
441 case SPGMR_GS_FAIL:
442 return(-1);
443 break;
444 case SPGMR_QRSOL_FAIL:
445 return(-1);
446 break;
447 }
448
449 return(0);
450
451 }
452
453 /*
454 * -----------------------------------------------------------------
455 * CVSpgmrFree
456 * -----------------------------------------------------------------
457 * This routine frees memory specific to the Spgmr linear solver.
458 * -----------------------------------------------------------------
459 */
460
CVSpgmrFree(CVodeMem cv_mem)461 static void CVSpgmrFree(CVodeMem cv_mem)
462 {
463 CVSpilsMem cvspils_mem;
464 SpgmrMem spgmr_mem;
465
466 cvspils_mem = (CVSpilsMem) lmem;
467
468 N_VDestroy(ytemp);
469 N_VDestroy(x);
470
471 spgmr_mem = (SpgmrMem) spils_mem;
472 SpgmrFree(spgmr_mem);
473
474 if (cvspils_mem->s_pfree != NULL) (cvspils_mem->s_pfree)(cv_mem);
475
476 free(cvspils_mem);
477 cv_mem->cv_lmem = NULL;
478 }
479
480
481 /*
482 * ================================================================
483 *
484 * PART II - backward problems
485 *
486 * ================================================================
487 */
488
489
490 /* Additional readability replacements */
491
492 #define pset_B (cvspilsB_mem->s_psetB)
493 #define psolve_B (cvspilsB_mem->s_psolveB)
494 #define jtimes_B (cvspilsB_mem->s_jtimesB)
495 #define P_data_B (cvspilsB_mem->s_P_dataB)
496
497 /*
498 * CVSpgmrB
499 *
500 * Wrapper for the backward phase
501 *
502 */
503
CVSpgmrB(void * cvode_mem,int which,int pretypeB,int maxlB)504 int CVSpgmrB(void *cvode_mem, int which, int pretypeB, int maxlB)
505 {
506 CVodeMem cv_mem;
507 CVadjMem ca_mem;
508 CVodeBMem cvB_mem;
509 void *cvodeB_mem;
510 CVSpilsMemB cvspilsB_mem;
511 int flag;
512
513 /* Check if cvode_mem exists */
514 if (cvode_mem == NULL) {
515 cvProcessError(NULL, CVSPILS_MEM_NULL, "CVSPGMR", "CVSpgmrB", MSGS_CVMEM_NULL);
516 return(CVSPILS_MEM_NULL);
517 }
518 cv_mem = (CVodeMem) cvode_mem;
519
520 /* Was ASA initialized? */
521 if (cv_mem->cv_adjMallocDone == FALSE) {
522 cvProcessError(cv_mem, CVSPILS_NO_ADJ, "CVSPGMR", "CVSpgmrB", MSGS_NO_ADJ);
523 return(CVSPILS_NO_ADJ);
524 }
525 ca_mem = cv_mem->cv_adj_mem;
526
527 /* Check which */
528 if ( which >= ca_mem->ca_nbckpbs ) {
529 cvProcessError(cv_mem, CVSPILS_ILL_INPUT, "CVSPGMR", "CVSpgmrB", MSGS_BAD_WHICH);
530 return(CVSPILS_ILL_INPUT);
531 }
532
533 /* Find the CVodeBMem entry in the linked list corresponding to which */
534 cvB_mem = ca_mem->cvB_mem;
535 while (cvB_mem != NULL) {
536 if ( which == cvB_mem->cv_index ) break;
537 cvB_mem = cvB_mem->cv_next;
538 }
539
540 cvodeB_mem = (void *) (cvB_mem->cv_mem);
541
542 /* Get memory for CVSpilsMemRecB */
543 cvspilsB_mem = NULL;
544 cvspilsB_mem = (CVSpilsMemB) malloc(sizeof(struct CVSpilsMemRecB));
545 if (cvspilsB_mem == NULL) {
546 cvProcessError(cv_mem, CVSPILS_MEM_FAIL, "CVSPGMR", "CVSpgmrB", MSGS_MEM_FAIL);
547 return(CVSPILS_MEM_FAIL);
548 }
549
550 pset_B = NULL;
551 psolve_B = NULL;
552 P_data_B = NULL;
553
554 /* initialize Jacobian function */
555 jtimes_B = NULL;
556
557 /* attach lmemB and lfreeB */
558 cvB_mem->cv_lmem = cvspilsB_mem;
559 cvB_mem->cv_lfree = CVSpgmrFreeB;
560
561 flag = CVSpgmr(cvodeB_mem, pretypeB, maxlB);
562
563 if (flag != CVSPILS_SUCCESS) {
564 free(cvspilsB_mem);
565 cvspilsB_mem = NULL;
566 }
567
568 return(flag);
569 }
570
571
572 /*
573 * CVSpgmrFreeB
574 */
575
576
CVSpgmrFreeB(CVodeBMem cvB_mem)577 static void CVSpgmrFreeB(CVodeBMem cvB_mem)
578 {
579 CVSpilsMemB cvspilsB_mem;
580
581 cvspilsB_mem = (CVSpilsMemB) (cvB_mem->cv_lmem);
582
583 free(cvspilsB_mem);
584 }
585