Actual source code: snesngmres.c
petsc-3.4.5 2014-06-29
1: #include <../src/snes/impls/ngmres/snesngmres.h> /*I "petscsnes.h" I*/
2: #include <petscblaslapack.h>
4: const char *const SNESNGMRESRestartTypes[] = {"NONE","PERIODIC","DIFFERENCE","SNESNGMRESRestartType","SNES_NGMRES_RESTART_",0};
5: const char *const 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);
18: return(0);
19: }
23: PetscErrorCode SNESDestroy_NGMRES(SNES snes)
24: {
26: SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
29: SNESReset_NGMRES(snes);
30: PetscFree5(ngmres->h,ngmres->beta,ngmres->xi,ngmres->fnorms,ngmres->q);
31: PetscFree(ngmres->s);
32: PetscFree(ngmres->xnorms);
33: #if PETSC_USE_COMPLEX
34: PetscFree(ngmres->rwork);
35: #endif
36: PetscFree(ngmres->work);
37: PetscFree(snes->data);
38: return(0);
39: }
43: PetscErrorCode SNESSetUp_NGMRES(SNES snes)
44: {
45: SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data;
46: const char *optionsprefix;
47: PetscInt msize,hsize;
51: SNESSetWorkVecs(snes,5);
52: if (!ngmres->Xdot) {VecDuplicateVecs(snes->vec_sol,ngmres->msize,&ngmres->Xdot);}
53: if (!ngmres->Fdot) {VecDuplicateVecs(snes->vec_sol,ngmres->msize,&ngmres->Fdot);}
54: if (!ngmres->setup_called) {
55: msize = ngmres->msize; /* restart size */
56: hsize = msize * msize;
58: /* explicit least squares minimization solve */
59: PetscMalloc5(hsize,PetscScalar,&ngmres->h,
60: msize,PetscScalar,&ngmres->beta,
61: msize,PetscScalar,&ngmres->xi,
62: msize,PetscReal, &ngmres->fnorms,
63: hsize,PetscScalar,&ngmres->q);
64: if (ngmres->singlereduction) {
65: PetscMalloc(msize*sizeof(PetscReal),&ngmres->xnorms);
66: }
67: ngmres->nrhs = 1;
68: ngmres->lda = msize;
69: ngmres->ldb = msize;
70: PetscMalloc(msize*sizeof(PetscScalar),&ngmres->s);
71: PetscMemzero(ngmres->h, hsize*sizeof(PetscScalar));
72: PetscMemzero(ngmres->q, hsize*sizeof(PetscScalar));
73: PetscMemzero(ngmres->xi, msize*sizeof(PetscScalar));
74: PetscMemzero(ngmres->beta,msize*sizeof(PetscScalar));
75: ngmres->lwork = 12*msize;
76: #if PETSC_USE_COMPLEX
77: PetscMalloc(sizeof(PetscReal)*ngmres->lwork,&ngmres->rwork);
78: #endif
79: PetscMalloc(sizeof(PetscScalar)*ngmres->lwork,&ngmres->work);
80: }
82: /* linesearch setup */
83: SNESGetOptionsPrefix(snes,&optionsprefix);
85: if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) {
86: SNESLineSearchCreate(PetscObjectComm((PetscObject)snes),&ngmres->additive_linesearch);
87: SNESLineSearchSetSNES(ngmres->additive_linesearch,snes);
88: SNESLineSearchSetType(ngmres->additive_linesearch,SNESLINESEARCHL2);
89: SNESLineSearchAppendOptionsPrefix(ngmres->additive_linesearch,"additive_");
90: SNESLineSearchAppendOptionsPrefix(ngmres->additive_linesearch,optionsprefix);
91: SNESLineSearchSetFromOptions(ngmres->additive_linesearch);
92: }
94: ngmres->setup_called = PETSC_TRUE;
95: return(0);
96: }
100: PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes)
101: {
102: SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data;
104: PetscBool debug;
105: 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,NULL);
111: PetscOptionsEnum("-snes_ngmres_restart_type","Restart type","SNESNGMRESSetRestartType",SNESNGMRESRestartTypes,
112: (PetscEnum)ngmres->restart_type,(PetscEnum*)&ngmres->restart_type,NULL);
113: PetscOptionsBool("-snes_ngmres_candidate", "Use candidate storage", "SNES",ngmres->candidate,&ngmres->candidate,NULL);
114: PetscOptionsBool("-snes_ngmres_approxfunc","Linearly approximate the function", "SNES",ngmres->approxfunc,&ngmres->approxfunc,NULL);
115: PetscOptionsInt("-snes_ngmres_m", "Number of directions", "SNES",ngmres->msize,&ngmres->msize,NULL);
116: PetscOptionsInt("-snes_ngmres_restart", "Iterations before forced restart", "SNES",ngmres->restart_periodic,&ngmres->restart_periodic,NULL);
117: PetscOptionsInt("-snes_ngmres_restart_it", "Tolerance iterations before restart","SNES",ngmres->restart_it,&ngmres->restart_it,NULL);
118: PetscOptionsBool("-snes_ngmres_monitor", "Monitor actions of NGMRES", "SNES",ngmres->monitor ? PETSC_TRUE : PETSC_FALSE,&debug,NULL);
119: if (debug) {
120: ngmres->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
121: }
122: PetscOptionsReal("-snes_ngmres_gammaA", "Residual selection constant", "SNES",ngmres->gammaA,&ngmres->gammaA,NULL);
123: PetscOptionsReal("-snes_ngmres_gammaC", "Residual restart constant", "SNES",ngmres->gammaC,&ngmres->gammaC,NULL);
124: PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES",ngmres->epsilonB,&ngmres->epsilonB,NULL);
125: PetscOptionsReal("-snes_ngmres_deltaB", "Difference residual selection constant", "SNES",ngmres->deltaB,&ngmres->deltaB,NULL);
126: PetscOptionsBool("-snes_ngmres_single_reduction", "Aggregate reductions", "SNES",ngmres->singlereduction,&ngmres->singlereduction,NULL);
127: PetscOptionsTail();
128: if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA;
130: /* set the default type of the line search if the user hasn't already. */
131: if (!snes->linesearch) {
132: SNESGetLineSearch(snes,&linesearch);
133: SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC);
134: }
135: return(0);
136: }
140: PetscErrorCode SNESView_NGMRES(SNES snes,PetscViewer viewer)
141: {
142: SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data;
143: PetscBool iascii;
147: PetscObjectTypeCompare((PetscObject) viewer,PETSCVIEWERASCII,&iascii);
148: 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: }
158: PetscErrorCode SNESSolve_NGMRES(SNES snes)
159: {
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: /* coefficients and RHS to the minimization problem */
169: PetscReal fnorm,fMnorm,fAnorm;
170: PetscInt k,k_restart,l,ivec,restart_count = 0;
172: /* solution selection data */
173: PetscBool selectRestart;
174: PetscReal dnorm,dminnorm = 0.0;
175: PetscReal fminnorm,xnorm,ynorm;
177: SNESConvergedReason reason;
178: PetscBool lssucceed;
179: PetscErrorCode ierr;
182: /* variable initialization */
183: snes->reason = SNES_CONVERGED_ITERATING;
184: X = snes->vec_sol;
185: F = snes->vec_func;
186: B = snes->vec_rhs;
187: XA = snes->vec_sol_update;
188: FA = snes->work[0];
189: D = snes->work[1];
191: /* work for the line search */
192: Y = snes->work[2];
193: XM = snes->work[3];
194: FM = snes->work[4];
196: PetscObjectAMSTakeAccess((PetscObject)snes);
197: snes->iter = 0;
198: snes->norm = 0.;
199: PetscObjectAMSGrantAccess((PetscObject)snes);
201: /* initialization */
203: /* r = F(x) */
204: if (!snes->vec_func_init_set) {
205: SNESComputeFunction(snes,X,F);
206: if (snes->domainerror) {
207: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
208: return(0);
209: }
210: } else snes->vec_func_init_set = PETSC_FALSE;
212: if (!snes->norm_init_set) {
213: VecNorm(F,NORM_2,&fnorm);
214: if (PetscIsInfOrNanReal(fnorm)) {
215: snes->reason = SNES_DIVERGED_FNORM_NAN;
216: return(0);
217: }
218: } else {
219: fnorm = snes->norm_init;
220: snes->norm_init_set = PETSC_FALSE;
221: }
222: fminnorm = fnorm;
224: /* q_{00} = nu */
225: SNESNGMRESUpdateSubspace_Private(snes,0,0,F,fnorm,X);
227: PetscObjectAMSTakeAccess((PetscObject)snes);
228: snes->norm = fnorm;
229: PetscObjectAMSGrantAccess((PetscObject)snes);
230: SNESLogConvergenceHistory(snes,fnorm,0);
231: SNESMonitor(snes,0,fnorm);
232: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
233: if (snes->reason) return(0);
235: k_restart = 1;
236: l = 1;
237: for (k=1; k < snes->max_its+1; k++) {
238: /* select which vector of the stored subspace will be updated */
239: ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */
241: /* Computation of x^M */
242: if (snes->pc && snes->pcside == PC_RIGHT) {
243: VecCopy(X,XM);
244: SNESSetInitialFunction(snes->pc,F);
245: SNESSetInitialFunctionNorm(snes->pc,fnorm);
247: PetscLogEventBegin(SNES_NPCSolve,snes->pc,XM,B,0);
248: SNESSolve(snes->pc,B,XM);
249: PetscLogEventEnd(SNES_NPCSolve,snes->pc,XM,B,0);
251: SNESGetConvergedReason(snes->pc,&reason);
252: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
253: snes->reason = SNES_DIVERGED_INNER;
254: return(0);
255: }
256: SNESGetFunction(snes->pc,&FPC,NULL,NULL);
257: VecCopy(FPC,FM);
258: SNESGetFunctionNorm(snes->pc,&fMnorm);
259: } else {
260: /* no preconditioner -- just take gradient descent with line search */
261: VecCopy(F,Y);
262: VecCopy(F,FM);
263: VecCopy(X,XM);
265: fMnorm = fnorm;
267: SNESLineSearchApply(snes->linesearch,XM,FM,&fMnorm,Y);
268: SNESLineSearchGetSuccess(snes->linesearch,&lssucceed);
269: if (!lssucceed) {
270: if (++snes->numFailures >= snes->maxFailures) {
271: snes->reason = SNES_DIVERGED_LINE_SEARCH;
272: return(0);
273: }
274: }
275: }
276: SNESNGMRESFormCombinedSolution_Private(snes,l,XM,FM,fMnorm,X,XA,FA);
277: /* r = F(x) */
278: if (fminnorm > fMnorm) fminnorm = fMnorm; /* the minimum norm is now of F^M */
280: /* differences for selection and restart */
281: if (ngmres->restart_type == SNES_NGMRES_RESTART_DIFFERENCE || ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) {
282: SNESNGMRESCalculateDifferences_Private(snes,l,X,F,XM,FM,XA,FA,D,&dnorm,&dminnorm,&fAnorm);
283: } else {
284: VecNorm(FA,NORM_2,&fAnorm);
285: }
286: if (PetscIsInfOrNanReal(fAnorm)) {
287: snes->reason = SNES_DIVERGED_FNORM_NAN;
288: return(0);
289: }
291: /* combination (additive) or selection (multiplicative) of the N-GMRES solution */
292: SNESNGMRESSelect_Private(snes,k_restart,XM,FM,fMnorm,XA,FA,fAnorm,dnorm,fminnorm,dminnorm,X,F,Y,&fnorm);
293: selectRestart = PETSC_FALSE;
294: if (ngmres->restart_type == SNES_NGMRES_RESTART_DIFFERENCE) {
295: SNESNGMRESSelectRestart_Private(snes,l,fAnorm,dnorm,fminnorm,dminnorm,&selectRestart);
296: /* if the restart conditions persist for more than restart_it iterations, restart. */
297: if (selectRestart) restart_count++;
298: else restart_count = 0;
299: } else if (ngmres->restart_type == SNES_NGMRES_RESTART_PERIODIC) {
300: if (k_restart > ngmres->restart_periodic) {
301: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor,"periodic restart after %D iterations\n",k_restart);
302: restart_count = ngmres->restart_it;
303: }
304: }
305: /* restart after restart conditions have persisted for a fixed number of iterations */
306: if (restart_count >= ngmres->restart_it) {
307: if (ngmres->monitor) {
308: PetscViewerASCIIPrintf(ngmres->monitor,"Restarted at iteration %d\n",k_restart);
309: }
310: restart_count = 0;
311: k_restart = 1;
312: l = 1;
313: /* q_{00} = nu */
314: if (ngmres->candidate) {
315: SNESNGMRESUpdateSubspace_Private(snes,0,0,FM,fMnorm,XM);
316: } else {
317: SNESNGMRESUpdateSubspace_Private(snes,0,0,F,fMnorm,X);
318: }
319: } else {
320: /* select the current size of the subspace */
321: if (l < ngmres->msize) l++;
322: k_restart++;
323: /* place the current entry in the list of previous entries */
324: if (ngmres->candidate) {
325: if (fminnorm > fMnorm) fminnorm = fMnorm;
326: SNESNGMRESUpdateSubspace_Private(snes,ivec,l,FM,fMnorm,XM);
327: } else {
328: if (fminnorm > fnorm) fminnorm = fnorm;
329: SNESNGMRESUpdateSubspace_Private(snes,ivec,l,F,fnorm,X);
330: }
331: }
333: PetscObjectAMSTakeAccess((PetscObject)snes);
334: snes->iter = k;
335: snes->norm = fnorm;
336: PetscObjectAMSGrantAccess((PetscObject)snes);
337: SNESLogConvergenceHistory(snes,snes->norm,snes->iter);
338: SNESMonitor(snes,snes->iter,snes->norm);
339: VecNormBegin(Y,NORM_2,&ynorm);
340: VecNormBegin(X,NORM_2,&xnorm);
341: VecNormEnd(Y,NORM_2,&ynorm);
342: VecNormEnd(X,NORM_2,&xnorm);
343: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
344: if (snes->reason) return(0);
345: }
346: snes->reason = SNES_DIVERGED_MAX_IT;
347: return(0);
348: }
352: /*@
353: SNESNGMRESSetRestartType - Sets the restart type for SNESNGMRES.
355: Logically Collective on SNES
357: Input Parameters:
358: + snes - the iterative context
359: - rtype - restart type
361: Options Database:
362: + -snes_ngmres_restart_type<difference,periodic,none> - set the restart type
363: - -snes_ngmres_restart[30] - sets the number of iterations before restart for periodic
365: Level: intermediate
367: SNESNGMRESRestartTypes:
368: + SNES_NGMRES_RESTART_NONE - never restart
369: . SNES_NGMRES_RESTART_DIFFERENCE - restart based upon difference criteria
370: - SNES_NGMRES_RESTART_PERIODIC - restart after a fixed number of iterations
372: Notes:
373: The default line search used is the L2 line search and it requires two additional function evaluations.
375: .keywords: SNES, SNESNGMRES, restart, type, set SNESLineSearch
376: @*/
377: PetscErrorCode SNESNGMRESSetRestartType(SNES snes,SNESNGMRESRestartType rtype)
378: {
383: PetscTryMethod(snes,"SNESNGMRESSetRestartType_C",(SNES,SNESNGMRESRestartType),(snes,rtype));
384: return(0);
385: }
389: /*@
390: SNESNGMRESSetSelectType - Sets the selection type for SNESNGMRES. This determines how the candidate solution and
391: combined solution are used to create the next iterate.
393: Logically Collective on SNES
395: Input Parameters:
396: + snes - the iterative context
397: - stype - selection type
399: Options Database:
400: . -snes_ngmres_select_type<difference,none,linesearch>
402: Level: intermediate
404: SNESNGMRESSelectTypes:
405: + SNES_NGMRES_SELECT_NONE - choose the combined solution all the time
406: . SNES_NGMRES_SELECT_DIFFERENCE - choose based upon the selection criteria
407: - SNES_NGMRES_SELECT_LINESEARCH - choose based upon line search combination
409: Notes:
410: The default line search used is the L2 line search and it requires two additional function evaluations.
412: .keywords: SNES, SNESNGMRES, selection, type, set SNESLineSearch
413: @*/
414: PetscErrorCode SNESNGMRESSetSelectType(SNES snes,SNESNGMRESSelectType stype)
415: {
420: PetscTryMethod(snes,"SNESNGMRESSetSelectType_C",(SNES,SNESNGMRESSelectType),(snes,stype));
421: return(0);
422: }
426: PetscErrorCode SNESNGMRESSetSelectType_NGMRES(SNES snes,SNESNGMRESSelectType stype)
427: {
428: SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
431: ngmres->select_type = stype;
432: return(0);
433: }
437: PetscErrorCode SNESNGMRESSetRestartType_NGMRES(SNES snes,SNESNGMRESRestartType rtype)
438: {
439: SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
442: ngmres->restart_type = rtype;
443: return(0);
444: }
446: /*MC
447: SNESNGMRES - The Nonlinear Generalized Minimum Residual method.
449: Level: beginner
451: Options Database:
452: + -snes_ngmres_select_type<difference,none,linesearch> - choose the select between candidate and combined solution
453: . -snes_ngmres_restart_type<difference,none,periodic> - choose the restart conditions
454: . -snes_ngmres_candidate - Use NGMRES variant which combines candidate solutions instead of actual solutions
455: . -snes_ngmres_m - Number of stored previous solutions and residuals
456: . -snes_ngmres_restart_it - Number of iterations the restart conditions hold before restart
457: . -snes_ngmres_gammaA - Residual tolerance for solution select between the candidate and combination
458: . -snes_ngmres_gammaC - Residual tolerance for restart
459: . -snes_ngmres_epsilonB - Difference tolerance between subsequent solutions triggering restart
460: . -snes_ngmres_deltaB - Difference tolerance between residuals triggering restart
461: . -snes_ngmres_monitor - Prints relevant information about the ngmres iteration
462: . -snes_linesearch_type <basic,l2,cp> - Line search type used for the default smoother
463: - -additive_snes_linesearch_type - linesearch type used to select between the candidate and combined solution with additive select type
465: Notes:
467: The N-GMRES method combines m previous solutions into a minimum-residual solution by solving a small linearized
468: optimization problem at each iteration.
470: References:
472: "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio,
473: SIAM Journal on Scientific Computing, 21(5), 2000.
475: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
476: M*/
480: PETSC_EXTERN PetscErrorCode SNESCreate_NGMRES(SNES snes)
481: {
482: SNES_NGMRES *ngmres;
486: snes->ops->destroy = SNESDestroy_NGMRES;
487: snes->ops->setup = SNESSetUp_NGMRES;
488: snes->ops->setfromoptions = SNESSetFromOptions_NGMRES;
489: snes->ops->view = SNESView_NGMRES;
490: snes->ops->solve = SNESSolve_NGMRES;
491: snes->ops->reset = SNESReset_NGMRES;
493: snes->usespc = PETSC_TRUE;
494: snes->usesksp = PETSC_FALSE;
496: PetscNewLog(snes,SNES_NGMRES,&ngmres);
497: snes->data = (void*) ngmres;
498: ngmres->msize = 30;
500: if (!snes->tolerancesset) {
501: snes->max_funcs = 30000;
502: snes->max_its = 10000;
503: }
505: ngmres->candidate = PETSC_FALSE;
507: ngmres->additive_linesearch = NULL;
508: ngmres->approxfunc = PETSC_FALSE;
509: ngmres->restart_it = 2;
510: ngmres->restart_periodic = 30;
511: ngmres->gammaA = 2.0;
512: ngmres->gammaC = 2.0;
513: ngmres->deltaB = 0.9;
514: ngmres->epsilonB = 0.1;
516: ngmres->restart_type = SNES_NGMRES_RESTART_DIFFERENCE;
517: ngmres->select_type = SNES_NGMRES_SELECT_DIFFERENCE;
519: PetscObjectComposeFunction((PetscObject)snes,"SNESNGMRESSetSelectType_C",SNESNGMRESSetSelectType_NGMRES);
520: PetscObjectComposeFunction((PetscObject)snes,"SNESNGMRESSetRestartType_C",SNESNGMRESSetRestartType_NGMRES);
521: return(0);
522: }