Actual source code: qn.c
petsc-3.10.5 2019-03-28
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,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: VecCopy(W,Y);
58: } else {
59: VecCopy(D,Y);
60: VecScale(Y,qn->scaling);
61: }
63: /* inward recursion starting at the first update and working forward */
64: for (i = 0; i < l-1; i++) {
65: j = (it+i-l)%l;
66: k = (it+i-l+1)%l;
67: VecNorm(U[j],NORM_2,&unorm);
68: VecDot(U[j],Y,&gdot);
69: unorm *= unorm;
70: udot = PetscRealPart(gdot);
71: a = (lambda[j]/lambda[k]);
72: b = -(1.-lambda[j]);
73: a *= udot/unorm;
74: b *= udot/unorm;
75: VecAXPBYPCZ(Y,a,b,1.,U[k],U[j]);
77: if (qn->monitor) {
78: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
79: PetscViewerASCIIPrintf(qn->monitor, "using vector %D and %D, gdot: %14.12e\n",k,j,(double)PetscRealPart(gdot));
80: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
81: }
82: }
83: if (l > 0) {
84: k = (it-1)%l;
85: VecDot(U[k],Y,&gdot);
86: VecNorm(U[k],NORM_2,&unorm);
87: unorm *= unorm;
88: udot = PetscRealPart(gdot);
89: a = unorm/(unorm-lambda[k]*udot);
90: b = -(1.-lambda[k])*udot/(unorm-lambda[k]*udot);
91: if (qn->monitor) {
92: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
93: PetscViewerASCIIPrintf(qn->monitor, "using vector %D: a: %14.12e b: %14.12e \n",k,(double)a,(double)b);
94: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
95: }
96: VecAXPBY(Y,b,a,U[k]);
97: }
98: l = m;
99: if (it+1<m)l=it+1;
100: k = it%l;
101: if (qn->monitor) {
102: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
103: PetscViewerASCIIPrintf(qn->monitor, "setting vector %D of %D\n",k,l);
104: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
105: }
106: return(0);
107: }
109: PetscErrorCode SNESQNApply_BadBroyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D,Vec Dold)
110: {
112: SNES_QN *qn = (SNES_QN*)snes->data;
113: Vec W = snes->work[3];
114: Vec *U = qn->U;
115: Vec *T = qn->V;
117: /* ksp thing for Jacobian scaling */
118: PetscInt h,k,j,i;
119: PetscInt m = qn->m;
120: PetscScalar gdot,udot;
121: PetscInt l = m;
124: if (it < m) l = it;
125: VecCopy(D,Y);
126: if (l > 0) {
127: k = (it-1)%l;
128: SNESLineSearchGetLambda(snes->linesearch,&qn->lambda[k]);
129: VecCopy(Dold,U[k]);
130: VecAXPY(U[k],-1.0,D);
131: VecCopy(Xold,T[k]);
132: VecAXPY(T[k],-1.0,X);
133: }
135: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
136: KSPSolve(snes->ksp,Y,W);
137: SNESCheckKSPSolve(snes);
138: VecCopy(W,Y);
139: } else {
140: VecScale(Y,qn->scaling);
141: }
143: /* inward recursion starting at the first update and working forward */
144: if (l > 0) {
145: for (i = 0; i < l-1; i++) {
146: j = (it+i-l)%l;
147: k = (it+i-l+1)%l;
148: h = (it-1)%l;
149: VecDotBegin(U[j],U[h],&gdot);
150: VecDotBegin(U[j],U[j],&udot);
151: VecDotEnd(U[j],U[h],&gdot);
152: VecDotEnd(U[j],U[j],&udot);
153: VecAXPY(Y,PetscRealPart(gdot)/PetscRealPart(udot),T[k]);
154: VecAXPY(Y,-(1.-qn->lambda[k])*PetscRealPart(gdot)/PetscRealPart(udot),T[j]);
155: if (qn->monitor) {
156: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
157: PetscViewerASCIIPrintf(qn->monitor, "it: %D k: %D gdot: %14.12e\n", it, k, (double)PetscRealPart(gdot));
158: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
159: }
160: }
161: }
162: return(0);
163: }
165: PetscErrorCode SNESQNApply_LBFGS(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D,Vec Dold)
166: {
168: SNES_QN *qn = (SNES_QN*)snes->data;
169: Vec W = snes->work[3];
170: Vec *dX = qn->U;
171: Vec *dF = qn->V;
172: PetscScalar *alpha = qn->alpha;
173: PetscScalar *beta = qn->beta;
174: PetscScalar *dXtdF = qn->dXtdF;
175: PetscScalar *dFtdX = qn->dFtdX;
176: PetscScalar *YtdX = qn->YtdX;
178: /* ksp thing for Jacobian scaling */
179: PetscInt k,i,j,g;
180: PetscInt m = qn->m;
181: PetscScalar t;
182: PetscInt l = m;
185: if (it < m) l = it;
186: VecCopy(D,Y);
187: if (it > 0) {
188: k = (it - 1) % l;
189: VecCopy(D,dF[k]);
190: VecAXPY(dF[k], -1.0, Dold);
191: VecCopy(X, dX[k]);
192: VecAXPY(dX[k], -1.0, Xold);
193: if (qn->singlereduction) {
194: VecMDotBegin(dF[k],l,dX,dXtdF);
195: VecMDotBegin(dX[k],l,dF,dFtdX);
196: VecMDotBegin(Y,l,dX,YtdX);
197: VecMDotEnd(dF[k],l,dX,dXtdF);
198: VecMDotEnd(dX[k],l,dF,dFtdX);
199: VecMDotEnd(Y,l,dX,YtdX);
200: for (j = 0; j < l; j++) {
201: H(k, j) = dFtdX[j];
202: H(j, k) = dXtdF[j];
203: }
204: /* copy back over to make the computation of alpha and beta easier */
205: for (j = 0; j < l; j++) dXtdF[j] = H(j, j);
206: } else {
207: VecDot(dX[k], dF[k], &dXtdF[k]);
208: }
209: if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
210: SNESLineSearchGetLambda(snes->linesearch,&qn->scaling);
211: }
212: }
214: /* outward recursion starting at iteration k's update and working back */
215: for (i=0; i<l; i++) {
216: k = (it-i-1)%l;
217: if (qn->singlereduction) {
218: /* construct t = dX[k] dot Y as Y_0 dot dX[k] + sum(-alpha[j]dX[k]dF[j]) */
219: t = YtdX[k];
220: for (j=0; j<i; j++) {
221: g = (it-j-1)%l;
222: t -= alpha[g]*H(k, g);
223: }
224: alpha[k] = t/H(k,k);
225: } else {
226: VecDot(dX[k],Y,&t);
227: alpha[k] = t/dXtdF[k];
228: }
229: if (qn->monitor) {
230: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
231: PetscViewerASCIIPrintf(qn->monitor, "it: %D k: %D alpha: %14.12e\n", it, k, (double)PetscRealPart(alpha[k]));
232: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
233: }
234: VecAXPY(Y,-alpha[k],dF[k]);
235: }
237: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
238: KSPSolve(snes->ksp,Y,W);
239: SNESCheckKSPSolve(snes);
240: VecCopy(W, Y);
241: } else {
242: VecScale(Y, qn->scaling);
243: }
244: if (qn->singlereduction) {
245: VecMDot(Y,l,dF,YtdX);
246: }
247: /* inward recursion starting at the first update and working forward */
248: for (i = 0; i < l; i++) {
249: k = (it + i - l) % l;
250: if (qn->singlereduction) {
251: t = YtdX[k];
252: for (j = 0; j < i; j++) {
253: g = (it + j - l) % l;
254: t += (alpha[g] - beta[g])*H(g, k);
255: }
256: beta[k] = t / H(k, k);
257: } else {
258: VecDot(dF[k], Y, &t);
259: beta[k] = t / dXtdF[k];
260: }
261: VecAXPY(Y, (alpha[k] - beta[k]), dX[k]);
262: if (qn->monitor) {
263: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
264: PetscViewerASCIIPrintf(qn->monitor, "it: %D k: %D alpha - beta: %14.12e\n", it, k, (double)PetscRealPart(alpha[k] - beta[k]));
265: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
266: }
267: }
268: return(0);
269: }
271: static PetscErrorCode SNESSolve_QN(SNES snes)
272: {
273: PetscErrorCode ierr;
274: SNES_QN *qn = (SNES_QN*) snes->data;
275: Vec X,Xold;
276: Vec F,W;
277: Vec Y,D,Dold;
278: PetscInt i, i_r;
279: PetscReal fnorm,xnorm,ynorm,gnorm;
280: SNESLineSearchReason lssucceed;
281: PetscBool powell,periodic;
282: PetscScalar DolddotD,DolddotDold;
283: SNESConvergedReason reason;
285: /* basically just a regular newton's method except for the application of the Jacobian */
288: 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);
290: PetscCitationsRegister(SNESCitation,&SNEScite);
291: F = snes->vec_func; /* residual vector */
292: Y = snes->vec_sol_update; /* search direction generated by J^-1D*/
293: W = snes->work[3];
294: X = snes->vec_sol; /* solution vector */
295: Xold = snes->work[0];
297: /* directions generated by the preconditioned problem with F_pre = F or x - M(x, b) */
298: D = snes->work[1];
299: Dold = snes->work[2];
301: snes->reason = SNES_CONVERGED_ITERATING;
303: PetscObjectSAWsTakeAccess((PetscObject)snes);
304: snes->iter = 0;
305: snes->norm = 0.;
306: PetscObjectSAWsGrantAccess((PetscObject)snes);
308: if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
309: SNESApplyNPC(snes,X,NULL,F);
310: SNESGetConvergedReason(snes->npc,&reason);
311: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
312: snes->reason = SNES_DIVERGED_INNER;
313: return(0);
314: }
315: VecNorm(F,NORM_2,&fnorm);
316: } else {
317: if (!snes->vec_func_init_set) {
318: SNESComputeFunction(snes,X,F);
319: } else snes->vec_func_init_set = PETSC_FALSE;
321: VecNorm(F,NORM_2,&fnorm);
322: SNESCheckFunctionNorm(snes,fnorm);
323: }
324: if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
325: SNESApplyNPC(snes,X,F,D);
326: SNESGetConvergedReason(snes->npc,&reason);
327: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
328: snes->reason = SNES_DIVERGED_INNER;
329: return(0);
330: }
331: } else {
332: VecCopy(F,D);
333: }
335: PetscObjectSAWsTakeAccess((PetscObject)snes);
336: snes->norm = fnorm;
337: PetscObjectSAWsGrantAccess((PetscObject)snes);
338: SNESLogConvergenceHistory(snes,fnorm,0);
339: SNESMonitor(snes,0,fnorm);
341: /* test convergence */
342: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
343: if (snes->reason) return(0);
345: if (snes->npc && snes->npcside== PC_RIGHT) {
346: PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,0,0);
347: SNESSolve(snes->npc,snes->vec_rhs,X);
348: PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,0,0);
349: SNESGetConvergedReason(snes->npc,&reason);
350: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
351: snes->reason = SNES_DIVERGED_INNER;
352: return(0);
353: }
354: SNESGetNPCFunction(snes,F,&fnorm);
355: VecCopy(F,D);
356: }
358: /* scale the initial update */
359: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
360: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
361: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
362: }
364: for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
365: if (qn->scale_type == SNES_QN_SCALE_SHANNO && i_r > 0) {
366: PetscScalar ff,xf;
367: VecCopy(Dold,Y);
368: VecCopy(Xold,W);
369: VecAXPY(Y,-1.0,D);
370: VecAXPY(W,-1.0,X);
371: VecDotBegin(Y,Y,&ff);
372: VecDotBegin(W,Y,&xf);
373: VecDotEnd(Y,Y,&ff);
374: VecDotEnd(W,Y,&xf);
375: qn->scaling = PetscRealPart(xf)/PetscRealPart(ff);
376: if (qn->monitor) {
377: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
378: PetscViewerASCIIPrintf(qn->monitor, "Shanno scaling %D %g\n", i,(double)qn->scaling);
379: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
380: }
381: }
382: switch (qn->type) {
383: case SNES_QN_BADBROYDEN:
384: SNESQNApply_BadBroyden(snes,i_r,Y,X,Xold,D,Dold);
385: break;
386: case SNES_QN_BROYDEN:
387: SNESQNApply_Broyden(snes,i_r,Y,X,Xold,D);
388: break;
389: case SNES_QN_LBFGS:
390: SNESQNApply_LBFGS(snes,i_r,Y,X,Xold,D,Dold);
391: break;
392: }
393: /* line search for lambda */
394: ynorm = 1; gnorm = fnorm;
395: VecCopy(D, Dold);
396: VecCopy(X, Xold);
397: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
398: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
399: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
400: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
401: if (lssucceed) {
402: if (++snes->numFailures >= snes->maxFailures) {
403: snes->reason = SNES_DIVERGED_LINE_SEARCH;
404: break;
405: }
406: }
407: if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
408: SNESLineSearchGetLambda(snes->linesearch, &qn->scaling);
409: }
411: /* convergence monitoring */
412: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);
414: if (snes->npc && snes->npcside== PC_RIGHT) {
415: PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,0,0);
416: SNESSolve(snes->npc,snes->vec_rhs,X);
417: PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,0,0);
418: SNESGetConvergedReason(snes->npc,&reason);
419: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
420: snes->reason = SNES_DIVERGED_INNER;
421: return(0);
422: }
423: SNESGetNPCFunction(snes,F,&fnorm);
424: }
426: SNESSetIterationNumber(snes, i+1);
427: snes->norm = fnorm;
429: SNESLogConvergenceHistory(snes,snes->norm,snes->iter);
430: SNESMonitor(snes,snes->iter,snes->norm);
431: /* set parameter for default relative tolerance convergence test */
432: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
433: if (snes->reason) return(0);
434: if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
435: SNESApplyNPC(snes,X,F,D);
436: SNESGetConvergedReason(snes->npc,&reason);
437: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
438: snes->reason = SNES_DIVERGED_INNER;
439: return(0);
440: }
441: } else {
442: VecCopy(F, D);
443: }
444: powell = PETSC_FALSE;
445: if (qn->restart_type == SNES_QN_RESTART_POWELL && i_r > 1) {
446: /* check restart by Powell's Criterion: |F^T H_0 Fold| > powell_gamma * |Fold^T H_0 Fold| */
447: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
448: MatMult(snes->jacobian_pre,Dold,W);
449: } else {
450: VecCopy(Dold,W);
451: }
452: VecDotBegin(W, Dold, &DolddotDold);
453: VecDotBegin(W, D, &DolddotD);
454: VecDotEnd(W, Dold, &DolddotDold);
455: VecDotEnd(W, D, &DolddotD);
456: if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE;
457: }
458: periodic = PETSC_FALSE;
459: if (qn->restart_type == SNES_QN_RESTART_PERIODIC) {
460: if (i_r>qn->m-1) periodic = PETSC_TRUE;
461: }
462: /* restart if either powell or periodic restart is satisfied. */
463: if (powell || periodic) {
464: if (qn->monitor) {
465: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
466: if (powell) {
467: 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);
468: } else {
469: PetscViewerASCIIPrintf(qn->monitor, "Periodic restart! i_r = %D\n", i_r);
470: }
471: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
472: }
473: i_r = -1;
474: /* general purpose update */
475: if (snes->ops->update) {
476: (*snes->ops->update)(snes, snes->iter);
477: }
478: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
479: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
480: }
481: }
482: /* general purpose update */
483: if (snes->ops->update) {
484: (*snes->ops->update)(snes, snes->iter);
485: }
486: }
487: if (i == snes->max_its) {
488: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);
489: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
490: }
491: return(0);
492: }
494: static PetscErrorCode SNESSetUp_QN(SNES snes)
495: {
496: SNES_QN *qn = (SNES_QN*)snes->data;
498: DM dm;
502: if (!snes->vec_sol) {
503: SNESGetDM(snes,&dm);
504: DMCreateGlobalVector(dm,&snes->vec_sol);
505: }
507: VecDuplicateVecs(snes->vec_sol, qn->m, &qn->U);
508: if (qn->type != SNES_QN_BROYDEN) VecDuplicateVecs(snes->vec_sol, qn->m, &qn->V);
509: PetscMalloc4(qn->m,&qn->alpha,qn->m,&qn->beta,qn->m,&qn->dXtdF,qn->m,&qn->lambda);
511: if (qn->singlereduction) {
512: PetscMalloc3(qn->m*qn->m,&qn->dXdFmat,qn->m,&qn->dFtdX,qn->m,&qn->YtdX);
513: }
514: SNESSetWorkVecs(snes,4);
515: /* set method defaults */
516: if (qn->scale_type == SNES_QN_SCALE_DEFAULT) {
517: if (qn->type == SNES_QN_BADBROYDEN) {
518: qn->scale_type = SNES_QN_SCALE_NONE;
519: } else {
520: qn->scale_type = SNES_QN_SCALE_SHANNO;
521: }
522: }
523: if (qn->restart_type == SNES_QN_RESTART_DEFAULT) {
524: if (qn->type == SNES_QN_LBFGS) {
525: qn->restart_type = SNES_QN_RESTART_POWELL;
526: } else {
527: qn->restart_type = SNES_QN_RESTART_PERIODIC;
528: }
529: }
531: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
532: SNESSetUpMatrices(snes);
533: }
534: if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) {snes->functype = SNES_FUNCTION_UNPRECONDITIONED;}
535: return(0);
536: }
538: static PetscErrorCode SNESReset_QN(SNES snes)
539: {
541: SNES_QN *qn;
544: if (snes->data) {
545: qn = (SNES_QN*)snes->data;
546: if (qn->U) {
547: VecDestroyVecs(qn->m, &qn->U);
548: }
549: if (qn->V) {
550: VecDestroyVecs(qn->m, &qn->V);
551: }
552: if (qn->singlereduction) {
553: PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);
554: }
555: PetscFree4(qn->alpha,qn->beta,qn->dXtdF,qn->lambda);
556: }
557: return(0);
558: }
560: static PetscErrorCode SNESDestroy_QN(SNES snes)
561: {
565: SNESReset_QN(snes);
566: PetscFree(snes->data);
567: PetscObjectComposeFunction((PetscObject)snes,"",NULL);
568: return(0);
569: }
571: static PetscErrorCode SNESSetFromOptions_QN(PetscOptionItems *PetscOptionsObject,SNES snes)
572: {
574: PetscErrorCode ierr;
575: SNES_QN *qn = (SNES_QN*)snes->data;
576: PetscBool monflg = PETSC_FALSE,flg;
577: SNESLineSearch linesearch;
578: SNESQNRestartType rtype = qn->restart_type;
579: SNESQNScaleType stype = qn->scale_type;
580: SNESQNType qtype = qn->type;
583: PetscOptionsHead(PetscOptionsObject,"SNES QN options");
584: PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,NULL);
585: PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL);
586: PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", monflg, &monflg, NULL);
587: PetscOptionsBool("-snes_qn_single_reduction", "Aggregate reductions", "SNESQN", qn->singlereduction, &qn->singlereduction, NULL);
588: PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);
589: if (flg) SNESQNSetScaleType(snes,stype);
591: PetscOptionsEnum("-snes_qn_restart_type","Restart type","SNESQNSetRestartType",SNESQNRestartTypes,(PetscEnum)rtype,(PetscEnum*)&rtype,&flg);
592: if (flg) SNESQNSetRestartType(snes,rtype);
594: PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes,(PetscEnum)qtype,(PetscEnum*)&qtype,&flg);
595: if (flg) {SNESQNSetType(snes,qtype);}
596: PetscOptionsTail();
597: if (!snes->linesearch) {
598: SNESGetLineSearch(snes, &linesearch);
599: if (qn->type == SNES_QN_LBFGS) {
600: SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);
601: } else if (qn->type == SNES_QN_BROYDEN) {
602: SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);
603: } else {
604: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
605: }
606: }
607: if (monflg) {
608: qn->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
609: }
610: return(0);
611: }
613: static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
614: {
615: SNES_QN *qn = (SNES_QN*)snes->data;
616: PetscBool iascii;
620: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
621: if (iascii) {
622: 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]);
623: PetscViewerASCIIPrintf(viewer," Stored subspace size: %D\n", qn->m);
624: if (qn->singlereduction) {
625: PetscViewerASCIIPrintf(viewer," Using the single reduction variant.\n");
626: }
627: }
628: return(0);
629: }
631: /*@
632: SNESQNSetRestartType - Sets the restart type for SNESQN.
634: Logically Collective on SNES
636: Input Parameters:
637: + snes - the iterative context
638: - rtype - restart type
640: Options Database:
641: + -snes_qn_restart_type <powell,periodic,none> - set the restart type
642: - -snes_qn_m <m> - sets the number of stored updates and the restart period for periodic
644: Level: intermediate
646: SNESQNRestartTypes:
647: + SNES_QN_RESTART_NONE - never restart
648: . SNES_QN_RESTART_POWELL - restart based upon descent criteria
649: - SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations
651: .keywords: SNES, SNESQN, restart, type, set SNESLineSearch, SNESQNRestartType
652: @*/
653: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
654: {
659: PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));
660: return(0);
661: }
663: /*@
664: SNESQNSetScaleType - Sets the scaling type for the inner inverse Jacobian in SNESQN.
666: Logically Collective on SNES
668: Input Parameters:
669: + snes - the iterative context
670: - stype - scale type
672: Options Database:
673: . -snes_qn_scale_type <shanno,none,linesearch,jacobian>
675: Level: intermediate
677: SNESQNScaleTypes:
678: + SNES_QN_SCALE_NONE - don't scale the problem
679: . SNES_QN_SCALE_SHANNO - use shanno scaling
680: . SNES_QN_SCALE_LINESEARCH - scale based upon line search lambda
681: - SNES_QN_SCALE_JACOBIAN - scale by solving a linear system coming from the Jacobian you provided with SNESSetJacobian() computed at the first iteration
682: of QN and at ever restart.
684: .keywords: scaling type
686: .seealso: SNES, SNESQN, SNESLineSearch, SNESQNScaleType, SNESetJacobian()
687: @*/
689: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
690: {
695: PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));
696: return(0);
697: }
699: PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
700: {
701: SNES_QN *qn = (SNES_QN*)snes->data;
704: qn->scale_type = stype;
705: return(0);
706: }
708: PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
709: {
710: SNES_QN *qn = (SNES_QN*)snes->data;
713: qn->restart_type = rtype;
714: return(0);
715: }
717: /*@
718: SNESQNSetType - Sets the quasi-Newton variant to be used in SNESQN.
720: Logically Collective on SNES
722: Input Parameters:
723: + snes - the iterative context
724: - qtype - variant type
726: Options Database:
727: . -snes_qn_type <lbfgs,broyden,badbroyden>
729: Level: beginner
731: SNESQNTypes:
732: + SNES_QN_LBFGS - LBFGS variant
733: . SNES_QN_BROYDEN - Broyden variant
734: - SNES_QN_BADBROYDEN - Bad Broyden variant
736: .keywords: SNES, SNESQN, type, set, SNESQNType
737: @*/
739: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
740: {
745: PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));
746: return(0);
747: }
749: PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
750: {
751: SNES_QN *qn = (SNES_QN*)snes->data;
754: qn->type = qtype;
755: return(0);
756: }
758: /* -------------------------------------------------------------------------- */
759: /*MC
760: SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.
762: Options Database:
764: + -snes_qn_m <m> - Number of past states saved for the L-Broyden methods.
765: + -snes_qn_restart_type <powell,periodic,none> - set the restart type
766: . -snes_qn_powell_gamma - Angle condition for restart.
767: . -snes_qn_powell_descent - Descent condition for restart.
768: . -snes_qn_type <lbfgs,broyden,badbroyden> - QN type
769: . -snes_qn_scale_type <shanno,none,linesearch,jacobian> - scaling performed on inner Jacobian
770: . -snes_linesearch_type <cp, l2, basic> - Type of line search.
771: - -snes_qn_monitor - Monitors the quasi-newton Jacobian.
773: Notes:
774: This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of F(x) = b using
775: previous change in F(x) and x to form the approximate inverse Jacobian using a series of multiplicative rank-one
776: updates.
778: When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of
779: these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this
780: iteration as the current iteration's values when constructing the approximate Jacobian. The second, composed,
781: perturbs the problem the Jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner.
783: Uses left nonlinear preconditioning by default.
785: References:
786: + 1. - Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995.
787: . 2. - R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi Newton Matrices and their use in Limited Memory Methods,
788: Technical Report, Northwestern University, June 1992.
789: . 3. - Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker, Experiments with Quasi-Newton Methods in Solving Stiff ODE
790: Systems, SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985.
791: - 4. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
792: SIAM Review, 57(4), 2015
794: Level: beginner
796: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR
798: M*/
799: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
800: {
802: SNES_QN *qn;
805: snes->ops->setup = SNESSetUp_QN;
806: snes->ops->solve = SNESSolve_QN;
807: snes->ops->destroy = SNESDestroy_QN;
808: snes->ops->setfromoptions = SNESSetFromOptions_QN;
809: snes->ops->view = SNESView_QN;
810: snes->ops->reset = SNESReset_QN;
812: snes->npcside= PC_LEFT;
814: snes->usesnpc = PETSC_TRUE;
815: snes->usesksp = PETSC_FALSE;
817: snes->alwayscomputesfinalresidual = PETSC_TRUE;
819: if (!snes->tolerancesset) {
820: snes->max_funcs = 30000;
821: snes->max_its = 10000;
822: }
824: PetscNewLog(snes,&qn);
825: snes->data = (void*) qn;
826: qn->m = 10;
827: qn->scaling = 1.0;
828: qn->U = NULL;
829: qn->V = NULL;
830: qn->dXtdF = NULL;
831: qn->dFtdX = NULL;
832: qn->dXdFmat = NULL;
833: qn->monitor = NULL;
834: qn->singlereduction = PETSC_TRUE;
835: qn->powell_gamma = 0.9999;
836: qn->scale_type = SNES_QN_SCALE_DEFAULT;
837: qn->restart_type = SNES_QN_RESTART_DEFAULT;
838: qn->type = SNES_QN_LBFGS;
840: PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);
841: PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);
842: PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);
843: return(0);
844: }