Actual source code: snesngmres.c
petsc-3.3-p7 2013-05-11
1: #include <../src/snes/impls/ngmres/snesngmres.h> /*I "petscsnes.h" I*/
2: #include <petscblaslapack.h>
4: const char *SNESNGMRESRestartTypes[] = {"NONE","PERIODIC","DIFFERENCE","SNESNGMRESRestartType","SNES_NGMRES_RESTART_",0};
5: const char *SNESNGMRESSelectTypes[] = {"NONE","DIFFERENCE","LINESEARCH","SNESNGMRESSelectType","SNES_NGMRES_SELECT_",0};
9: PetscErrorCode SNESReset_NGMRES(SNES snes)
10: {
11: SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data;
15: VecDestroyVecs(ngmres->msize, &ngmres->Fdot);
16: VecDestroyVecs(ngmres->msize, &ngmres->Xdot);
17: SNESLineSearchDestroy(&ngmres->additive_linesearch);
19: return(0);
20: }
24: PetscErrorCode SNESDestroy_NGMRES(SNES snes)
25: {
27: SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
30: SNESReset_NGMRES(snes);
31: PetscFree5(ngmres->h, ngmres->beta, ngmres->xi, ngmres->fnorms, ngmres->q);
32: PetscFree(ngmres->s);
33: PetscFree(ngmres->xnorms);
34: #if PETSC_USE_COMPLEX
35: PetscFree(ngmres->rwork);
36: #endif
37: PetscFree(ngmres->work);
38: PetscFree(snes->data);
39: return(0);
40: }
44: PetscErrorCode SNESSetUp_NGMRES(SNES snes)
45: {
46: SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data;
47: const char *optionsprefix;
48: PetscInt msize,hsize;
52: SNESDefaultGetWork(snes,5);
53: if (!ngmres->Xdot) {VecDuplicateVecs(snes->vec_sol,ngmres->msize,&ngmres->Xdot);}
54: if (!ngmres->Fdot) {VecDuplicateVecs(snes->vec_sol,ngmres->msize,&ngmres->Fdot);}
55: if (!ngmres->setup_called) {
56: msize = ngmres->msize; /* restart size */
57: hsize = msize * msize;
59: /* explicit least squares minimization solve */
60: PetscMalloc5(hsize,PetscScalar,&ngmres->h,
61: msize,PetscScalar,&ngmres->beta,
62: msize,PetscScalar,&ngmres->xi,
63: msize,PetscReal, &ngmres->fnorms,
64: hsize,PetscScalar,&ngmres->q);
65: if (ngmres->singlereduction) {
66: PetscMalloc(msize*sizeof(PetscReal),&ngmres->xnorms);
67: }
68: ngmres->nrhs = 1;
69: ngmres->lda = msize;
70: ngmres->ldb = msize;
71: PetscMalloc(msize*sizeof(PetscScalar),&ngmres->s);
72: PetscMemzero(ngmres->h, hsize*sizeof(PetscScalar));
73: PetscMemzero(ngmres->q, hsize*sizeof(PetscScalar));
74: PetscMemzero(ngmres->xi, msize*sizeof(PetscScalar));
75: PetscMemzero(ngmres->beta,msize*sizeof(PetscScalar));
76: ngmres->lwork = 12*msize;
77: #if PETSC_USE_COMPLEX
78: PetscMalloc(sizeof(PetscReal)*ngmres->lwork,&ngmres->rwork);
79: #endif
80: PetscMalloc(sizeof(PetscScalar)*ngmres->lwork,&ngmres->work);
81: }
83: /* linesearch setup */
84: SNESGetOptionsPrefix(snes, &optionsprefix);
86: if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) {
87: SNESLineSearchCreate(((PetscObject)snes)->comm, &ngmres->additive_linesearch);
88: SNESLineSearchSetSNES(ngmres->additive_linesearch, snes);
89: SNESLineSearchSetType(ngmres->additive_linesearch, SNESLINESEARCHL2);
90: SNESLineSearchAppendOptionsPrefix(ngmres->additive_linesearch, "additive_");
91: SNESLineSearchAppendOptionsPrefix(ngmres->additive_linesearch, optionsprefix);
92: SNESLineSearchSetFromOptions(ngmres->additive_linesearch);
93: }
95: ngmres->setup_called = PETSC_TRUE;
96: return(0);
97: }
101: PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes)
102: {
103: SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data;
105: PetscBool debug;
106: SNESLineSearch linesearch;
108: PetscOptionsHead("SNES NGMRES options");
109: PetscOptionsEnum("-snes_ngmres_select_type","Select type","SNESNGMRESSetSelectType",SNESNGMRESSelectTypes,
110: (PetscEnum)ngmres->select_type,(PetscEnum*)&ngmres->select_type,PETSC_NULL);
111: PetscOptionsEnum("-snes_ngmres_restart_type","Restart type","SNESNGMRESSetRestartType",SNESNGMRESRestartTypes,
112: (PetscEnum)ngmres->restart_type,(PetscEnum*)&ngmres->restart_type,PETSC_NULL);
113: PetscOptionsBool("-snes_ngmres_anderson", "Use Anderson mixing storage", "SNES", ngmres->anderson, &ngmres->anderson, PETSC_NULL);
114: PetscOptionsInt("-snes_ngmres_m", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, PETSC_NULL);
115: PetscOptionsInt("-snes_ngmres_restart", "Iterations before forced restart", "SNES", ngmres->restart_periodic, &ngmres->restart_periodic, PETSC_NULL);
116: PetscOptionsInt("-snes_ngmres_restart_it","Tolerance iterations before restart","SNES", ngmres->restart_it, &ngmres->restart_it, PETSC_NULL);
117: PetscOptionsBool("-snes_ngmres_monitor", "Monitor actions of NGMRES", "SNES", ngmres->monitor ? PETSC_TRUE: PETSC_FALSE, &debug, PETSC_NULL);
118: if (debug) {
119: ngmres->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);
120: }
121: PetscOptionsReal("-snes_ngmres_gammaA", "Residual selection constant", "SNES", ngmres->gammaA, &ngmres->gammaA, PETSC_NULL);
122: PetscOptionsReal("-snes_ngmres_gammaC", " Residual restart constant", "SNES", ngmres->gammaC, &ngmres->gammaC, PETSC_NULL);
123: PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES", ngmres->epsilonB, &ngmres->epsilonB, PETSC_NULL);
124: PetscOptionsReal("-snes_ngmres_deltaB", "Difference residual selection constant", "SNES", ngmres->deltaB, &ngmres->deltaB, PETSC_NULL);
125: PetscOptionsBool("-snes_ngmres_single_reduction", "Aggregate reductions", "SNES", ngmres->singlereduction, &ngmres->singlereduction, PETSC_NULL);
126: PetscOptionsTail();
127: if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA;
129: /* set the default type of the line search if the user hasn't already. */
130: if (!snes->linesearch) {
131: SNESGetSNESLineSearch(snes, &linesearch);
132: SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);
133: }
134: return(0);
135: }
139: PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer)
140: {
141: SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data;
142: PetscBool iascii;
146: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
147: if (iascii) {
149: PetscViewerASCIIPrintf(viewer, " Number of stored past updates: %d\n", ngmres->msize);
150: PetscViewerASCIIPrintf(viewer, " Residual selection: gammaA=%1.0e, gammaC=%1.0e\n", ngmres->gammaA, ngmres->gammaC);
151: PetscViewerASCIIPrintf(viewer, " Difference restart: epsilonB=%1.0e, deltaB=%1.0e\n", ngmres->epsilonB, ngmres->deltaB);
152: }
153: return(0);
154: }
159: PetscErrorCode SNESSolve_NGMRES(SNES snes)
160: {
161: SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data;
162: /* present solution, residual, and preconditioned residual */
163: Vec X, F, B, D, Y;
165: /* candidate linear combination answers */
166: Vec XA, FA, XM, FM, FPC;
168: /* previous iterations to construct the subspace */
169: Vec *Fdot = ngmres->Fdot;
170: Vec *Xdot = ngmres->Xdot;
172: /* coefficients and RHS to the minimization problem */
173: PetscScalar *beta = ngmres->beta;
174: PetscScalar *xi = ngmres->xi;
175: PetscReal fnorm, fMnorm, fAnorm;
176: PetscReal nu;
177: PetscScalar alph_total = 0.;
178: PetscInt i, j, k, k_restart, l, ivec, restart_count = 0;
180: /* solution selection data */
181: PetscBool selectA, selectRestart;
182: PetscReal dnorm, dminnorm = 0.0, dcurnorm;
183: PetscReal fminnorm,xnorm,ynorm;
185: SNESConvergedReason reason;
186: PetscBool lssucceed,changed_y,changed_w;
187: PetscErrorCode ierr;
190: /* variable initialization */
191: snes->reason = SNES_CONVERGED_ITERATING;
192: X = snes->vec_sol;
193: F = snes->vec_func;
194: B = snes->vec_rhs;
195: XA = snes->vec_sol_update;
196: FA = snes->work[0];
197: D = snes->work[1];
199: /* work for the line search */
200: Y = snes->work[2];
201: XM = snes->work[3];
202: FM = snes->work[4];
204: PetscObjectTakeAccess(snes);
205: snes->iter = 0;
206: snes->norm = 0.;
207: PetscObjectGrantAccess(snes);
209: /* initialization */
211: /* r = F(x) */
212: if (!snes->vec_func_init_set) {
213: SNESComputeFunction(snes, X, F);
214: if (snes->domainerror) {
215: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
216: return(0);
217: }
218: } else {
219: snes->vec_func_init_set = PETSC_FALSE;
220: }
222: if (!snes->norm_init_set) {
223: VecNorm(F, NORM_2, &fnorm);
224: if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation");
225: } else {
226: fnorm = snes->norm_init;
227: snes->norm_init_set = PETSC_FALSE;
228: }
229: fminnorm = fnorm;
230: /* nu = (r, r) */
231: nu = fnorm*fnorm;
233: /* q_{00} = nu */
234: Q(0,0) = nu;
235: ngmres->fnorms[0] = fnorm;
236: /* Fdot[0] = F */
237: VecCopy(X, Xdot[0]);
238: VecCopy(F, Fdot[0]);
240: PetscObjectTakeAccess(snes);
241: snes->norm = fnorm;
242: PetscObjectGrantAccess(snes);
243: SNESLogConvHistory(snes, fnorm, 0);
244: SNESMonitor(snes, 0, fnorm);
245: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
246: if (snes->reason) return(0);
248: k_restart = 1;
249: l = 1;
250: for (k=1; k < snes->max_its+1; k++) {
251: /* select which vector of the stored subspace will be updated */
252: ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */
254: /* Computation of x^M */
255: if (snes->pc) {
256: VecCopy(X, XM);
257: SNESSetInitialFunction(snes->pc, F);
258: SNESSetInitialFunctionNorm(snes->pc, fnorm);
259: SNESSolve(snes->pc, B, XM);
260: SNESGetConvergedReason(snes->pc,&reason);
261: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
262: snes->reason = SNES_DIVERGED_INNER;
263: return(0);
264: }
265: SNESGetFunction(snes->pc, &FPC, PETSC_NULL, PETSC_NULL);
266: VecCopy(FPC, FM);
267: SNESGetFunctionNorm(snes->pc, &fMnorm);
268: } else {
269: /* no preconditioner -- just take gradient descent with line search */
270: VecCopy(F, Y);
271: VecCopy(F, FM);
272: VecCopy(X, XM);
273: fMnorm = fnorm;
274: SNESLineSearchApply(snes->linesearch,XM,FM,&fMnorm,Y);
275: SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);
276: if (!lssucceed) {
277: if (++snes->numFailures >= snes->maxFailures) {
278: snes->reason = SNES_DIVERGED_LINE_SEARCH;
279: return(0);
280: }
281: }
282: }
284: /* r = F(x) */
285: nu = fMnorm*fMnorm;
286: if (fminnorm > fMnorm) fminnorm = fMnorm; /* the minimum norm is now of F^M */
288: /* construct the right hand side and xi factors */
289: VecMDot(FM, l, Fdot, xi);
291: for (i = 0; i < l; i++) {
292: beta[i] = nu - xi[i];
293: }
295: /* construct h */
296: for (j = 0; j < l; j++) {
297: for (i = 0; i < l; i++) {
298: H(i, j) = Q(i, j) - xi[i] - xi[j] + nu;
299: }
300: }
302: if(l == 1) {
303: /* simply set alpha[0] = beta[0] / H[0, 0] */
304: if (H(0, 0) != 0.) {
305: beta[0] = beta[0] / H(0, 0);
306: } else {
307: beta[0] = 0.;
308: }
309: } else {
310: #ifdef PETSC_MISSING_LAPACK_GELSS
311: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine.");
312: #else
313: ngmres->m = PetscBLASIntCast(l);
314: ngmres->n = PetscBLASIntCast(l);
315: ngmres->info = PetscBLASIntCast(0);
316: ngmres->rcond = -1.;
317: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
318: #ifdef PETSC_USE_COMPLEX
319: LAPACKgelss_(&ngmres->m,
320: &ngmres->n,
321: &ngmres->nrhs,
322: ngmres->h,
323: &ngmres->lda,
324: ngmres->beta,
325: &ngmres->ldb,
326: ngmres->s,
327: &ngmres->rcond,
328: &ngmres->rank,
329: ngmres->work,
330: &ngmres->lwork,
331: ngmres->rwork,
332: &ngmres->info);
333: #else
334: LAPACKgelss_(&ngmres->m,
335: &ngmres->n,
336: &ngmres->nrhs,
337: ngmres->h,
338: &ngmres->lda,
339: ngmres->beta,
340: &ngmres->ldb,
341: ngmres->s,
342: &ngmres->rcond,
343: &ngmres->rank,
344: ngmres->work,
345: &ngmres->lwork,
346: &ngmres->info);
347: #endif
348: PetscFPTrapPop();
349: if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS");
350: if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge");
351: #endif
352: }
354: alph_total = 0.;
355: for (i = 0; i < l; i++) {
356: alph_total += beta[i];
357: }
359: VecCopy(XM, XA);
360: VecScale(XA, 1. - alph_total);
362: VecMAXPY(XA, l, beta, Xdot);
364: /* check the validity of the step */
365: VecCopy(XA,Y);
366: VecAXPY(Y,-1.0,X);
367: SNESLineSearchPostCheck(snes->linesearch,X,Y,XA,&changed_y,&changed_w);
368: SNESComputeFunction(snes, XA, FA);
370: /* differences for selection and restart */
371: if (ngmres->restart_type == SNES_NGMRES_RESTART_DIFFERENCE || ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) {
372: if (ngmres->singlereduction) {
373: dminnorm = -1.0;
374: ierr=VecCopy(XA, D);
375: ierr=VecAXPY(D, -1.0, XM);
376: for(i=0;i<l;i++) {
377: ierr=VecAXPY(Xdot[i],-1.0,XA);
378: }
379: VecNormBegin(FA, NORM_2, &fAnorm);
380: VecNormBegin(D, NORM_2, &dnorm);
381: for (i=0;i<l;i++) {
382: VecNormBegin(Xdot[i], NORM_2, &ngmres->xnorms[i]);
383: }
384: VecNormEnd(FA, NORM_2, &fAnorm);
385: VecNormEnd(D, NORM_2, &dnorm);
386: for (i=0;i<l;i++) {
387: VecNormEnd(Xdot[i], NORM_2, &ngmres->xnorms[i]);
388: }
389: for(i=0;i<l;i++) {
390: dcurnorm = ngmres->xnorms[i];
391: if ((dcurnorm < dminnorm) || (dminnorm < 0.0)) dminnorm = dcurnorm;
392: ierr=VecAXPY(Xdot[i],1.0,XA);
393: }
394: } else {
395: ierr=VecCopy(XA, D);
396: ierr=VecAXPY(D, -1.0, XM);
397: ierr=VecNormBegin(D, NORM_2, &dnorm);
398: ierr=VecNormBegin(FA, NORM_2, &fAnorm);
399: ierr=VecNormEnd(D, NORM_2, &dnorm);
400: ierr=VecNormEnd(FA, NORM_2, &fAnorm);
401: dminnorm = -1.0;
402: for(i=0;i<l;i++) {
403: ierr=VecCopy(XA, D);
404: ierr=VecAXPY(D, -1.0, Xdot[i]);
405: ierr=VecNorm(D, NORM_2, &dcurnorm);
406: if((dcurnorm < dminnorm) || (dminnorm < 0.0)) dminnorm = dcurnorm;
407: }
408: }
409: } else {
410: VecNorm(FA, NORM_2, &fAnorm);
411: }
412: /* combination (additive) or selection (multiplicative) of the N-GMRES solution */
413: if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) {
414: /* X = X + \lambda(XA - X) */
415: if (ngmres->monitor) {
416: PetscViewerASCIIPrintf(ngmres->monitor, "||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fMnorm);
417: }
418: VecCopy(FM, F);
419: VecCopy(XM, X);
420: VecCopy(XA, Y);
421: VecAYPX(Y, -1.0, X);
422: fnorm = fMnorm;
423: SNESLineSearchApply(ngmres->additive_linesearch,X,F,&fnorm,Y);
424: SNESLineSearchGetSuccess(ngmres->additive_linesearch, &lssucceed);
425: if (!lssucceed) {
426: if (++snes->numFailures >= snes->maxFailures) {
427: snes->reason = SNES_DIVERGED_LINE_SEARCH;
428: return(0);
429: }
430: }
431: if (ngmres->monitor) {
432: PetscViewerASCIIPrintf(ngmres->monitor, "Additive solution: ||F||_2 = %e\n", fnorm);
433: }
434: } else if (ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) {
435: selectA = PETSC_TRUE;
436: /* Conditions for choosing the accelerated answer */
437: /* Criterion A -- the norm of the function isn't increased above the minimum by too much */
438: if (fAnorm >= ngmres->gammaA*fminnorm) {
439: selectA = PETSC_FALSE;
440: }
441: /* Criterion B -- the choice of x^A isn't too close to some other choice */
442: if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) {
443: } else {
444: selectA=PETSC_FALSE;
445: }
446: if (selectA) {
447: if (ngmres->monitor) {
448: PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fMnorm);
449: }
450: /* copy it over */
451: fnorm = fAnorm;
452: nu = fnorm*fnorm;
453: VecCopy(FA, F);
454: VecCopy(XA, X);
455: } else {
456: if (ngmres->monitor) {
457: PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fMnorm);
458: }
459: fnorm = fMnorm;
460: nu = fnorm*fnorm;
461: VecCopy(XM, Y);
462: VecAXPY(Y,-1.0,X);
463: VecCopy(FM, F);
464: VecCopy(XM, X);
465: }
466: } else { /* none */
467: fnorm = fAnorm;
468: nu = fnorm*fnorm;
469: VecCopy(FA, F);
470: VecCopy(XA, X);
471: }
473: selectRestart = PETSC_FALSE;
474: if (ngmres->restart_type == SNES_NGMRES_RESTART_DIFFERENCE) {
475: /* difference stagnation restart */
476: if((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm))) {
477: if (ngmres->monitor) {
478: PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*dnorm, dminnorm);
479: }
480: selectRestart = PETSC_TRUE;
481: }
482: /* residual stagnation restart */
483: if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) {
484: if (ngmres->monitor) {
485: PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", PetscSqrtReal(fAnorm), ngmres->gammaC*PetscSqrtReal(fminnorm));
486: }
487: selectRestart = PETSC_TRUE;
488: }
489: /* if the restart conditions persist for more than restart_it iterations, restart. */
490: if (selectRestart) {
491: restart_count++;
492: } else {
493: restart_count = 0;
494: }
495: } else if (ngmres->restart_type == SNES_NGMRES_RESTART_PERIODIC) {
496: if (k_restart > ngmres->restart_periodic) {
497: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "periodic restart after %D iterations\n", k_restart);
498: restart_count = ngmres->restart_it;
499: }
500: }
501: /* restart after restart conditions have persisted for a fixed number of iterations */
502: if (restart_count >= ngmres->restart_it) {
503: if (ngmres->monitor){
504: PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);
505: }
506: restart_count = 0;
507: k_restart = 1;
508: l = 1;
509: /* q_{00} = nu */
510: if (ngmres->anderson) {
511: ngmres->fnorms[0] = fMnorm;
512: nu = fMnorm*fMnorm;
513: Q(0,0) = nu;
514: /* Fdot[0] = F */
515: VecCopy(XM, Xdot[0]);
516: VecCopy(FM, Fdot[0]);
517: } else {
518: ngmres->fnorms[0] = fnorm;
519: nu = fnorm*fnorm;
520: Q(0,0) = nu;
521: /* Fdot[0] = F */
522: VecCopy(X, Xdot[0]);
523: VecCopy(F, Fdot[0]);
524: }
525: } else {
526: /* select the current size of the subspace */
527: if (l < ngmres->msize) l++;
528: k_restart++;
529: /* place the current entry in the list of previous entries */
530: if (ngmres->anderson) {
531: VecCopy(FM, Fdot[ivec]);
532: VecCopy(XM, Xdot[ivec]);
533: ngmres->fnorms[ivec] = fMnorm;
534: if (fminnorm > fMnorm) fminnorm = fMnorm; /* the minimum norm is now of FM */
535: xi[ivec] = fMnorm*fMnorm;
536: for (i = 0; i < l; i++) {
537: Q(i, ivec) = xi[i];
538: Q(ivec, i) = xi[i];
539: }
540: } else {
541: VecCopy(F, Fdot[ivec]);
542: VecCopy(X, Xdot[ivec]);
543: ngmres->fnorms[ivec] = fnorm;
544: if (fminnorm > fnorm) fminnorm = fnorm; /* the minimum norm is now of FA */
545: VecMDot(F, l, Fdot, xi);
546: for (i = 0; i < l; i++) {
547: Q(i, ivec) = xi[i];
548: Q(ivec, i) = xi[i];
549: }
550: }
551: }
553: PetscObjectTakeAccess(snes);
554: snes->iter = k;
555: snes->norm = fnorm;
556: PetscObjectGrantAccess(snes);
557: SNESLogConvHistory(snes, snes->norm, snes->iter);
558: SNESMonitor(snes, snes->iter, snes->norm);
559: VecNormBegin(Y,NORM_2,&ynorm);
560: VecNormBegin(X,NORM_2,&xnorm);
561: VecNormEnd(Y,NORM_2,&ynorm);
562: VecNormEnd(X,NORM_2,&xnorm);
563: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
564: if (snes->reason) return(0);
565: }
566: snes->reason = SNES_DIVERGED_MAX_IT;
567: return(0);
568: }
572: /*@
573: SNESNGMRESSetRestartType - Sets the restart type for SNESNGMRES.
575: Logically Collective on SNES
577: Input Parameters:
578: + snes - the iterative context
579: - rtype - restart type
581: Options Database:
582: + -snes_ngmres_restart_type<difference,periodic,none> - set the restart type
583: - -snes_ngmres_restart[30] - sets the number of iterations before restart for periodic
585: Level: intermediate
587: SNESNGMRESRestartTypes:
588: + SNES_NGMRES_RESTART_NONE - never restart
589: . SNES_NGMRES_RESTART_DIFFERENCE - restart based upon difference criteria
590: - SNES_NGMRES_RESTART_PERIODIC - restart after a fixed number of iterations
592: Notes:
593: The default line search used is the L2 line search and it requires two additional function evaluations.
595: .keywords: SNES, SNESNGMRES, restart, type, set SNESLineSearch
596: @*/
597: PetscErrorCode SNESNGMRESSetRestartType(SNES snes, SNESNGMRESRestartType rtype) {
601: PetscTryMethod(snes,"SNESNGMRESSetRestartType_C",(SNES,SNESNGMRESRestartType),(snes,rtype));
602: return(0);
603: }
607: /*@
608: SNESNGMRESSetSelectType - Sets the selection type for SNESNGMRES. This determines how the candidate solution and
609: combined solution are used to create the next iterate.
611: Logically Collective on SNES
613: Input Parameters:
614: + snes - the iterative context
615: - stype - selection type
617: Options Database:
618: . -snes_ngmres_select_type<difference,none,linesearch>
620: Level: intermediate
622: SNESNGMRESSelectTypes:
623: + SNES_NGMRES_SELECT_NONE - choose the combined solution all the time
624: . SNES_NGMRES_SELECT_DIFFERENCE - choose based upon the selection criteria
625: - SNES_NGMRES_SELECT_LINESEARCH - choose based upon line search combination
627: Notes:
628: The default line search used is the L2 line search and it requires two additional function evaluations.
630: .keywords: SNES, SNESNGMRES, selection, type, set SNESLineSearch
631: @*/
633: PetscErrorCode SNESNGMRESSetSelectType(SNES snes, SNESNGMRESSelectType stype) {
637: PetscTryMethod(snes,"SNESNGMRESSetSelectType_C",(SNES,SNESNGMRESSelectType),(snes,stype));
638: return(0);
639: }
642: EXTERN_C_BEGIN
646: PetscErrorCode SNESNGMRESSetSelectType_NGMRES(SNES snes, SNESNGMRESSelectType stype) {
647: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
649: ngmres->select_type = stype;
650: return(0);
651: }
656: PetscErrorCode SNESNGMRESSetRestartType_NGMRES(SNES snes, SNESNGMRESRestartType rtype) {
657: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
659: ngmres->restart_type = rtype;
660: return(0);
661: }
662: EXTERN_C_END
665: /*MC
666: SNESNGMRES - The Nonlinear Generalized Minimum Residual method.
668: Level: beginner
670: Options Database:
671: + -snes_ngmres_select_type<difference,none,linesearch> - choose the select between candidate and combined solution
672: + -snes_ngmres_restart_type<difference,none,periodic> - choose the restart conditions
673: . -snes_ngmres_anderson - Use Anderson mixing NGMRES variant which combines candidate solutions instead of actual solutions
674: . -snes_ngmres_m - Number of stored previous solutions and residuals
675: . -snes_ngmres_restart_it - Number of iterations the restart conditions hold before restart
676: . -snes_ngmres_gammaA - Residual tolerance for solution select between the candidate and combination
677: . -snes_ngmres_gammaC - Residual tolerance for restart
678: . -snes_ngmres_epsilonB - Difference tolerance between subsequent solutions triggering restart
679: . -snes_ngmres_deltaB - Difference tolerance between residuals triggering restart
680: . -snes_ngmres_monitor - Prints relevant information about the ngmres iteration
681: . -snes_linesearch_type <basic,basicnonorms,quadratic,critical> - Line search type used for the default smoother
682: - -additive_snes_linesearch_type - linesearch type used to select between the candidate and combined solution with additive select type
684: Notes:
686: The N-GMRES method combines m previous solutions into a minimum-residual solution by solving a small linearized
687: optimization problem at each iteration.
689: References:
691: "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio,
692: SIAM Journal on Scientific Computing, 21(5), 2000.
694: This is also the same as the algorithm called Anderson acceleration i