Actual source code: qn.c
petsc-3.8.4 2018-03-24
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: }