Actual source code: qn.c
petsc-3.7.5 2017-01-01
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: Uses left nonlinear preconditioning by default.
820: References:
821: + 1. - Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995.
822: . 2. - R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi Newton Matrices and their use in Limited Memory Methods,
823: Technical Report, Northwestern University, June 1992.
824: . 3. - Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker, Experiments with Quasi-Newton Methods in Solving Stiff ODE
825: Systems, SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985.
826: - 4. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
827: SIAM Review, 57(4), 2015
829: Level: beginner
831: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR
833: M*/
836: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
837: {
839: SNES_QN *qn;
842: snes->ops->setup = SNESSetUp_QN;
843: snes->ops->solve = SNESSolve_QN;
844: snes->ops->destroy = SNESDestroy_QN;
845: snes->ops->setfromoptions = SNESSetFromOptions_QN;
846: snes->ops->view = SNESView_QN;
847: snes->ops->reset = SNESReset_QN;
849: snes->pcside = PC_LEFT;
851: snes->usespc = PETSC_TRUE;
852: snes->usesksp = PETSC_FALSE;
854: if (!snes->tolerancesset) {
855: snes->max_funcs = 30000;
856: snes->max_its = 10000;
857: }
859: PetscNewLog(snes,&qn);
860: snes->data = (void*) qn;
861: qn->m = 10;
862: qn->scaling = 1.0;
863: qn->U = NULL;
864: qn->V = NULL;
865: qn->dXtdF = NULL;
866: qn->dFtdX = NULL;
867: qn->dXdFmat = NULL;
868: qn->monitor = NULL;
869: qn->singlereduction = PETSC_TRUE;
870: qn->powell_gamma = 0.9999;
871: qn->scale_type = SNES_QN_SCALE_DEFAULT;
872: qn->restart_type = SNES_QN_RESTART_DEFAULT;
873: qn->type = SNES_QN_LBFGS;
875: PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);
876: PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);
877: PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);
878: return(0);
879: }