Actual source code: qn.c

petsc-3.4.4 2014-03-13
  1: #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/

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

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

  9: typedef enum {SNES_QN_LBFGS      = 0,
 10:               SNES_QN_BROYDEN    = 1,
 11:               SNES_QN_BADBROYDEN = 2
 12:              } SNESQNType;

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

 33: PetscErrorCode SNESQNApply_Broyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold, Vec D, Vec Dold)
 34: {
 35:   PetscErrorCode     ierr;
 36:   SNES_QN            *qn = (SNES_QN*)snes->data;
 37:   Vec                W   = snes->work[3];
 38:   Vec                *U  = qn->U;
 39:   Vec                *V  = qn->V;
 40:   KSPConvergedReason kspreason;
 41:   PetscInt           k,i,lits;
 42:   PetscInt           m = qn->m;
 43:   PetscScalar        gdot;
 44:   PetscInt           l = m;

 47:   if (it < m) l = it;
 48:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
 49:     KSPSolve(snes->ksp,D,W);
 50:     KSPGetConvergedReason(snes->ksp,&kspreason);
 51:     if (kspreason < 0) {
 52:       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
 53:         PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
 54:         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
 55:         return(0);
 56:       }
 57:     }
 58:     KSPGetIterationNumber(snes->ksp,&lits);
 59:     snes->linear_its += lits;
 60:     VecCopy(W,Y);
 61:   } else {
 62:     VecCopy(D,Y);
 63:     VecScale(Y,qn->scaling);
 64:   }

 66:   /* inward recursion starting at the first update and working forward */
 67:   if (it > 1) {
 68:     for (i = 0; i < l-1; i++) {
 69:       k    = (it+i-l)%l;
 70:       VecDot(U[k],Y,&gdot);
 71:       VecAXPY(Y,gdot,V[k]);
 72:       if (qn->monitor) {
 73:         PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
 74:         PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));
 75:         PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
 76:       }
 77:     }
 78:   }
 79:   if (it < m) l = it;

 81:   /* set W to be the last step, post-linesearch */
 82:   VecCopy(Xold,W);
 83:   VecAXPY(W,-1.0,X);

 85:   if (l > 0) {
 86:     k    = (it-1)%l;
 87:     VecCopy(W,U[k]);
 88:     VecAXPY(W,-1.0,Y);
 89:     VecDot(U[k],W,&gdot);
 90:     VecCopy(Y,V[k]);
 91:     VecScale(V[k],1.0/gdot);
 92:     VecDot(U[k],Y,&gdot);
 93:     VecAXPY(Y,gdot,V[k]);
 94:     if (qn->monitor) {
 95:       PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
 96:       PetscViewerASCIIPrintf(qn->monitor, "update: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));
 97:       PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
 98:     }
 99:   }
100:   return(0);
101: }

105: PetscErrorCode SNESQNApply_BadBroyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D,Vec Dold)
106: {
108:   SNES_QN        *qn = (SNES_QN*)snes->data;
109:   Vec            W   = snes->work[3];
110:   Vec            *U  = qn->U;
111:   Vec            *T  = qn->V;

113:   /* ksp thing for jacobian scaling */
114:   KSPConvergedReason kspreason;
115:   PetscInt           k, i, lits;
116:   PetscInt           m = qn->m;
117:   PetscScalar        gdot;
118:   PetscInt           l = m;

121:   if (it < m) l = it;
122:   VecCopy(D,Y);
123:   if (l > 0) {
124:     k    = (it-1)%l;
125:     VecCopy(Dold,U[k]);
126:     VecAXPY(U[k],-1.0,D);
127:     VecDot(U[k],U[k],&gdot);
128:     VecCopy(D,T[k]);
129:     VecScale(T[k],1.0/gdot);
130:     VecDot(U[k],Y,&gdot);
131:     VecAXPY(Y,gdot,T[k]);
132:     if (qn->monitor) {
133:       PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
134:       PetscViewerASCIIPrintf(qn->monitor, "update: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));
135:       PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
136:     }
137:   }

139:   /* inward recursion starting at the first update and working forward */
140:   if (it > 1) {
141:     for (i = 0; i < l-1; i++) {
142:       k    = (it+i-l)%l;
143:       VecDot(U[k],Y,&gdot);
144:       VecAXPY(Y,gdot,T[k]);
145:       if (qn->monitor) {
146:         PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
147:         PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));
148:         PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
149:       }
150:     }
151:   }

153:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
154:     KSPSolve(snes->ksp,Y,W);
155:     KSPGetConvergedReason(snes->ksp,&kspreason);
156:     if (kspreason < 0) {
157:       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
158:         PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
159:         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
160:         return(0);
161:       }
162:     }
163:     KSPGetIterationNumber(snes->ksp,&lits);
164:     snes->linear_its += lits;
165:     VecCopy(W,Y);
166:   } else {
167:     VecScale(Y,qn->scaling);
168:   }
169:   return(0);
170: }

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

187:   /* ksp thing for jacobian scaling */
188:   KSPConvergedReason kspreason;
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:   if (it > 0) {
197:     k    = (it - 1) % l;
198:     VecCopy(D, dF[k]);
199:     VecAXPY(dF[k], -1.0, Dold);
200:     VecCopy(X, dX[k]);
201:     VecAXPY(dX[k], -1.0, Xold);
202:     if (qn->singlereduction) {
203:       VecMDot(dF[k], l, dX, dXtdF);
204:       VecMDot(dX[k], l, dF, dFtdX);
205:       for (j = 0; j < l; j++) {
206:         H(k, j) = dFtdX[j];
207:         H(j, k) = dXtdF[j];
208:       }
209:       /* copy back over to make the computation of alpha and beta easier */
210:       for (j = 0; j < l; j++) dXtdF[j] = H(j, j);
211:     } else {
212:       VecDot(dX[k], dF[k], &dXtdF[k]);
213:     }
214:     if (qn->scale_type == SNES_QN_SCALE_SHANNO) {
215:       PetscReal dFtdF;
216:       VecDotRealPart(dF[k],dF[k],&dFtdF);
217:       qn->scaling = PetscRealPart(dXtdF[k])/dFtdF;
218:     } else if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
219:       SNESLineSearchGetLambda(snes->linesearch,&qn->scaling);
220:     }
221:   }

223:   VecCopy(D,Y);

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

251:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
252:     KSPSolve(snes->ksp,Y,W);
253:     KSPGetConvergedReason(snes->ksp,&kspreason);
254:     if (kspreason < 0) {
255:       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
256:         PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
257:         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
258:         return(0);
259:       }
260:     }
261:     KSPGetIterationNumber(snes->ksp,&lits);
262:     snes->linear_its += lits;
263:     VecCopy(W, Y);
264:   } else {
265:     VecScale(Y, qn->scaling);
266:   }
267:   if (qn->singlereduction) {
268:     VecMDot(Y,l,dF,YtdX);
269:   }
270:   /* inward recursion starting at the first update and working forward */
271:   for (i = 0; i < l; i++) {
272:     k = (it + i - l) % l;
273:     if (qn->singlereduction) {
274:       t = YtdX[k];
275:       for (j = 0; j < i; j++) {
276:         g  = (it + j - l) % l;
277:         t += (alpha[g] - beta[g])*H(k, g);
278:       }
279:       beta[k] = t / H(k, k);
280:     } else {
281:       VecDot(dF[k], Y, &t);
282:       beta[k] = t / dXtdF[k];
283:     }
284:     VecAXPY(Y, (alpha[k] - beta[k]), dX[k]);
285:     if (qn->monitor) {
286:       PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
287:       PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d alpha - beta: %14.12e\n", it, k, PetscRealPart(alpha[k] - beta[k]));
288:       PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
289:     }
290:   }
291:   return(0);
292: }

296: static PetscErrorCode SNESSolve_QN(SNES snes)
297: {
298:   PetscErrorCode      ierr;
299:   SNES_QN             *qn = (SNES_QN*) snes->data;
300:   Vec                 X,Xold;
301:   Vec                 F,B;
302:   Vec                 Y,FPC,D,Dold;
303:   SNESConvergedReason reason;
304:   PetscInt            i, i_r;
305:   PetscReal           fnorm,xnorm,ynorm,gnorm;
306:   PetscBool           lssucceed,powell,periodic;
307:   PetscScalar         DolddotD,DolddotDold,DdotD,YdotD;
308:   MatStructure        flg = DIFFERENT_NONZERO_PATTERN;

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

313:   F = snes->vec_func;                   /* residual vector */
314:   Y = snes->vec_sol_update;             /* search direction generated by J^-1D*/
315:   B = snes->vec_rhs;

317:   X    = snes->vec_sol;                 /* solution vector */
318:   Xold = snes->work[0];

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

324:   snes->reason = SNES_CONVERGED_ITERATING;

326:   PetscObjectAMSTakeAccess((PetscObject)snes);
327:   snes->iter = 0;
328:   snes->norm = 0.;
329:   PetscObjectAMSGrantAccess((PetscObject)snes);
330:   if (!snes->vec_func_init_set) {
331:     SNESComputeFunction(snes,X,F);
332:     if (snes->domainerror) {
333:       snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
334:       return(0);
335:     }
336:   } else snes->vec_func_init_set = PETSC_FALSE;

338:   if (!snes->norm_init_set) {
339:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
340:     if (PetscIsInfOrNanReal(fnorm)) {
341:       snes->reason = SNES_DIVERGED_FNORM_NAN;
342:       return(0);
343:     }
344:   } else {
345:     fnorm               = snes->norm_init;
346:     snes->norm_init_set = PETSC_FALSE;
347:   }

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

355:   /* set parameter for default relative tolerance convergence test */
356:   snes->ttol = fnorm*snes->rtol;
357:   /* test convergence */
358:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
359:   if (snes->reason) return(0);

361:   /* composed solve */
362:   if (snes->pc && snes->pcside == PC_RIGHT) {
363:     SNESSetInitialFunction(snes->pc, F);
364:     SNESSetInitialFunctionNorm(snes->pc, fnorm);

366:     PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,B,0);
367:     SNESSolve(snes->pc, B, X);
368:     PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,B,0);

370:     SNESGetConvergedReason(snes->pc,&reason);
371:     if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) {
372:       snes->reason = SNES_DIVERGED_INNER;
373:       return(0);
374:     }
375:     SNESGetFunction(snes->pc, &FPC, NULL, NULL);
376:     VecCopy(FPC, F);
377:     SNESGetFunctionNorm(snes->pc, &fnorm);
378:     VecCopy(F, Y);
379:   } else {
380:     VecCopy(F, Y);
381:   }
382:   VecCopy(Y, D);

384:   /* scale the initial update */
385:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
386:     SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
387:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);
388:   }

390:   for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
391:     switch (qn->type) {
392:     case SNES_QN_BADBROYDEN:
393:       SNESQNApply_BadBroyden(snes,i_r,Y,X,Xold,D,Dold);
394:       break;
395:     case SNES_QN_BROYDEN:
396:       SNESQNApply_Broyden(snes,i_r,Y,X,Xold,D,Dold);
397:       break;
398:     case SNES_QN_LBFGS:
399:       SNESQNApply_LBFGS(snes,i_r,Y,X,Xold,D,Dold);
400:       break;
401:     }
402:     /* line search for lambda */
403:     ynorm = 1; gnorm = fnorm;
404:     VecCopy(D, Dold);
405:     VecCopy(X, Xold);
406:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
407:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
408:     if (snes->domainerror) {
409:       snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
410:       return(0);
411:     }
412:     SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);
413:     if (!lssucceed) {
414:       if (++snes->numFailures >= snes->maxFailures) {
415:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
416:         break;
417:       }
418:     }
419:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
420:     if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
421:       SNESLineSearchGetLambda(snes->linesearch, &qn->scaling);
422:     }

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

427:     SNESSetIterationNumber(snes, i+1);
428:     SNESSetFunctionNorm(snes, fnorm);

430:     SNESLogConvergenceHistory(snes,snes->norm,snes->iter);
431:     SNESMonitor(snes,snes->iter,snes->norm);
432:     /* set parameter for default relative tolerance convergence test */
433:     (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
434:     if (snes->reason) return(0);

436:     if (snes->pc && snes->pcside == PC_RIGHT) {
437:       SNESSetInitialFunction(snes->pc, F);
438:       SNESSetInitialFunctionNorm(snes->pc, fnorm);
439:       PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,B,0);
440:       SNESSolve(snes->pc, B, X);
441:       PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,B,0);
442:       SNESGetConvergedReason(snes->pc,&reason);
443:       if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) {
444:         snes->reason = SNES_DIVERGED_INNER;
445:         return(0);
446:       }
447:       SNESGetFunction(snes->pc, &FPC, NULL, NULL);
448:       VecCopy(FPC, F);
449:       SNESGetFunctionNorm(snes->pc, &fnorm);
450:       VecCopy(F, D);
451:     } else {
452:       VecCopy(F, D);
453:     }

455:     powell = PETSC_FALSE;
456:     if (qn->restart_type == SNES_QN_RESTART_POWELL) {
457:       /* check restart by Powell's Criterion: |F^T H_0 Fold| > 0.2 * |Fold^T H_0 Fold| */
458:       VecDotBegin(Dold, Dold, &DolddotDold);
459:       VecDotBegin(Dold, D, &DolddotD);
460:       VecDotBegin(D, D, &DdotD);
461:       VecDotBegin(Y, D, &YdotD);
462:       VecDotEnd(Dold, Dold, &DolddotDold);
463:       VecDotEnd(Dold, D, &DolddotD);
464:       VecDotEnd(D, D, &DdotD);
465:       VecDotEnd(Y, D, &YdotD);
466:       if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE;
467:     }
468:     periodic = PETSC_FALSE;
469:     if (qn->restart_type == SNES_QN_RESTART_PERIODIC) {
470:       if (i_r>qn->m-1) periodic = PETSC_TRUE;
471:     }
472:     /* restart if either powell or periodic restart is satisfied. */
473:     if (powell || periodic) {
474:       if (qn->monitor) {
475:         PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
476:         PetscViewerASCIIPrintf(qn->monitor, "restart! |%14.12e| > %4.2f*|%14.12e| or i_r = %d\n", PetscRealPart(DolddotD), qn->powell_gamma, PetscRealPart(DolddotDold), i_r);
477:         PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
478:       }
479:       i_r = -1;
480:       /* general purpose update */
481:       if (snes->ops->update) {
482:         (*snes->ops->update)(snes, snes->iter);
483:       }
484:       if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
485:         SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
486:         KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);
487:       }
488:     }
489:     /* general purpose update */
490:     if (snes->ops->update) {
491:       (*snes->ops->update)(snes, snes->iter);
492:     }
493:   }
494:   if (i == snes->max_its) {
495:     PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);
496:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
497:   }
498:   return(0);
499: }

503: static PetscErrorCode SNESSetUp_QN(SNES snes)
504: {
505:   SNES_QN        *qn = (SNES_QN*)snes->data;

509:   VecDuplicateVecs(snes->vec_sol, qn->m, &qn->U);
510:   VecDuplicateVecs(snes->vec_sol, qn->m, &qn->V);
511:   PetscMalloc3(qn->m, PetscScalar, &qn->alpha,
512:                       qn->m, PetscScalar, &qn->beta,
513:                       qn->m, PetscScalar, &qn->dXtdF);

515:   if (qn->singlereduction) {
516:     PetscMalloc3(qn->m*qn->m, PetscScalar, &qn->dXdFmat,
517:                         qn->m, PetscScalar, &qn->dFtdX,
518:                         qn->m, PetscScalar, &qn->YtdX);
519:   }
520:   SNESSetWorkVecs(snes,4);

522:   /* set up the line search */
523:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
524:     SNESSetUpMatrices(snes);
525:   }
526:   return(0);
527: }

531: static PetscErrorCode SNESReset_QN(SNES snes)
532: {
534:   SNES_QN        *qn;

537:   if (snes->data) {
538:     qn = (SNES_QN*)snes->data;
539:     if (qn->U) {
540:       VecDestroyVecs(qn->m, &qn->U);
541:     }
542:     if (qn->V) {
543:       VecDestroyVecs(qn->m, &qn->V);
544:     }
545:     if (qn->singlereduction) {
546:       PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);
547:     }
548:     PetscFree3(qn->alpha, qn->beta, qn->dXtdF);
549:   }
550:   return(0);
551: }

555: static PetscErrorCode SNESDestroy_QN(SNES snes)
556: {

560:   SNESReset_QN(snes);
561:   PetscFree(snes->data);
562:   PetscObjectComposeFunction((PetscObject)snes,"",NULL);
563:   return(0);
564: }

568: static PetscErrorCode SNESSetFromOptions_QN(SNES snes)
569: {

571:   PetscErrorCode    ierr;
572:   SNES_QN           *qn    = (SNES_QN*)snes->data;
573:   PetscBool         monflg = PETSC_FALSE,flg;
574:   SNESLineSearch    linesearch;
575:   SNESQNRestartType rtype = qn->restart_type;
576:   SNESQNScaleType   stype = qn->scale_type;

579:   PetscOptionsHead("SNES QN options");
580:   PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,NULL);
581:   PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance",          "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL);
582:   PetscOptionsReal("-snes_qn_powell_downhill","Powell descent tolerance",        "SNESQN", qn->powell_downhill, &qn->powell_downhill, NULL);
583:   PetscOptionsBool("-snes_qn_monitor",         "Monitor for the QN methods",      "SNESQN", monflg, &monflg, NULL);
584:   PetscOptionsBool("-snes_qn_single_reduction", "Aggregate reductions",           "SNESQN", qn->singlereduction, &qn->singlereduction, NULL);
585:   PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);
586:   if (flg) SNESQNSetScaleType(snes,stype);

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

591:   PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes,
592:                           (PetscEnum)qn->type,(PetscEnum*)&qn->type,NULL);
593:   PetscOptionsTail();
594:   if (!snes->linesearch) {
595:     SNESGetLineSearch(snes, &linesearch);
596:     SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);
597:   }
598:   if (monflg) {
599:     qn->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
600:   }
601:   return(0);
602: }


607: /*@
608:     SNESQNSetRestartType - Sets the restart type for SNESQN.

610:     Logically Collective on SNES

612:     Input Parameters:
613: +   snes - the iterative context
614: -   rtype - restart type

616:     Options Database:
617: +   -snes_qn_restart_type<powell,periodic,none> - set the restart type
618: -   -snes_qn_restart[10] - sets the number of iterations before restart for periodic

620:     Level: intermediate

622:     SNESQNRestartTypes:
623: +   SNES_QN_RESTART_NONE - never restart
624: .   SNES_QN_RESTART_POWELL - restart based upon descent criteria
625: -   SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations

627:     Notes:
628:     The default line search used is the L2 line search and it requires two additional function evaluations.

630: .keywords: SNES, SNESQN, restart, type, set SNESLineSearch
631: @*/
632: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
633: {

638:   PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));
639:   return(0);
640: }

644: /*@
645:     SNESQNSetScaleType - Sets the scaling type for the inner inverse jacobian in SNESQN.

647:     Logically Collective on SNES

649:     Input Parameters:
650: +   snes - the iterative context
651: -   stype - scale type

653:     Options Database:
654: .   -snes_qn_scale_type<shanno,none,linesearch,jacobian>

656:     Level: intermediate

658:     SNESQNSelectTypes:
659: +   SNES_QN_SCALE_NONE - don't scale the problem
660: .   SNES_QN_SCALE_SHANNO - use shanno scaling
661: .   SNES_QN_SCALE_LINESEARCH - scale based upon line search lambda
662: -   SNES_QN_SCALE_JACOBIAN - scale by inverting a previously computed Jacobian.

664: .keywords: SNES, SNESQN, scaling, type, set SNESLineSearch
665: @*/

667: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
668: {

673:   PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));
674:   return(0);
675: }

679: PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
680: {
681:   SNES_QN *qn = (SNES_QN*)snes->data;

684:   qn->scale_type = stype;
685:   return(0);
686: }

690: PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
691: {
692:   SNES_QN *qn = (SNES_QN*)snes->data;

695:   qn->restart_type = rtype;
696:   return(0);
697: }

699: /* -------------------------------------------------------------------------- */
700: /*MC
701:       SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.

703:       Options Database:

705: +     -snes_qn_m - Number of past states saved for the L-Broyden methods.
706: .     -snes_qn_powell_angle - Angle condition for restart.
707: .     -snes_qn_powell_descent - Descent condition for restart.
708: .     -snes_linesearch_type <cp, l2, basic> - Type of line search.
709: -     -snes_qn_monitor - Monitors the quasi-newton jacobian.

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

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

720:       References:

722:       L-Broyden Methods: a generalization of the L-BFGS method to the limited memory Broyden family, M. B. Reed,
723:       International Journal of Computer Mathematics, vol. 86, 2009.

725:       Experiments with Quasi-Newton Methods in Solving Stiff ODE Systems, Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker
726:       SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985.


729:       Level: beginner

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

733: M*/
736: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
737: {
739:   SNES_QN        *qn;

742:   snes->ops->setup          = SNESSetUp_QN;
743:   snes->ops->solve          = SNESSolve_QN;
744:   snes->ops->destroy        = SNESDestroy_QN;
745:   snes->ops->setfromoptions = SNESSetFromOptions_QN;
746:   snes->ops->view           = 0;
747:   snes->ops->reset          = SNESReset_QN;

749:   snes->usespc  = PETSC_TRUE;
750:   snes->usesksp = PETSC_FALSE;

752:   if (!snes->tolerancesset) {
753:     snes->max_funcs = 30000;
754:     snes->max_its   = 10000;
755:   }

757:   PetscNewLog(snes,SNES_QN,&qn);
758:   snes->data          = (void*) qn;
759:   qn->m               = 10;
760:   qn->scaling         = 1.0;
761:   qn->U               = NULL;
762:   qn->V               = NULL;
763:   qn->dXtdF           = NULL;
764:   qn->dFtdX           = NULL;
765:   qn->dXdFmat         = NULL;
766:   qn->monitor         = NULL;
767:   qn->singlereduction = PETSC_FALSE;
768:   qn->powell_gamma    = 0.9999;
769:   qn->scale_type      = SNES_QN_SCALE_SHANNO;
770:   qn->restart_type    = SNES_QN_RESTART_POWELL;
771:   qn->type            = SNES_QN_LBFGS;

773:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);
774:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);
775:   return(0);
776: }