Actual source code: qn.c

petsc-3.7.3 2016-07-24
Report Typos and Errors
  1: #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
  2: #include <petscdm.h>

  4: #define H(i,j)  qn->dXdFmat[i*qn->m + j]

  6: const char *const SNESQNScaleTypes[] =        {"DEFAULT","NONE","SHANNO","LINESEARCH","JACOBIAN","SNESQNScaleType","SNES_QN_SCALING_",0};
  7: const char *const SNESQNRestartTypes[] =      {"DEFAULT","NONE","POWELL","PERIODIC","SNESQNRestartType","SNES_QN_RESTART_",0};
  8: const char *const SNESQNTypes[] =             {"LBFGS","BROYDEN","BADBROYDEN","SNESQNType","SNES_QN_",0};

 10: typedef struct {
 11:   Vec               *U;                   /* Stored past states (vary from method to method) */
 12:   Vec               *V;                   /* Stored past states (vary from method to method) */
 13:   PetscInt          m;                    /* The number of kept previous steps */
 14:   PetscReal         *lambda;              /* The line search history of the method */
 15:   PetscReal         *norm;                /* norms of the steps */
 16:   PetscScalar       *alpha, *beta;
 17:   PetscScalar       *dXtdF, *dFtdX, *YtdX;
 18:   PetscBool         singlereduction;      /* Aggregated reduction implementation */
 19:   PetscScalar       *dXdFmat;             /* A matrix of values for dX_i dot dF_j */
 20:   PetscViewer       monitor;
 21:   PetscReal         powell_gamma;         /* Powell angle restart condition */
 22:   PetscReal         scaling;              /* scaling of H0 */
 23:   SNESQNType        type;                 /* the type of quasi-newton method used */
 24:   SNESQNScaleType   scale_type;           /* the type of scaling used */
 25:   SNESQNRestartType restart_type;         /* determine the frequency and type of restart conditions */
 26: } SNES_QN;

 30: PetscErrorCode SNESQNApply_Broyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D)
 31: {
 32:   PetscErrorCode     ierr;
 33:   SNES_QN            *qn = (SNES_QN*)snes->data;
 34:   Vec                W   = snes->work[3];
 35:   Vec                *U  = qn->U;
 36:   PetscInt           m = qn->m;
 37:   PetscInt           k,i,j,lits,l = m;
 38:   PetscReal          unorm,a,b;
 39:   PetscReal          *lambda=qn->lambda;
 40:   PetscScalar        gdot;
 41:   PetscReal          udot;

 44:   if (it < m) l = it;
 45:   if (it > 0) {
 46:     k = (it-1)%l;
 47:     SNESLineSearchGetLambda(snes->linesearch,&lambda[k]);
 48:     VecCopy(Xold,U[k]);
 49:     VecAXPY(U[k],-1.0,X);
 50:     if (qn->monitor) {
 51:       PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
 52:       PetscViewerASCIIPrintf(qn->monitor, "scaling vector %D of %D by lambda: %14.12e \n",k,l,(double)lambda[k]);
 53:       PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
 54:     }
 55:   }
 56:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
 57:     KSPSolve(snes->ksp,D,W);
 58:     SNESCheckKSPSolve(snes);
 59:     KSPGetIterationNumber(snes->ksp,&lits);
 60:     snes->linear_its += lits;
 61:     VecCopy(W,Y);
 62:   } else {
 63:     VecCopy(D,Y);
 64:     VecScale(Y,qn->scaling);
 65:   }

 67:   /* inward recursion starting at the first update and working forward */
 68:   for (i = 0; i < l-1; i++) {
 69:     j = (it+i-l)%l;
 70:     k = (it+i-l+1)%l;
 71:     VecNorm(U[j],NORM_2,&unorm);
 72:     VecDot(U[j],Y,&gdot);
 73:     unorm *= unorm;
 74:     udot = PetscRealPart(gdot);
 75:     a = (lambda[j]/lambda[k]);
 76:     b = -(1.-lambda[j]);
 77:     a *= udot/unorm;
 78:     b *= udot/unorm;
 79:     VecAXPBYPCZ(Y,a,b,1.,U[k],U[j]);

 81:     if (qn->monitor) {
 82:       PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
 83:       PetscViewerASCIIPrintf(qn->monitor, "using vector %D and %D, gdot: %14.12e\n",k,j,(double)PetscRealPart(gdot));
 84:       PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
 85:     }
 86:   }
 87:   if (l > 0) {
 88:     k = (it-1)%l;
 89:     VecDot(U[k],Y,&gdot);
 90:     VecNorm(U[k],NORM_2,&unorm);
 91:     unorm *= unorm;
 92:     udot = PetscRealPart(gdot);
 93:     a = unorm/(unorm-lambda[k]*udot);
 94:     b = -(1.-lambda[k])*udot/(unorm-lambda[k]*udot);
 95:     if (qn->monitor) {
 96:       PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
 97:       PetscViewerASCIIPrintf(qn->monitor, "using vector %D: a: %14.12e b: %14.12e \n",k,(double)a,(double)b);
 98:       PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
 99:     }
100:     VecAXPBY(Y,b,a,U[k]);
101:   }
102:   l = m;
103:   if (it+1<m)l=it+1;
104:   k = it%l;
105:   if (qn->monitor) {
106:     PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
107:     PetscViewerASCIIPrintf(qn->monitor, "setting vector %D of %D\n",k,l);
108:     PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
109:   }
110:   return(0);
111: }

115: PetscErrorCode SNESQNApply_BadBroyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D,Vec Dold)
116: {
118:   SNES_QN        *qn = (SNES_QN*)snes->data;
119:   Vec            W   = snes->work[3];
120:   Vec            *U  = qn->U;
121:   Vec            *T  = qn->V;

123:   /* ksp thing for Jacobian scaling */
124:   PetscInt           h,k,j,i,lits;
125:   PetscInt           m = qn->m;
126:   PetscScalar        gdot,udot;
127:   PetscInt           l = m;

130:   if (it < m) l = it;
131:   VecCopy(D,Y);
132:   if (l > 0) {
133:     k    = (it-1)%l;
134:     SNESLineSearchGetLambda(snes->linesearch,&qn->lambda[k]);
135:     VecCopy(Dold,U[k]);
136:     VecAXPY(U[k],-1.0,D);
137:     VecCopy(Xold,T[k]);
138:     VecAXPY(T[k],-1.0,X);
139:   }

141:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
142:     KSPSolve(snes->ksp,Y,W);
143:     SNESCheckKSPSolve(snes);
144:     KSPGetIterationNumber(snes->ksp,&lits);
145:     snes->linear_its += lits;
146:     VecCopy(W,Y);
147:   } else {
148:     VecScale(Y,qn->scaling);
149:   }

151:   /* inward recursion starting at the first update and working forward */
152:   if (l > 0) {
153:     for (i = 0; i < l-1; i++) {
154:       j    = (it+i-l)%l;
155:       k    = (it+i-l+1)%l;
156:       h    = (it-1)%l;
157:       VecDotBegin(U[j],U[h],&gdot);
158:       VecDotBegin(U[j],U[j],&udot);
159:       VecDotEnd(U[j],U[h],&gdot);
160:       VecDotEnd(U[j],U[j],&udot);
161:       VecAXPY(Y,PetscRealPart(gdot)/PetscRealPart(udot),T[k]);
162:       VecAXPY(Y,-(1.-qn->lambda[k])*PetscRealPart(gdot)/PetscRealPart(udot),T[j]);
163:       if (qn->monitor) {
164:         PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
165:         PetscViewerASCIIPrintf(qn->monitor, "it: %D k: %D gdot: %14.12e\n", it, k, (double)PetscRealPart(gdot));
166:         PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
167:       }
168:     }
169:   }
170:   return(0);
171: }

175: PetscErrorCode SNESQNApply_LBFGS(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D,Vec Dold)
176: {
178:   SNES_QN        *qn    = (SNES_QN*)snes->data;
179:   Vec            W      = snes->work[3];
180:   Vec            *dX    = qn->U;
181:   Vec            *dF    = qn->V;
182:   PetscScalar    *alpha = qn->alpha;
183:   PetscScalar    *beta  = qn->beta;
184:   PetscScalar    *dXtdF = qn->dXtdF;
185:   PetscScalar    *dFtdX = qn->dFtdX;
186:   PetscScalar    *YtdX  = qn->YtdX;

188:   /* ksp thing for Jacobian scaling */
189:   PetscInt           k,i,j,g,lits;
190:   PetscInt           m = qn->m;
191:   PetscScalar        t;
192:   PetscInt           l = m;

195:   if (it < m) l = it;
196:   VecCopy(D,Y);
197:   if (it > 0) {
198:     k    = (it - 1) % l;
199:     VecCopy(D,dF[k]);
200:     VecAXPY(dF[k], -1.0, Dold);
201:     VecCopy(X, dX[k]);
202:     VecAXPY(dX[k], -1.0, Xold);
203:     if (qn->singlereduction) {
204:       VecMDotBegin(dF[k],l,dX,dXtdF);
205:       VecMDotBegin(dX[k],l,dF,dFtdX);
206:       VecMDotBegin(Y,l,dX,YtdX);
207:       VecMDotEnd(dF[k],l,dX,dXtdF);
208:       VecMDotEnd(dX[k],l,dF,dFtdX);
209:       VecMDotEnd(Y,l,dX,YtdX);
210:       for (j = 0; j < l; j++) {
211:         H(k, j) = dFtdX[j];
212:         H(j, k) = dXtdF[j];
213:       }
214:       /* copy back over to make the computation of alpha and beta easier */
215:       for (j = 0; j < l; j++) dXtdF[j] = H(j, j);
216:     } else {
217:       VecDot(dX[k], dF[k], &dXtdF[k]);
218:     }
219:     if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
220:       SNESLineSearchGetLambda(snes->linesearch,&qn->scaling);
221:     }
222:   }

224:   /* outward recursion starting at iteration k's update and working back */
225:   for (i=0; i<l; i++) {
226:     k = (it-i-1)%l;
227:     if (qn->singlereduction) {
228:       /* construct t = dX[k] dot Y as Y_0 dot dX[k] + sum(-alpha[j]dX[k]dF[j]) */
229:       t = YtdX[k];
230:       for (j=0; j<i; j++) {
231:         g  = (it-j-1)%l;
232:         t -= alpha[g]*H(k, g);
233:       }
234:       alpha[k] = t/H(k,k);
235:     } else {
236:       VecDot(dX[k],Y,&t);
237:       alpha[k] = t/dXtdF[k];
238:     }
239:     if (qn->monitor) {
240:       PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
241:       PetscViewerASCIIPrintf(qn->monitor, "it: %D k: %D alpha:        %14.12e\n", it, k, (double)PetscRealPart(alpha[k]));
242:       PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
243:     }
244:     VecAXPY(Y,-alpha[k],dF[k]);
245:   }

247:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
248:     KSPSolve(snes->ksp,Y,W);
249:     SNESCheckKSPSolve(snes);
250:     KSPGetIterationNumber(snes->ksp,&lits);
251:     snes->linear_its += lits;
252:     VecCopy(W, Y);
253:   } else {
254:     VecScale(Y, qn->scaling);
255:   }
256:   if (qn->singlereduction) {
257:     VecMDot(Y,l,dF,YtdX);
258:   }
259:   /* inward recursion starting at the first update and working forward */
260:   for (i = 0; i < l; i++) {
261:     k = (it + i - l) % l;
262:     if (qn->singlereduction) {
263:       t = YtdX[k];
264:       for (j = 0; j < i; j++) {
265:         g  = (it + j - l) % l;
266:         t += (alpha[g] - beta[g])*H(g, k);
267:       }
268:       beta[k] = t / H(k, k);
269:     } else {
270:       VecDot(dF[k], Y, &t);
271:       beta[k] = t / dXtdF[k];
272:     }
273:     VecAXPY(Y, (alpha[k] - beta[k]), dX[k]);
274:     if (qn->monitor) {
275:       PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
276:       PetscViewerASCIIPrintf(qn->monitor, "it: %D k: %D alpha - beta: %14.12e\n", it, k, (double)PetscRealPart(alpha[k] - beta[k]));
277:       PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
278:     }
279:   }
280:   return(0);
281: }

285: static PetscErrorCode SNESSolve_QN(SNES snes)
286: {
287:   PetscErrorCode       ierr;
288:   SNES_QN              *qn = (SNES_QN*) snes->data;
289:   Vec                  X,Xold;
290:   Vec                  F,W;
291:   Vec                  Y,D,Dold;
292:   PetscInt             i, i_r;
293:   PetscReal            fnorm,xnorm,ynorm,gnorm;
294:   SNESLineSearchReason lssucceed;
295:   PetscBool            powell,periodic;
296:   PetscScalar          DolddotD,DolddotDold;
297:   SNESConvergedReason  reason;

299:   /* basically just a regular newton's method except for the application of the Jacobian */

302:   if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

304:   PetscCitationsRegister(SNESCitation,&SNEScite);
305:   F    = snes->vec_func;                /* residual vector */
306:   Y    = snes->vec_sol_update;          /* search direction generated by J^-1D*/
307:   W    = snes->work[3];
308:   X    = snes->vec_sol;                 /* solution vector */
309:   Xold = snes->work[0];

311:   /* directions generated by the preconditioned problem with F_pre = F or x - M(x, b) */
312:   D    = snes->work[1];
313:   Dold = snes->work[2];

315:   snes->reason = SNES_CONVERGED_ITERATING;

317:   PetscObjectSAWsTakeAccess((PetscObject)snes);
318:   snes->iter = 0;
319:   snes->norm = 0.;
320:   PetscObjectSAWsGrantAccess((PetscObject)snes);

322:   if (snes->pc && snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
323:     SNESApplyNPC(snes,X,NULL,F);
324:     SNESGetConvergedReason(snes->pc,&reason);
325:     if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
326:       snes->reason = SNES_DIVERGED_INNER;
327:       return(0);
328:     }
329:     VecNorm(F,NORM_2,&fnorm);
330:   } else {
331:     if (!snes->vec_func_init_set) {
332:       SNESComputeFunction(snes,X,F);
333:     } else snes->vec_func_init_set = PETSC_FALSE;

335:     VecNorm(F,NORM_2,&fnorm);
336:     SNESCheckFunctionNorm(snes,fnorm);
337:   }
338:   if (snes->pc && snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
339:       SNESApplyNPC(snes,X,F,D);
340:       SNESGetConvergedReason(snes->pc,&reason);
341:       if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
342:         snes->reason = SNES_DIVERGED_INNER;
343:         return(0);
344:       }
345:   } else {
346:     VecCopy(F,D);
347:   }

349:   PetscObjectSAWsTakeAccess((PetscObject)snes);
350:   snes->norm = fnorm;
351:   PetscObjectSAWsGrantAccess((PetscObject)snes);
352:   SNESLogConvergenceHistory(snes,fnorm,0);
353:   SNESMonitor(snes,0,fnorm);

355:   /* test convergence */
356:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
357:   if (snes->reason) return(0);

359:   if (snes->pc && snes->pcside == PC_RIGHT) {
360:     PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,0,0);
361:     SNESSolve(snes->pc,snes->vec_rhs,X);
362:     PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,0,0);
363:     SNESGetConvergedReason(snes->pc,&reason);
364:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
365:       snes->reason = SNES_DIVERGED_INNER;
366:       return(0);
367:     }
368:     SNESGetNPCFunction(snes,F,&fnorm);
369:     VecCopy(F,D);
370:   }

372:   /* scale the initial update */
373:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
374:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
375:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
376:   }

378:   for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
379:     if (qn->scale_type == SNES_QN_SCALE_SHANNO && i_r > 0) {
380:       PetscScalar ff,xf;
381:       VecCopy(Dold,Y);
382:       VecCopy(Xold,W);
383:       VecAXPY(Y,-1.0,D);
384:       VecAXPY(W,-1.0,X);
385:       VecDotBegin(Y,Y,&ff);
386:       VecDotBegin(W,Y,&xf);
387:       VecDotEnd(Y,Y,&ff);
388:       VecDotEnd(W,Y,&xf);
389:       qn->scaling = PetscRealPart(xf)/PetscRealPart(ff);
390:       if (qn->monitor) {
391:         PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
392:         PetscViewerASCIIPrintf(qn->monitor, "Shanno scaling %D %g\n", i,(double)qn->scaling);
393:         PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
394:       }
395:     }
396:     switch (qn->type) {
397:     case SNES_QN_BADBROYDEN:
398:       SNESQNApply_BadBroyden(snes,i_r,Y,X,Xold,D,Dold);
399:       break;
400:     case SNES_QN_BROYDEN:
401:       SNESQNApply_Broyden(snes,i_r,Y,X,Xold,D);
402:       break;
403:     case SNES_QN_LBFGS:
404:       SNESQNApply_LBFGS(snes,i_r,Y,X,Xold,D,Dold);
405:       break;
406:     }
407:     /* line search for lambda */
408:     ynorm = 1; gnorm = fnorm;
409:     VecCopy(D, Dold);
410:     VecCopy(X, Xold);
411:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
412:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
413:     SNESLineSearchGetReason(snes->linesearch, &lssucceed);
414:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
415:     if (lssucceed) {
416:       if (++snes->numFailures >= snes->maxFailures) {
417:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
418:         break;
419:       }
420:     }
421:     if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
422:       SNESLineSearchGetLambda(snes->linesearch, &qn->scaling);
423:     }

425:     /* convergence monitoring */
426:     PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);

428:     if (snes->pc && snes->pcside == PC_RIGHT) {
429:       PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,0,0);
430:       SNESSolve(snes->pc,snes->vec_rhs,X);
431:       PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,0,0);
432:       SNESGetConvergedReason(snes->pc,&reason);
433:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
434:         snes->reason = SNES_DIVERGED_INNER;
435:         return(0);
436:       }
437:       SNESGetNPCFunction(snes,F,&fnorm);
438:     }

440:     SNESSetIterationNumber(snes, i+1);
441:     snes->norm = fnorm;

443:     SNESLogConvergenceHistory(snes,snes->norm,snes->iter);
444:     SNESMonitor(snes,snes->iter,snes->norm);
445:     /* set parameter for default relative tolerance convergence test */
446:     (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
447:     if (snes->reason) return(0);
448:     if (snes->pc && snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
449:       SNESApplyNPC(snes,X,F,D);
450:       SNESGetConvergedReason(snes->pc,&reason);
451:       if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
452:         snes->reason = SNES_DIVERGED_INNER;
453:         return(0);
454:       }
455:     } else {
456:       VecCopy(F, D);
457:     }
458:     powell = PETSC_FALSE;
459:     if (qn->restart_type == SNES_QN_RESTART_POWELL && i_r > 1) {
460:       /* check restart by Powell's Criterion: |F^T H_0 Fold| > powell_gamma * |Fold^T H_0 Fold| */
461:       if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
462:         MatMult(snes->jacobian_pre,Dold,W);
463:       } else {
464:         VecCopy(Dold,W);
465:       }
466:       VecDotBegin(W, Dold, &DolddotDold);
467:       VecDotBegin(W, D, &DolddotD);
468:       VecDotEnd(W, Dold, &DolddotDold);
469:       VecDotEnd(W, D, &DolddotD);
470:       if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE;
471:     }
472:     periodic = PETSC_FALSE;
473:     if (qn->restart_type == SNES_QN_RESTART_PERIODIC) {
474:       if (i_r>qn->m-1) periodic = PETSC_TRUE;
475:     }
476:     /* restart if either powell or periodic restart is satisfied. */
477:     if (powell || periodic) {
478:       if (qn->monitor) {
479:         PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
480:         if (powell) {
481:           PetscViewerASCIIPrintf(qn->monitor, "Powell restart! |%14.12e| > %6.4f*|%14.12e| i_r = %D\n", (double)PetscRealPart(DolddotD), (double)qn->powell_gamma, (double)PetscRealPart(DolddotDold),i_r);
482:         } else {
483:           PetscViewerASCIIPrintf(qn->monitor, "Periodic restart! i_r = %D\n", i_r);
484:         }
485:         PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
486:       }
487:       i_r = -1;
488:       /* general purpose update */
489:       if (snes->ops->update) {
490:         (*snes->ops->update)(snes, snes->iter);
491:       }
492:       if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
493:         SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
494:       }
495:     }
496:     /* general purpose update */
497:     if (snes->ops->update) {
498:       (*snes->ops->update)(snes, snes->iter);
499:     }
500:   }
501:   if (i == snes->max_its) {
502:     PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);
503:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
504:   }
505:   return(0);
506: }

510: static PetscErrorCode SNESSetUp_QN(SNES snes)
511: {
512:   SNES_QN        *qn = (SNES_QN*)snes->data;
514:   DM             dm;


518:   if (!snes->vec_sol) {
519:     SNESGetDM(snes,&dm);
520:     DMCreateGlobalVector(dm,&snes->vec_sol);
521:   }

523:   VecDuplicateVecs(snes->vec_sol, qn->m, &qn->U);
524:   if (qn->type != SNES_QN_BROYDEN) VecDuplicateVecs(snes->vec_sol, qn->m, &qn->V);
525:   PetscMalloc4(qn->m,&qn->alpha,qn->m,&qn->beta,qn->m,&qn->dXtdF,qn->m,&qn->lambda);

527:   if (qn->singlereduction) {
528:     PetscMalloc3(qn->m*qn->m,&qn->dXdFmat,qn->m,&qn->dFtdX,qn->m,&qn->YtdX);
529:   }
530:   SNESSetWorkVecs(snes,4);
531:   /* set method defaults */
532:   if (qn->scale_type == SNES_QN_SCALE_DEFAULT) {
533:     if (qn->type == SNES_QN_BADBROYDEN) {
534:       qn->scale_type = SNES_QN_SCALE_NONE;
535:     } else {
536:       qn->scale_type = SNES_QN_SCALE_SHANNO;
537:     }
538:   }
539:   if (qn->restart_type == SNES_QN_RESTART_DEFAULT) {
540:     if (qn->type == SNES_QN_LBFGS) {
541:       qn->restart_type = SNES_QN_RESTART_POWELL;
542:     } else {
543:       qn->restart_type = SNES_QN_RESTART_PERIODIC;
544:     }
545:   }

547:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
548:     SNESSetUpMatrices(snes);
549:   }
550:   if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) {snes->functype = SNES_FUNCTION_UNPRECONDITIONED;}
551:   return(0);
552: }

556: static PetscErrorCode SNESReset_QN(SNES snes)
557: {
559:   SNES_QN        *qn;

562:   if (snes->data) {
563:     qn = (SNES_QN*)snes->data;
564:     if (qn->U) {
565:       VecDestroyVecs(qn->m, &qn->U);
566:     }
567:     if (qn->V) {
568:       VecDestroyVecs(qn->m, &qn->V);
569:     }
570:     if (qn->singlereduction) {
571:       PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);
572:     }
573:     PetscFree4(qn->alpha,qn->beta,qn->dXtdF,qn->lambda);
574:   }
575:   return(0);
576: }

580: static PetscErrorCode SNESDestroy_QN(SNES snes)
581: {

585:   SNESReset_QN(snes);
586:   PetscFree(snes->data);
587:   PetscObjectComposeFunction((PetscObject)snes,"",NULL);
588:   return(0);
589: }

593: static PetscErrorCode SNESSetFromOptions_QN(PetscOptionItems *PetscOptionsObject,SNES snes)
594: {

596:   PetscErrorCode    ierr;
597:   SNES_QN           *qn    = (SNES_QN*)snes->data;
598:   PetscBool         monflg = PETSC_FALSE,flg;
599:   SNESLineSearch    linesearch;
600:   SNESQNRestartType rtype = qn->restart_type;
601:   SNESQNScaleType   stype = qn->scale_type;
602:   SNESQNType        qtype = qn->type;

605:   PetscOptionsHead(PetscOptionsObject,"SNES QN options");
606:   PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,NULL);
607:   PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance",          "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL);
608:   PetscOptionsBool("-snes_qn_monitor",         "Monitor for the QN methods",      "SNESQN", monflg, &monflg, NULL);
609:   PetscOptionsBool("-snes_qn_single_reduction", "Aggregate reductions",           "SNESQN", qn->singlereduction, &qn->singlereduction, NULL);
610:   PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);
611:   if (flg) SNESQNSetScaleType(snes,stype);

613:   PetscOptionsEnum("-snes_qn_restart_type","Restart type","SNESQNSetRestartType",SNESQNRestartTypes,(PetscEnum)rtype,(PetscEnum*)&rtype,&flg);
614:   if (flg) SNESQNSetRestartType(snes,rtype);

616:   PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes,(PetscEnum)qtype,(PetscEnum*)&qtype,&flg);
617:   if (flg) {SNESQNSetType(snes,qtype);}
618:   PetscOptionsTail();
619:   if (!snes->linesearch) {
620:     SNESGetLineSearch(snes, &linesearch);
621:     if (qn->type == SNES_QN_LBFGS) {
622:       SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);
623:     } else if (qn->type == SNES_QN_BROYDEN) {
624:       SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);
625:     } else {
626:       SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
627:     }
628:   }
629:   if (monflg) {
630:     qn->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
631:   }
632:   return(0);
633: }

637: static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
638: {
639:   SNES_QN        *qn    = (SNES_QN*)snes->data;
640:   PetscBool      iascii;

644:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
645:   if (iascii) {
646:     PetscViewerASCIIPrintf(viewer,"  QN type is %s, restart type is %s, scale type is %s\n",SNESQNTypes[qn->type],SNESQNRestartTypes[qn->restart_type],SNESQNScaleTypes[qn->scale_type]);
647:     PetscViewerASCIIPrintf(viewer,"  Stored subspace size: %D\n", qn->m);
648:     if (qn->singlereduction) {
649:       PetscViewerASCIIPrintf(viewer,"  Using the single reduction variant.\n");
650:     }
651:   }
652:   return(0);
653: }

657: /*@
658:     SNESQNSetRestartType - Sets the restart type for SNESQN.

660:     Logically Collective on SNES

662:     Input Parameters:
663: +   snes - the iterative context
664: -   rtype - restart type

666:     Options Database:
667: +   -snes_qn_restart_type <powell,periodic,none> - set the restart type
668: -   -snes_qn_m <m> - sets the number of stored updates and the restart period for periodic

670:     Level: intermediate

672:     SNESQNRestartTypes:
673: +   SNES_QN_RESTART_NONE - never restart
674: .   SNES_QN_RESTART_POWELL - restart based upon descent criteria
675: -   SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations

677: .keywords: SNES, SNESQN, restart, type, set SNESLineSearch, SNESQNRestartType
678: @*/
679: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
680: {

685:   PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));
686:   return(0);
687: }

691: /*@
692:     SNESQNSetScaleType - Sets the scaling type for the inner inverse Jacobian in SNESQN.

694:     Logically Collective on SNES

696:     Input Parameters:
697: +   snes - the iterative context
698: -   stype - scale type

700:     Options Database:
701: .   -snes_qn_scale_type <shanno,none,linesearch,jacobian>

703:     Level: intermediate

705:     SNESQNScaleTypes:
706: +   SNES_QN_SCALE_NONE - don't scale the problem
707: .   SNES_QN_SCALE_SHANNO - use shanno scaling
708: .   SNES_QN_SCALE_LINESEARCH - scale based upon line search lambda
709: -   SNES_QN_SCALE_JACOBIAN - scale by solving a linear system coming from the Jacobian you provided with SNESSetJacobian() computed at the first iteration 
710:                              of QN and at ever restart.

712: .keywords: scaling type

714: .seealso: SNES, SNESQN, SNESLineSearch, SNESQNScaleType, SNESetJacobian()
715: @*/

717: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
718: {

723:   PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));
724:   return(0);
725: }

729: PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
730: {
731:   SNES_QN *qn = (SNES_QN*)snes->data;

734:   qn->scale_type = stype;
735:   return(0);
736: }

740: PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
741: {
742:   SNES_QN *qn = (SNES_QN*)snes->data;

745:   qn->restart_type = rtype;
746:   return(0);
747: }

751: /*@
752:     SNESQNSetType - Sets the quasi-Newton variant to be used in SNESQN.

754:     Logically Collective on SNES

756:     Input Parameters:
757: +   snes - the iterative context
758: -   qtype - variant type

760:     Options Database:
761: .   -snes_qn_type <lbfgs,broyden,badbroyden>

763:     Level: beginner

765:     SNESQNTypes:
766: +   SNES_QN_LBFGS - LBFGS variant
767: .   SNES_QN_BROYDEN - Broyden variant
768: -   SNES_QN_BADBROYDEN - Bad Broyden variant

770: .keywords: SNES, SNESQN, type, set, SNESQNType
771: @*/

773: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
774: {

779:   PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));
780:   return(0);
781: }

785: PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
786: {
787:   SNES_QN *qn = (SNES_QN*)snes->data;

790:   qn->type = qtype;
791:   return(0);
792: }

794: /* -------------------------------------------------------------------------- */
795: /*MC
796:       SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.

798:       Options Database:

800: +     -snes_qn_m <m> - Number of past states saved for the L-Broyden methods.
801: +     -snes_qn_restart_type <powell,periodic,none> - set the restart type
802: .     -snes_qn_powell_gamma - Angle condition for restart.
803: .     -snes_qn_powell_descent - Descent condition for restart.
804: .     -snes_qn_type <lbfgs,broyden,badbroyden> - QN type
805: .     -snes_qn_scale_type <shanno,none,linesearch,jacobian> - scaling performed on inner Jacobian
806: .     -snes_linesearch_type <cp, l2, basic> - Type of line search.
807: -     -snes_qn_monitor - Monitors the quasi-newton Jacobian.

809:       Notes: This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of F(x) = b using
810:       previous change in F(x) and x to form the approximate inverse Jacobian using a series of multiplicative rank-one
811:       updates.

813:       When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied.  The first of
814:       these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this
815:       iteration as the current iteration's values when constructing the approximate Jacobian.  The second, composed,
816:       perturbs the problem the Jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner.

818:       References:
819: +   1. -   Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995.
820: .   2. -   R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi Newton Matrices and their use in Limited Memory Methods,
821:       Technical Report, Northwestern University, June 1992.
822: .   3. -   Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker, Experiments with Quasi-Newton Methods in Solving Stiff ODE
823:       Systems, SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985.
824: -   4. -   Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers", 
825:        SIAM Review, 57(4), 2015

827:       Level: beginner

829: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR

831: M*/
834: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
835: {
837:   SNES_QN        *qn;

840:   snes->ops->setup          = SNESSetUp_QN;
841:   snes->ops->solve          = SNESSolve_QN;
842:   snes->ops->destroy        = SNESDestroy_QN;
843:   snes->ops->setfromoptions = SNESSetFromOptions_QN;
844:   snes->ops->view           = SNESView_QN;
845:   snes->ops->reset          = SNESReset_QN;

847:   snes->pcside = PC_LEFT;

849:   snes->usespc  = PETSC_TRUE;
850:   snes->usesksp = PETSC_FALSE;

852:   if (!snes->tolerancesset) {
853:     snes->max_funcs = 30000;
854:     snes->max_its   = 10000;
855:   }

857:   PetscNewLog(snes,&qn);
858:   snes->data          = (void*) qn;
859:   qn->m               = 10;
860:   qn->scaling         = 1.0;
861:   qn->U               = NULL;
862:   qn->V               = NULL;
863:   qn->dXtdF           = NULL;
864:   qn->dFtdX           = NULL;
865:   qn->dXdFmat         = NULL;
866:   qn->monitor         = NULL;
867:   qn->singlereduction = PETSC_TRUE;
868:   qn->powell_gamma    = 0.9999;
869:   qn->scale_type      = SNES_QN_SCALE_DEFAULT;
870:   qn->restart_type    = SNES_QN_RESTART_DEFAULT;
871:   qn->type            = SNES_QN_LBFGS;

873:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);
874:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);
875:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);
876:   return(0);
877: }