Actual source code: qn.c

petsc-3.8.4 2018-03-24
Report Typos and Errors
  1:  #include <petsc/private/snesimpl.h>
  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;

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

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

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

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

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

119:   /* ksp thing for Jacobian scaling */
120:   PetscInt           h,k,j,i,lits;
121:   PetscInt           m = qn->m;
122:   PetscScalar        gdot,udot;
123:   PetscInt           l = m;

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

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

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

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

182:   /* ksp thing for Jacobian scaling */
183:   PetscInt           k,i,j,g,lits;
184:   PetscInt           m = qn->m;
185:   PetscScalar        t;
186:   PetscInt           l = m;

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

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

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

277: static PetscErrorCode SNESSolve_QN(SNES snes)
278: {
279:   PetscErrorCode       ierr;
280:   SNES_QN              *qn = (SNES_QN*) snes->data;
281:   Vec                  X,Xold;
282:   Vec                  F,W;
283:   Vec                  Y,D,Dold;
284:   PetscInt             i, i_r;
285:   PetscReal            fnorm,xnorm,ynorm,gnorm;
286:   SNESLineSearchReason lssucceed;
287:   PetscBool            powell,periodic;
288:   PetscScalar          DolddotD,DolddotDold;
289:   SNESConvergedReason  reason;

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

294:   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);

296:   PetscCitationsRegister(SNESCitation,&SNEScite);
297:   F    = snes->vec_func;                /* residual vector */
298:   Y    = snes->vec_sol_update;          /* search direction generated by J^-1D*/
299:   W    = snes->work[3];
300:   X    = snes->vec_sol;                 /* solution vector */
301:   Xold = snes->work[0];

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

307:   snes->reason = SNES_CONVERGED_ITERATING;

309:   PetscObjectSAWsTakeAccess((PetscObject)snes);
310:   snes->iter = 0;
311:   snes->norm = 0.;
312:   PetscObjectSAWsGrantAccess((PetscObject)snes);

314:   if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
315:     SNESApplyNPC(snes,X,NULL,F);
316:     SNESGetConvergedReason(snes->npc,&reason);
317:     if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
318:       snes->reason = SNES_DIVERGED_INNER;
319:       return(0);
320:     }
321:     VecNorm(F,NORM_2,&fnorm);
322:   } else {
323:     if (!snes->vec_func_init_set) {
324:       SNESComputeFunction(snes,X,F);
325:     } else snes->vec_func_init_set = PETSC_FALSE;

327:     VecNorm(F,NORM_2,&fnorm);
328:     SNESCheckFunctionNorm(snes,fnorm);
329:   }
330:   if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
331:       SNESApplyNPC(snes,X,F,D);
332:       SNESGetConvergedReason(snes->npc,&reason);
333:       if (reason < 0  && reason != SNES_DIVERGED_MAX_IT) {
334:         snes->reason = SNES_DIVERGED_INNER;
335:         return(0);
336:       }
337:   } else {
338:     VecCopy(F,D);
339:   }

341:   PetscObjectSAWsTakeAccess((PetscObject)snes);
342:   snes->norm = fnorm;
343:   PetscObjectSAWsGrantAccess((PetscObject)snes);
344:   SNESLogConvergenceHistory(snes,fnorm,0);
345:   SNESMonitor(snes,0,fnorm);

347:   /* test convergence */
348:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
349:   if (snes->reason) return(0);

351:   if (snes->npc && snes->npcside== PC_RIGHT) {
352:     PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,0,0);
353:     SNESSolve(snes->npc,snes->vec_rhs,X);
354:     PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,0,0);
355:     SNESGetConvergedReason(snes->npc,&reason);
356:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
357:       snes->reason = SNES_DIVERGED_INNER;
358:       return(0);
359:     }
360:     SNESGetNPCFunction(snes,F,&fnorm);
361:     VecCopy(F,D);
362:   }

364:   /* scale the initial update */
365:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
366:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
367:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
368:   }

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

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

420:     if (snes->npc && snes->npcside== PC_RIGHT) {
421:       PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,0,0);
422:       SNESSolve(snes->npc,snes->vec_rhs,X);
423:       PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,0,0);
424:       SNESGetConvergedReason(snes->npc,&reason);
425:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
426:         snes->reason = SNES_DIVERGED_INNER;
427:         return(0);
428:       }
429:       SNESGetNPCFunction(snes,F,&fnorm);
430:     }

432:     SNESSetIterationNumber(snes, i+1);
433:     snes->norm = fnorm;

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

500: static PetscErrorCode SNESSetUp_QN(SNES snes)
501: {
502:   SNES_QN        *qn = (SNES_QN*)snes->data;
504:   DM             dm;


508:   if (!snes->vec_sol) {
509:     SNESGetDM(snes,&dm);
510:     DMCreateGlobalVector(dm,&snes->vec_sol);
511:   }

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

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

537:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
538:     SNESSetUpMatrices(snes);
539:   }
540:   if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) {snes->functype = SNES_FUNCTION_UNPRECONDITIONED;}
541:   return(0);
542: }

544: static PetscErrorCode SNESReset_QN(SNES snes)
545: {
547:   SNES_QN        *qn;

550:   if (snes->data) {
551:     qn = (SNES_QN*)snes->data;
552:     if (qn->U) {
553:       VecDestroyVecs(qn->m, &qn->U);
554:     }
555:     if (qn->V) {
556:       VecDestroyVecs(qn->m, &qn->V);
557:     }
558:     if (qn->singlereduction) {
559:       PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);
560:     }
561:     PetscFree4(qn->alpha,qn->beta,qn->dXtdF,qn->lambda);
562:   }
563:   return(0);
564: }

566: static PetscErrorCode SNESDestroy_QN(SNES snes)
567: {

571:   SNESReset_QN(snes);
572:   PetscFree(snes->data);
573:   PetscObjectComposeFunction((PetscObject)snes,"",NULL);
574:   return(0);
575: }

577: static PetscErrorCode SNESSetFromOptions_QN(PetscOptionItems *PetscOptionsObject,SNES snes)
578: {

580:   PetscErrorCode    ierr;
581:   SNES_QN           *qn    = (SNES_QN*)snes->data;
582:   PetscBool         monflg = PETSC_FALSE,flg;
583:   SNESLineSearch    linesearch;
584:   SNESQNRestartType rtype = qn->restart_type;
585:   SNESQNScaleType   stype = qn->scale_type;
586:   SNESQNType        qtype = qn->type;

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

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

600:   PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes,(PetscEnum)qtype,(PetscEnum*)&qtype,&flg);
601:   if (flg) {SNESQNSetType(snes,qtype);}
602:   PetscOptionsTail();
603:   if (!snes->linesearch) {
604:     SNESGetLineSearch(snes, &linesearch);
605:     if (qn->type == SNES_QN_LBFGS) {
606:       SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);
607:     } else if (qn->type == SNES_QN_BROYDEN) {
608:       SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);
609:     } else {
610:       SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
611:     }
612:   }
613:   if (monflg) {
614:     qn->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
615:   }
616:   return(0);
617: }

619: static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
620: {
621:   SNES_QN        *qn    = (SNES_QN*)snes->data;
622:   PetscBool      iascii;

626:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
627:   if (iascii) {
628:     PetscViewerASCIIPrintf(viewer,"  type is %s, restart type is %s, scale type is %s\n",SNESQNTypes[qn->type],SNESQNRestartTypes[qn->restart_type],SNESQNScaleTypes[qn->scale_type]);
629:     PetscViewerASCIIPrintf(viewer,"  Stored subspace size: %D\n", qn->m);
630:     if (qn->singlereduction) {
631:       PetscViewerASCIIPrintf(viewer,"  Using the single reduction variant.\n");
632:     }
633:   }
634:   return(0);
635: }

637: /*@
638:     SNESQNSetRestartType - Sets the restart type for SNESQN.

640:     Logically Collective on SNES

642:     Input Parameters:
643: +   snes - the iterative context
644: -   rtype - restart type

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

650:     Level: intermediate

652:     SNESQNRestartTypes:
653: +   SNES_QN_RESTART_NONE - never restart
654: .   SNES_QN_RESTART_POWELL - restart based upon descent criteria
655: -   SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations

657: .keywords: SNES, SNESQN, restart, type, set SNESLineSearch, SNESQNRestartType
658: @*/
659: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
660: {

665:   PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));
666:   return(0);
667: }

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

672:     Logically Collective on SNES

674:     Input Parameters:
675: +   snes - the iterative context
676: -   stype - scale type

678:     Options Database:
679: .   -snes_qn_scale_type <shanno,none,linesearch,jacobian>

681:     Level: intermediate

683:     SNESQNScaleTypes:
684: +   SNES_QN_SCALE_NONE - don't scale the problem
685: .   SNES_QN_SCALE_SHANNO - use shanno scaling
686: .   SNES_QN_SCALE_LINESEARCH - scale based upon line search lambda
687: -   SNES_QN_SCALE_JACOBIAN - scale by solving a linear system coming from the Jacobian you provided with SNESSetJacobian() computed at the first iteration 
688:                              of QN and at ever restart.

690: .keywords: scaling type

692: .seealso: SNES, SNESQN, SNESLineSearch, SNESQNScaleType, SNESetJacobian()
693: @*/

695: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
696: {

701:   PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));
702:   return(0);
703: }

705: PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
706: {
707:   SNES_QN *qn = (SNES_QN*)snes->data;

710:   qn->scale_type = stype;
711:   return(0);
712: }

714: PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
715: {
716:   SNES_QN *qn = (SNES_QN*)snes->data;

719:   qn->restart_type = rtype;
720:   return(0);
721: }

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

726:     Logically Collective on SNES

728:     Input Parameters:
729: +   snes - the iterative context
730: -   qtype - variant type

732:     Options Database:
733: .   -snes_qn_type <lbfgs,broyden,badbroyden>

735:     Level: beginner

737:     SNESQNTypes:
738: +   SNES_QN_LBFGS - LBFGS variant
739: .   SNES_QN_BROYDEN - Broyden variant
740: -   SNES_QN_BADBROYDEN - Bad Broyden variant

742: .keywords: SNES, SNESQN, type, set, SNESQNType
743: @*/

745: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
746: {

751:   PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));
752:   return(0);
753: }

755: PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
756: {
757:   SNES_QN *qn = (SNES_QN*)snes->data;

760:   qn->type = qtype;
761:   return(0);
762: }

764: /* -------------------------------------------------------------------------- */
765: /*MC
766:       SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.

768:       Options Database:

770: +     -snes_qn_m <m> - Number of past states saved for the L-Broyden methods.
771: +     -snes_qn_restart_type <powell,periodic,none> - set the restart type
772: .     -snes_qn_powell_gamma - Angle condition for restart.
773: .     -snes_qn_powell_descent - Descent condition for restart.
774: .     -snes_qn_type <lbfgs,broyden,badbroyden> - QN type
775: .     -snes_qn_scale_type <shanno,none,linesearch,jacobian> - scaling performed on inner Jacobian
776: .     -snes_linesearch_type <cp, l2, basic> - Type of line search.
777: -     -snes_qn_monitor - Monitors the quasi-newton Jacobian.

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

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

788:       Uses left nonlinear preconditioning by default.

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

799:       Level: beginner

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

803: M*/
804: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
805: {
807:   SNES_QN        *qn;

810:   snes->ops->setup          = SNESSetUp_QN;
811:   snes->ops->solve          = SNESSolve_QN;
812:   snes->ops->destroy        = SNESDestroy_QN;
813:   snes->ops->setfromoptions = SNESSetFromOptions_QN;
814:   snes->ops->view           = SNESView_QN;
815:   snes->ops->reset          = SNESReset_QN;

817:   snes->npcside= PC_LEFT;

819:   snes->usesnpc = PETSC_TRUE;
820:   snes->usesksp = PETSC_FALSE;

822:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

824:   if (!snes->tolerancesset) {
825:     snes->max_funcs = 30000;
826:     snes->max_its   = 10000;
827:   }

829:   PetscNewLog(snes,&qn);
830:   snes->data          = (void*) qn;
831:   qn->m               = 10;
832:   qn->scaling         = 1.0;
833:   qn->U               = NULL;
834:   qn->V               = NULL;
835:   qn->dXtdF           = NULL;
836:   qn->dFtdX           = NULL;
837:   qn->dXdFmat         = NULL;
838:   qn->monitor         = NULL;
839:   qn->singlereduction = PETSC_TRUE;
840:   qn->powell_gamma    = 0.9999;
841:   qn->scale_type      = SNES_QN_SCALE_DEFAULT;
842:   qn->restart_type    = SNES_QN_RESTART_DEFAULT;
843:   qn->type            = SNES_QN_LBFGS;

845:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);
846:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);
847:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);
848:   return(0);
849: }