Actual source code: anderson.c

  1: #include <../src/snes/impls/ngmres/snesngmres.h>

  3: static PetscErrorCode SNESSetFromOptions_Anderson(SNES snes, PetscOptionItems *PetscOptionsObject)
  4: {
  5:   SNES_NGMRES *ngmres  = (SNES_NGMRES *)snes->data;
  6:   PetscBool    monitor = PETSC_FALSE;

  8:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES NGMRES options");
  9:   PetscOptionsInt("-snes_anderson_m", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, NULL);
 10:   PetscOptionsReal("-snes_anderson_beta", "Mixing parameter", "SNES", ngmres->andersonBeta, &ngmres->andersonBeta, NULL);
 11:   PetscOptionsInt("-snes_anderson_restart", "Iterations before forced restart", "SNES", ngmres->restart_periodic, &ngmres->restart_periodic, NULL);
 12:   PetscOptionsInt("-snes_anderson_restart_it", "Tolerance iterations before restart", "SNES", ngmres->restart_it, &ngmres->restart_it, NULL);
 13:   PetscOptionsEnum("-snes_anderson_restart_type", "Restart type", "SNESNGMRESSetRestartType", SNESNGMRESRestartTypes, (PetscEnum)ngmres->restart_type, (PetscEnum *)&ngmres->restart_type, NULL);
 14:   PetscOptionsBool("-snes_anderson_monitor", "Monitor steps of Anderson Mixing", "SNES", ngmres->monitor ? PETSC_TRUE : PETSC_FALSE, &monitor, NULL);
 15:   if (monitor) ngmres->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
 16:   PetscOptionsHeadEnd();
 17:   return 0;
 18: }

 20: static PetscErrorCode SNESSolve_Anderson(SNES snes)
 21: {
 22:   SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
 23:   /* present solution, residual, and preconditioned residual */
 24:   Vec X, F, B, D;
 25:   /* candidate linear combination answers */
 26:   Vec XA, FA, XM, FM;

 28:   /* coefficients and RHS to the minimization problem */
 29:   PetscReal           fnorm, fMnorm, fAnorm;
 30:   PetscReal           xnorm, ynorm;
 31:   PetscReal           dnorm, dminnorm = 0.0, fminnorm;
 32:   PetscInt            restart_count = 0;
 33:   PetscInt            k, k_restart, l, ivec;
 34:   PetscBool           selectRestart;
 35:   SNESConvergedReason reason;


 39:   PetscCitationsRegister(SNESCitation, &SNEScite);
 40:   /* variable initialization */
 41:   snes->reason = SNES_CONVERGED_ITERATING;
 42:   X            = snes->vec_sol;
 43:   F            = snes->vec_func;
 44:   B            = snes->vec_rhs;
 45:   XA           = snes->vec_sol_update;
 46:   FA           = snes->work[0];
 47:   D            = snes->work[1];

 49:   /* work for the line search */
 50:   XM = snes->work[3];
 51:   FM = snes->work[4];

 53:   PetscObjectSAWsTakeAccess((PetscObject)snes);
 54:   snes->iter = 0;
 55:   snes->norm = 0.;
 56:   PetscObjectSAWsGrantAccess((PetscObject)snes);

 58:   /* initialization */

 60:   /* r = F(x) */

 62:   if (snes->npc && snes->npcside == PC_LEFT) {
 63:     SNESApplyNPC(snes, X, NULL, F);
 64:     SNESGetConvergedReason(snes->npc, &reason);
 65:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 66:       snes->reason = SNES_DIVERGED_INNER;
 67:       return 0;
 68:     }
 69:     VecNorm(F, NORM_2, &fnorm);
 70:   } else {
 71:     if (!snes->vec_func_init_set) {
 72:       SNESComputeFunction(snes, X, F);
 73:     } else snes->vec_func_init_set = PETSC_FALSE;

 75:     VecNorm(F, NORM_2, &fnorm);
 76:     SNESCheckFunctionNorm(snes, fnorm);
 77:   }
 78:   fminnorm = fnorm;

 80:   PetscObjectSAWsTakeAccess((PetscObject)snes);
 81:   snes->norm = fnorm;
 82:   PetscObjectSAWsGrantAccess((PetscObject)snes);
 83:   SNESLogConvergenceHistory(snes, fnorm, 0);
 84:   SNESMonitor(snes, 0, fnorm);
 85:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
 86:   if (snes->reason) return 0;

 88:   k_restart = 0;
 89:   l         = 0;
 90:   ivec      = 0;
 91:   for (k = 1; k < snes->max_its + 1; k++) {
 92:     /* select which vector of the stored subspace will be updated */
 93:     if (snes->npc && snes->npcside == PC_RIGHT) {
 94:       VecCopy(X, XM);
 95:       SNESSetInitialFunction(snes->npc, F);

 97:       PetscLogEventBegin(SNES_NPCSolve, snes->npc, XM, B, 0);
 98:       SNESSolve(snes->npc, B, XM);
 99:       PetscLogEventEnd(SNES_NPCSolve, snes->npc, XM, B, 0);

101:       SNESGetConvergedReason(snes->npc, &reason);
102:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
103:         snes->reason = SNES_DIVERGED_INNER;
104:         return 0;
105:       }
106:       SNESGetNPCFunction(snes, FM, &fMnorm);
107:       if (ngmres->andersonBeta != 1.0) VecAXPBY(XM, (1.0 - ngmres->andersonBeta), ngmres->andersonBeta, X);
108:     } else {
109:       VecCopy(F, FM);
110:       VecCopy(X, XM);
111:       VecAXPY(XM, -ngmres->andersonBeta, FM);
112:       fMnorm = fnorm;
113:     }

115:     SNESNGMRESFormCombinedSolution_Private(snes, ivec, l, XM, FM, fMnorm, X, XA, FA);
116:     ivec = k_restart % ngmres->msize;
117:     if (ngmres->restart_type == SNES_NGMRES_RESTART_DIFFERENCE) {
118:       SNESNGMRESNorms_Private(snes, l, X, F, XM, FM, XA, FA, D, &dnorm, &dminnorm, NULL, NULL, NULL, &xnorm, &fAnorm, &ynorm);
119:       SNESNGMRESSelectRestart_Private(snes, l, fMnorm, fnorm, dnorm, fminnorm, dminnorm, &selectRestart);
120:       /* if the restart conditions persist for more than restart_it iterations, restart. */
121:       if (selectRestart) restart_count++;
122:       else restart_count = 0;
123:     } else if (ngmres->restart_type == SNES_NGMRES_RESTART_PERIODIC) {
124:       SNESNGMRESNorms_Private(snes, l, X, F, XM, FM, XA, FA, D, NULL, NULL, NULL, NULL, NULL, &xnorm, &fAnorm, &ynorm);
125:       if (k_restart > ngmres->restart_periodic) {
126:         if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "periodic restart after %" PetscInt_FMT " iterations\n", k_restart);
127:         restart_count = ngmres->restart_it;
128:       }
129:     } else {
130:       SNESNGMRESNorms_Private(snes, l, X, F, XM, FM, XA, FA, D, NULL, NULL, NULL, NULL, NULL, &xnorm, &fAnorm, &ynorm);
131:     }
132:     /* restart after restart conditions have persisted for a fixed number of iterations */
133:     if (restart_count >= ngmres->restart_it) {
134:       if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %" PetscInt_FMT "\n", k_restart);
135:       restart_count = 0;
136:       k_restart     = 0;
137:       l             = 0;
138:       ivec          = 0;
139:     } else {
140:       if (l < ngmres->msize) l++;
141:       k_restart++;
142:       SNESNGMRESUpdateSubspace_Private(snes, ivec, l, FM, fnorm, XM);
143:     }

145:     fnorm = fAnorm;
146:     if (fminnorm > fnorm) fminnorm = fnorm;

148:     VecCopy(XA, X);
149:     VecCopy(FA, F);

151:     PetscObjectSAWsTakeAccess((PetscObject)snes);
152:     snes->iter  = k;
153:     snes->norm  = fnorm;
154:     snes->xnorm = xnorm;
155:     snes->ynorm = ynorm;
156:     PetscObjectSAWsGrantAccess((PetscObject)snes);
157:     SNESLogConvergenceHistory(snes, snes->norm, snes->iter);
158:     SNESMonitor(snes, snes->iter, snes->norm);
159:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
160:     if (snes->reason) return 0;
161:   }
162:   snes->reason = SNES_DIVERGED_MAX_IT;
163:   return 0;
164: }

166: /*MC
167:   SNESANDERSON - Anderson Mixing nonlinear solver

169:    Level: beginner

171:    Options Database Keys:
172: +  -snes_anderson_m                - Number of stored previous solutions and residuals
173: .  -snes_anderson_beta             - Anderson mixing parameter
174: .  -snes_anderson_restart_type     - Type of restart (see SNESNGMRES)
175: .  -snes_anderson_restart_it       - Number of iterations of restart conditions before restart
176: .  -snes_anderson_restart          - Number of iterations before periodic restart
177: -  -snes_anderson_monitor          - Prints relevant information about the ngmres iteration

179:    Notes:
180:    The Anderson Mixing method combines m previous solutions into a minimum-residual solution by solving a small linearized
181:    optimization problem at each iteration.

183:    Very similar to the `SNESNGMRES` algorithm.

185:    References:
186: +  * -  D. G. Anderson. Iterative procedures for nonlinear integral equations.
187:     J. Assoc. Comput. Mach., 12, 1965."
188: -  * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu,"Composing Scalable Nonlinear Algebraic Solvers",
189:    SIAM Review, 57(4), 2015

191: .seealso: `SNESNGMRES`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESType`
192: M*/

194: PETSC_EXTERN PetscErrorCode SNESCreate_Anderson(SNES snes)
195: {
196:   SNES_NGMRES   *ngmres;
197:   SNESLineSearch linesearch;

199:   snes->ops->destroy        = SNESDestroy_NGMRES;
200:   snes->ops->setup          = SNESSetUp_NGMRES;
201:   snes->ops->setfromoptions = SNESSetFromOptions_Anderson;
202:   snes->ops->view           = SNESView_NGMRES;
203:   snes->ops->solve          = SNESSolve_Anderson;
204:   snes->ops->reset          = SNESReset_NGMRES;

206:   snes->usesnpc = PETSC_TRUE;
207:   snes->usesksp = PETSC_FALSE;
208:   snes->npcside = PC_RIGHT;

210:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

212:   PetscNew(&ngmres);
213:   snes->data    = (void *)ngmres;
214:   ngmres->msize = 30;

216:   if (!snes->tolerancesset) {
217:     snes->max_funcs = 30000;
218:     snes->max_its   = 10000;
219:   }

221:   SNESGetLineSearch(snes, &linesearch);
222:   if (!((PetscObject)linesearch)->type_name) SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);

224:   ngmres->additive_linesearch = NULL;
225:   ngmres->approxfunc          = PETSC_FALSE;
226:   ngmres->restart_type        = SNES_NGMRES_RESTART_NONE;
227:   ngmres->restart_it          = 2;
228:   ngmres->restart_periodic    = 30;
229:   ngmres->gammaA              = 2.0;
230:   ngmres->gammaC              = 2.0;
231:   ngmres->deltaB              = 0.9;
232:   ngmres->epsilonB            = 0.1;

234:   ngmres->andersonBeta = 1.0;
235:   return 0;
236: }