Actual source code: qn.c

petsc-3.6.4 2016-04-12
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         powell_downhill;      /* Powell descent restart condition */
 23:   PetscReal         scaling;              /* scaling of H0 */
 24:   SNESQNType        type;                 /* the type of quasi-newton method used */
 25:   SNESQNScaleType   scale_type;           /* the type of scaling used */
 26:   SNESQNRestartType restart_type;         /* determine the frequency and type of restart conditions */
 27: } SNES_QN;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

319:   snes->reason = SNES_CONVERGED_ITERATING;

321:   PetscObjectSAWsTakeAccess((PetscObject)snes);
322:   snes->iter = 0;
323:   snes->norm = 0.;
324:   PetscObjectSAWsGrantAccess((PetscObject)snes);

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

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

353:   PetscObjectSAWsTakeAccess((PetscObject)snes);
354:   snes->norm = fnorm;
355:   PetscObjectSAWsGrantAccess((PetscObject)snes);
356:   SNESLogConvergenceHistory(snes,fnorm,0);
357:   SNESMonitor(snes,0,fnorm);

359:   /* test convergence */
360:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
361:   if (snes->reason) return(0);

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

376:   /* scale the initial update */
377:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
378:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
379:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
380:   }

382:   for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
383:     if (qn->scale_type == SNES_QN_SCALE_SHANNO && i_r > 0) {
384:       PetscScalar ff,xf;
385:       VecCopy(Dold,Y);
386:       VecCopy(Xold,W);
387:       VecAXPY(Y,-1.0,D);
388:       VecAXPY(W,-1.0,X);
389:       VecDotBegin(Y,Y,&ff);
390:       VecDotBegin(W,Y,&xf);
391:       VecDotEnd(Y,Y,&ff);
392:       VecDotEnd(W,Y,&xf);
393:       qn->scaling = PetscRealPart(xf)/PetscRealPart(ff);
394:     }
395:     switch (qn->type) {
396:     case SNES_QN_BADBROYDEN:
397:       SNESQNApply_BadBroyden(snes,i_r,Y,X,Xold,D,Dold);
398:       break;
399:     case SNES_QN_BROYDEN:
400:       SNESQNApply_Broyden(snes,i_r,Y,X,Xold,D);
401:       break;
402:     case SNES_QN_LBFGS:
403:       SNESQNApply_LBFGS(snes,i_r,Y,X,Xold,D,Dold);
404:       break;
405:     }
406:     /* line search for lambda */
407:     ynorm = 1; gnorm = fnorm;
408:     VecCopy(D, Dold);
409:     VecCopy(X, Xold);
410:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
411:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
412:     SNESLineSearchGetReason(snes->linesearch, &lssucceed);
413:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
414:     if (lssucceed) {
415:       if (++snes->numFailures >= snes->maxFailures) {
416:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
417:         break;
418:       }
419:     }
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:     if (snes->pc && snes->pcside == PC_RIGHT) {
428:       PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,0,0);
429:       SNESSolve(snes->pc,snes->vec_rhs,X);
430:       PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,0,0);
431:       SNESGetConvergedReason(snes->pc,&reason);
432:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
433:         snes->reason = SNES_DIVERGED_INNER;
434:         return(0);
435:       }
436:       SNESGetNPCFunction(snes,F,&fnorm);
437:     }

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

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

505: static PetscErrorCode SNESSetUp_QN(SNES snes)
506: {
507:   SNES_QN        *qn = (SNES_QN*)snes->data;
509:   DM             dm;


513:   if (!snes->vec_sol) {
514:     SNESGetDM(snes,&dm);
515:     DMCreateGlobalVector(dm,&snes->vec_sol);
516:   }

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

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

542:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
543:     SNESSetUpMatrices(snes);
544:   }
545:   if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) {snes->functype = SNES_FUNCTION_UNPRECONDITIONED;}
546:   return(0);
547: }

551: static PetscErrorCode SNESReset_QN(SNES snes)
552: {
554:   SNES_QN        *qn;

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

575: static PetscErrorCode SNESDestroy_QN(SNES snes)
576: {

580:   SNESReset_QN(snes);
581:   PetscFree(snes->data);
582:   PetscObjectComposeFunction((PetscObject)snes,"",NULL);
583:   return(0);
584: }

588: static PetscErrorCode SNESSetFromOptions_QN(PetscOptions *PetscOptionsObject,SNES snes)
589: {

591:   PetscErrorCode    ierr;
592:   SNES_QN           *qn    = (SNES_QN*)snes->data;
593:   PetscBool         monflg = PETSC_FALSE,flg;
594:   SNESLineSearch    linesearch;
595:   SNESQNRestartType rtype = qn->restart_type;
596:   SNESQNScaleType   stype = qn->scale_type;
597:   SNESQNType        qtype = qn->type;

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

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

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

633: static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
634: {
635:   SNES_QN        *qn    = (SNES_QN*)snes->data;
636:   PetscBool      iascii;

640:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
641:   if (iascii) {
642:     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]);
643:     PetscViewerASCIIPrintf(viewer,"  Stored subspace size: %d\n", qn->m);
644:     if (qn->singlereduction) {
645:       PetscViewerASCIIPrintf(viewer,"  Using the single reduction variant.\n");
646:     }
647:   }
648:   return(0);
649: }

653: /*@
654:     SNESQNSetRestartType - Sets the restart type for SNESQN.

656:     Logically Collective on SNES

658:     Input Parameters:
659: +   snes - the iterative context
660: -   rtype - restart type

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

666:     Level: intermediate

668:     SNESQNRestartTypes:
669: +   SNES_QN_RESTART_NONE - never restart
670: .   SNES_QN_RESTART_POWELL - restart based upon descent criteria
671: -   SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations

673: .keywords: SNES, SNESQN, restart, type, set SNESLineSearch, SNESQNRestartType
674: @*/
675: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
676: {

681:   PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));
682:   return(0);
683: }

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

690:     Logically Collective on SNES

692:     Input Parameters:
693: +   snes - the iterative context
694: -   stype - scale type

696:     Options Database:
697: .   -snes_qn_scale_type <shanno,none,linesearch,jacobian>

699:     Level: intermediate

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

708: .keywords: scaling type

710: .seealso: SNES, SNESQN, SNESLineSearch, SNESQNScaleType, SNESetJacobian()
711: @*/

713: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
714: {

719:   PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));
720:   return(0);
721: }

725: PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
726: {
727:   SNES_QN *qn = (SNES_QN*)snes->data;

730:   qn->scale_type = stype;
731:   return(0);
732: }

736: PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
737: {
738:   SNES_QN *qn = (SNES_QN*)snes->data;

741:   qn->restart_type = rtype;
742:   return(0);
743: }

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

750:     Logically Collective on SNES

752:     Input Parameters:
753: +   snes - the iterative context
754: -   qtype - variant type

756:     Options Database:
757: .   -snes_qn_type <lbfgs,broyden,badbroyden>

759:     Level: beginner

761:     SNESQNTypes:
762: +   SNES_QN_LBFGS - LBFGS variant
763: .   SNES_QN_BROYDEN - Broyden variant
764: -   SNES_QN_BADBROYDEN - Bad Broyden variant

766: .keywords: SNES, SNESQN, type, set, SNESQNType
767: @*/

769: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
770: {

775:   PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));
776:   return(0);
777: }

781: PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
782: {
783:   SNES_QN *qn = (SNES_QN*)snes->data;

786:   qn->type = qtype;
787:   return(0);
788: }

790: /* -------------------------------------------------------------------------- */
791: /*MC
792:       SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.

794:       Options Database:

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

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

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

814:       References:

816:       Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995.

818:       R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi-Newton Matrices and their use in Limited Memory Methods,
819:       Technical Report, Northwestern University, June 1992.

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


825:       Level: beginner

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

829: M*/
832: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
833: {
835:   SNES_QN        *qn;

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

845:   snes->pcside = PC_LEFT;

847:   snes->usespc  = PETSC_TRUE;
848:   snes->usesksp = PETSC_FALSE;

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

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

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