Actual source code: qn.c
petsc-3.4.5 2014-06-29
1: #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/
3: #define H(i,j) qn->dXdFmat[i*qn->m + j]
5: const char *const SNESQNScaleTypes[] = {"NONE","SHANNO","LINESEARCH","JACOBIAN","SNESQNScaleType","SNES_QN_SCALING_",0};
6: const char *const SNESQNRestartTypes[] = {"NONE","POWELL","PERIODIC","SNESQNRestartType","SNES_QN_RESTART_",0};
7: const char *const SNESQNTypes[] = {"LBFGS","BROYDEN","BADBROYDEN","SNESQNType","SNES_QN_RESTART_",0};
9: typedef enum {SNES_QN_LBFGS = 0,
10: SNES_QN_BROYDEN = 1,
11: SNES_QN_BADBROYDEN = 2
12: } SNESQNType;
14: typedef struct {
15: Vec *U; /* Stored past states (vary from method to method) */
16: Vec *V; /* Stored past states (vary from method to method) */
17: PetscInt m; /* The number of kept previous steps */
18: PetscScalar *alpha, *beta;
19: PetscScalar *dXtdF, *dFtdX, *YtdX;
20: PetscBool singlereduction; /* Aggregated reduction implementation */
21: PetscScalar *dXdFmat; /* A matrix of values for dX_i dot dF_j */
22: PetscViewer monitor;
23: PetscReal powell_gamma; /* Powell angle restart condition */
24: PetscReal powell_downhill; /* Powell descent restart condition */
25: PetscReal scaling; /* scaling of H0 */
26: SNESQNType type; /* the type of quasi-newton method used */
27: SNESQNScaleType scale_type; /* the type of scaling used */
28: SNESQNRestartType restart_type; /* determine the frequency and type of restart conditions */
29: } SNES_QN;
33: PetscErrorCode SNESQNApply_Broyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold, Vec D, Vec Dold)
34: {
35: PetscErrorCode ierr;
36: SNES_QN *qn = (SNES_QN*)snes->data;
37: Vec W = snes->work[3];
38: Vec *U = qn->U;
39: Vec *V = qn->V;
40: KSPConvergedReason kspreason;
41: PetscInt k,i,lits;
42: PetscInt m = qn->m;
43: PetscScalar gdot;
44: PetscInt l = m;
47: if (it < m) l = it;
48: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
49: KSPSolve(snes->ksp,D,W);
50: KSPGetConvergedReason(snes->ksp,&kspreason);
51: if (kspreason < 0) {
52: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
53: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
54: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
55: return(0);
56: }
57: }
58: KSPGetIterationNumber(snes->ksp,&lits);
59: snes->linear_its += lits;
60: VecCopy(W,Y);
61: } else {
62: VecCopy(D,Y);
63: VecScale(Y,qn->scaling);
64: }
66: /* inward recursion starting at the first update and working forward */
67: if (it > 1) {
68: for (i = 0; i < l-1; i++) {
69: k = (it+i-l)%l;
70: VecDot(U[k],Y,&gdot);
71: VecAXPY(Y,gdot,V[k]);
72: if (qn->monitor) {
73: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
74: PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));
75: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
76: }
77: }
78: }
79: if (it < m) l = it;
81: /* set W to be the last step, post-linesearch */
82: VecCopy(Xold,W);
83: VecAXPY(W,-1.0,X);
85: if (l > 0) {
86: k = (it-1)%l;
87: VecCopy(W,U[k]);
88: VecAXPY(W,-1.0,Y);
89: VecDot(U[k],W,&gdot);
90: VecCopy(Y,V[k]);
91: VecScale(V[k],1.0/gdot);
92: VecDot(U[k],Y,&gdot);
93: VecAXPY(Y,gdot,V[k]);
94: if (qn->monitor) {
95: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
96: PetscViewerASCIIPrintf(qn->monitor, "update: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));
97: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
98: }
99: }
100: return(0);
101: }
105: PetscErrorCode SNESQNApply_BadBroyden(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D,Vec Dold)
106: {
108: SNES_QN *qn = (SNES_QN*)snes->data;
109: Vec W = snes->work[3];
110: Vec *U = qn->U;
111: Vec *T = qn->V;
113: /* ksp thing for jacobian scaling */
114: KSPConvergedReason kspreason;
115: PetscInt k, i, lits;
116: PetscInt m = qn->m;
117: PetscScalar gdot;
118: PetscInt l = m;
121: if (it < m) l = it;
122: VecCopy(D,Y);
123: if (l > 0) {
124: k = (it-1)%l;
125: VecCopy(Dold,U[k]);
126: VecAXPY(U[k],-1.0,D);
127: VecDot(U[k],U[k],&gdot);
128: VecCopy(D,T[k]);
129: VecScale(T[k],1.0/gdot);
130: VecDot(U[k],Y,&gdot);
131: VecAXPY(Y,gdot,T[k]);
132: if (qn->monitor) {
133: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
134: PetscViewerASCIIPrintf(qn->monitor, "update: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));
135: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
136: }
137: }
139: /* inward recursion starting at the first update and working forward */
140: if (it > 1) {
141: for (i = 0; i < l-1; i++) {
142: k = (it+i-l)%l;
143: VecDot(U[k],Y,&gdot);
144: VecAXPY(Y,gdot,T[k]);
145: if (qn->monitor) {
146: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
147: PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d gdot: %14.12e\n", it, k, PetscRealPart(gdot));
148: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
149: }
150: }
151: }
153: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
154: KSPSolve(snes->ksp,Y,W);
155: KSPGetConvergedReason(snes->ksp,&kspreason);
156: if (kspreason < 0) {
157: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
158: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
159: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
160: return(0);
161: }
162: }
163: KSPGetIterationNumber(snes->ksp,&lits);
164: snes->linear_its += lits;
165: VecCopy(W,Y);
166: } else {
167: VecScale(Y,qn->scaling);
168: }
169: return(0);
170: }
174: PetscErrorCode SNESQNApply_LBFGS(SNES snes,PetscInt it,Vec Y,Vec X,Vec Xold,Vec D,Vec Dold)
175: {
177: SNES_QN *qn = (SNES_QN*)snes->data;
178: Vec W = snes->work[3];
179: Vec *dX = qn->U;
180: Vec *dF = qn->V;
181: PetscScalar *alpha = qn->alpha;
182: PetscScalar *beta = qn->beta;
183: PetscScalar *dXtdF = qn->dXtdF;
184: PetscScalar *dFtdX = qn->dFtdX;
185: PetscScalar *YtdX = qn->YtdX;
187: /* ksp thing for jacobian scaling */
188: KSPConvergedReason kspreason;
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: if (it > 0) {
197: k = (it - 1) % l;
198: VecCopy(D, dF[k]);
199: VecAXPY(dF[k], -1.0, Dold);
200: VecCopy(X, dX[k]);
201: VecAXPY(dX[k], -1.0, Xold);
202: if (qn->singlereduction) {
203: VecMDot(dF[k], l, dX, dXtdF);
204: VecMDot(dX[k], l, dF, dFtdX);
205: for (j = 0; j < l; j++) {
206: H(k, j) = dFtdX[j];
207: H(j, k) = dXtdF[j];
208: }
209: /* copy back over to make the computation of alpha and beta easier */
210: for (j = 0; j < l; j++) dXtdF[j] = H(j, j);
211: } else {
212: VecDot(dX[k], dF[k], &dXtdF[k]);
213: }
214: if (qn->scale_type == SNES_QN_SCALE_SHANNO) {
215: PetscReal dFtdF;
216: VecDotRealPart(dF[k],dF[k],&dFtdF);
217: qn->scaling = PetscRealPart(dXtdF[k])/dFtdF;
218: } else if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
219: SNESLineSearchGetLambda(snes->linesearch,&qn->scaling);
220: }
221: }
223: VecCopy(D,Y);
225: if (qn->singlereduction) {
226: VecMDot(Y,l,dX,YtdX);
227: }
228: /* outward recursion starting at iteration k's update and working back */
229: for (i=0; i<l; i++) {
230: k = (it-i-1)%l;
231: if (qn->singlereduction) {
232: /* construct t = dX[k] dot Y as Y_0 dot dX[k] + sum(-alpha[j]dX[k]dF[j]) */
233: t = YtdX[k];
234: for (j=0; j<i; j++) {
235: g = (it-j-1)%l;
236: t += -alpha[g]*H(g, k);
237: }
238: alpha[k] = t/H(k,k);
239: } else {
240: VecDot(dX[k],Y,&t);
241: alpha[k] = t/dXtdF[k];
242: }
243: if (qn->monitor) {
244: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
245: PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d alpha: %14.12e\n", it, k, PetscRealPart(alpha[k]));
246: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
247: }
248: VecAXPY(Y,-alpha[k],dF[k]);
249: }
251: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
252: KSPSolve(snes->ksp,Y,W);
253: KSPGetConvergedReason(snes->ksp,&kspreason);
254: if (kspreason < 0) {
255: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
256: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
257: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
258: return(0);
259: }
260: }
261: KSPGetIterationNumber(snes->ksp,&lits);
262: snes->linear_its += lits;
263: VecCopy(W, Y);
264: } else {
265: VecScale(Y, qn->scaling);
266: }
267: if (qn->singlereduction) {
268: VecMDot(Y,l,dF,YtdX);
269: }
270: /* inward recursion starting at the first update and working forward */
271: for (i = 0; i < l; i++) {
272: k = (it + i - l) % l;
273: if (qn->singlereduction) {
274: t = YtdX[k];
275: for (j = 0; j < i; j++) {
276: g = (it + j - l) % l;
277: t += (alpha[g] - beta[g])*H(k, g);
278: }
279: beta[k] = t / H(k, k);
280: } else {
281: VecDot(dF[k], Y, &t);
282: beta[k] = t / dXtdF[k];
283: }
284: VecAXPY(Y, (alpha[k] - beta[k]), dX[k]);
285: if (qn->monitor) {
286: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
287: PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d alpha - beta: %14.12e\n", it, k, PetscRealPart(alpha[k] - beta[k]));
288: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
289: }
290: }
291: return(0);
292: }
296: static PetscErrorCode SNESSolve_QN(SNES snes)
297: {
298: PetscErrorCode ierr;
299: SNES_QN *qn = (SNES_QN*) snes->data;
300: Vec X,Xold;
301: Vec F,B;
302: Vec Y,FPC,D,Dold;
303: SNESConvergedReason reason;
304: PetscInt i, i_r;
305: PetscReal fnorm,xnorm,ynorm,gnorm;
306: PetscBool lssucceed,powell,periodic;
307: PetscScalar DolddotD,DolddotDold,DdotD,YdotD;
308: MatStructure flg = DIFFERENT_NONZERO_PATTERN;
310: /* basically just a regular newton's method except for the application of the jacobian */
313: F = snes->vec_func; /* residual vector */
314: Y = snes->vec_sol_update; /* search direction generated by J^-1D*/
315: B = snes->vec_rhs;
317: X = snes->vec_sol; /* solution vector */
318: Xold = snes->work[0];
320: /* directions generated by the preconditioned problem with F_pre = F or x - M(x, b) */
321: D = snes->work[1];
322: Dold = snes->work[2];
324: snes->reason = SNES_CONVERGED_ITERATING;
326: PetscObjectAMSTakeAccess((PetscObject)snes);
327: snes->iter = 0;
328: snes->norm = 0.;
329: PetscObjectAMSGrantAccess((PetscObject)snes);
330: if (!snes->vec_func_init_set) {
331: SNESComputeFunction(snes,X,F);
332: if (snes->domainerror) {
333: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
334: return(0);
335: }
336: } else snes->vec_func_init_set = PETSC_FALSE;
338: if (!snes->norm_init_set) {
339: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
340: if (PetscIsInfOrNanReal(fnorm)) {
341: snes->reason = SNES_DIVERGED_FNORM_NAN;
342: return(0);
343: }
344: } else {
345: fnorm = snes->norm_init;
346: snes->norm_init_set = PETSC_FALSE;
347: }
349: PetscObjectAMSTakeAccess((PetscObject)snes);
350: snes->norm = fnorm;
351: PetscObjectAMSGrantAccess((PetscObject)snes);
352: SNESLogConvergenceHistory(snes,fnorm,0);
353: SNESMonitor(snes,0,fnorm);
355: /* set parameter for default relative tolerance convergence test */
356: snes->ttol = fnorm*snes->rtol;
357: /* test convergence */
358: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
359: if (snes->reason) return(0);
361: /* composed solve */
362: if (snes->pc && snes->pcside == PC_RIGHT) {
363: SNESSetInitialFunction(snes->pc, F);
364: SNESSetInitialFunctionNorm(snes->pc, fnorm);
366: PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,B,0);
367: SNESSolve(snes->pc, B, X);
368: PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,B,0);
370: SNESGetConvergedReason(snes->pc,&reason);
371: if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) {
372: snes->reason = SNES_DIVERGED_INNER;
373: return(0);
374: }
375: SNESGetFunction(snes->pc, &FPC, NULL, NULL);
376: VecCopy(FPC, F);
377: SNESGetFunctionNorm(snes->pc, &fnorm);
378: VecCopy(F, Y);
379: } else {
380: VecCopy(F, Y);
381: }
382: VecCopy(Y, D);
384: /* scale the initial update */
385: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
386: SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
387: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);
388: }
390: for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
391: switch (qn->type) {
392: case SNES_QN_BADBROYDEN:
393: SNESQNApply_BadBroyden(snes,i_r,Y,X,Xold,D,Dold);
394: break;
395: case SNES_QN_BROYDEN:
396: SNESQNApply_Broyden(snes,i_r,Y,X,Xold,D,Dold);
397: break;
398: case SNES_QN_LBFGS:
399: SNESQNApply_LBFGS(snes,i_r,Y,X,Xold,D,Dold);
400: break;
401: }
402: /* line search for lambda */
403: ynorm = 1; gnorm = fnorm;
404: VecCopy(D, Dold);
405: VecCopy(X, Xold);
406: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
407: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
408: if (snes->domainerror) {
409: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
410: return(0);
411: }
412: SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);
413: if (!lssucceed) {
414: if (++snes->numFailures >= snes->maxFailures) {
415: snes->reason = SNES_DIVERGED_LINE_SEARCH;
416: break;
417: }
418: }
419: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
420: if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
421: SNESLineSearchGetLambda(snes->linesearch, &qn->scaling);
422: }
424: /* convergence monitoring */
425: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);
427: SNESSetIterationNumber(snes, i+1);
428: SNESSetFunctionNorm(snes, fnorm);
430: SNESLogConvergenceHistory(snes,snes->norm,snes->iter);
431: SNESMonitor(snes,snes->iter,snes->norm);
432: /* set parameter for default relative tolerance convergence test */
433: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
434: if (snes->reason) return(0);
436: if (snes->pc && snes->pcside == PC_RIGHT) {
437: SNESSetInitialFunction(snes->pc, F);
438: SNESSetInitialFunctionNorm(snes->pc, fnorm);
439: PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,B,0);
440: SNESSolve(snes->pc, B, X);
441: PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,B,0);
442: SNESGetConvergedReason(snes->pc,&reason);
443: if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) {
444: snes->reason = SNES_DIVERGED_INNER;
445: return(0);
446: }
447: SNESGetFunction(snes->pc, &FPC, NULL, NULL);
448: VecCopy(FPC, F);
449: SNESGetFunctionNorm(snes->pc, &fnorm);
450: VecCopy(F, D);
451: } else {
452: VecCopy(F, D);
453: }
455: powell = PETSC_FALSE;
456: if (qn->restart_type == SNES_QN_RESTART_POWELL) {
457: /* check restart by Powell's Criterion: |F^T H_0 Fold| > 0.2 * |Fold^T H_0 Fold| */
458: VecDotBegin(Dold, Dold, &DolddotDold);
459: VecDotBegin(Dold, D, &DolddotD);
460: VecDotBegin(D, D, &DdotD);
461: VecDotBegin(Y, D, &YdotD);
462: VecDotEnd(Dold, Dold, &DolddotDold);
463: VecDotEnd(Dold, D, &DolddotD);
464: VecDotEnd(D, D, &DdotD);
465: VecDotEnd(Y, D, &YdotD);
466: if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE;
467: }
468: periodic = PETSC_FALSE;
469: if (qn->restart_type == SNES_QN_RESTART_PERIODIC) {
470: if (i_r>qn->m-1) periodic = PETSC_TRUE;
471: }
472: /* restart if either powell or periodic restart is satisfied. */
473: if (powell || periodic) {
474: if (qn->monitor) {
475: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
476: PetscViewerASCIIPrintf(qn->monitor, "restart! |%14.12e| > %4.2f*|%14.12e| or i_r = %d\n", PetscRealPart(DolddotD), qn->powell_gamma, PetscRealPart(DolddotDold), i_r);
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,&flg);
486: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);
487: }
488: }
489: /* general purpose update */
490: if (snes->ops->update) {
491: (*snes->ops->update)(snes, snes->iter);
492: }
493: }
494: if (i == snes->max_its) {
495: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);
496: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
497: }
498: return(0);
499: }
503: static PetscErrorCode SNESSetUp_QN(SNES snes)
504: {
505: SNES_QN *qn = (SNES_QN*)snes->data;
509: VecDuplicateVecs(snes->vec_sol, qn->m, &qn->U);
510: VecDuplicateVecs(snes->vec_sol, qn->m, &qn->V);
511: PetscMalloc3(qn->m, PetscScalar, &qn->alpha,
512: qn->m, PetscScalar, &qn->beta,
513: qn->m, PetscScalar, &qn->dXtdF);
515: if (qn->singlereduction) {
516: PetscMalloc3(qn->m*qn->m, PetscScalar, &qn->dXdFmat,
517: qn->m, PetscScalar, &qn->dFtdX,
518: qn->m, PetscScalar, &qn->YtdX);
519: }
520: SNESSetWorkVecs(snes,4);
522: /* set up the line search */
523: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
524: SNESSetUpMatrices(snes);
525: }
526: return(0);
527: }
531: static PetscErrorCode SNESReset_QN(SNES snes)
532: {
534: SNES_QN *qn;
537: if (snes->data) {
538: qn = (SNES_QN*)snes->data;
539: if (qn->U) {
540: VecDestroyVecs(qn->m, &qn->U);
541: }
542: if (qn->V) {
543: VecDestroyVecs(qn->m, &qn->V);
544: }
545: if (qn->singlereduction) {
546: PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);
547: }
548: PetscFree3(qn->alpha, qn->beta, qn->dXtdF);
549: }
550: return(0);
551: }
555: static PetscErrorCode SNESDestroy_QN(SNES snes)
556: {
560: SNESReset_QN(snes);
561: PetscFree(snes->data);
562: PetscObjectComposeFunction((PetscObject)snes,"",NULL);
563: return(0);
564: }
568: static PetscErrorCode SNESSetFromOptions_QN(SNES snes)
569: {
571: PetscErrorCode ierr;
572: SNES_QN *qn = (SNES_QN*)snes->data;
573: PetscBool monflg = PETSC_FALSE,flg;
574: SNESLineSearch linesearch;
575: SNESQNRestartType rtype = qn->restart_type;
576: SNESQNScaleType stype = qn->scale_type;
579: PetscOptionsHead("SNES QN options");
580: PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,NULL);
581: PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL);
582: PetscOptionsReal("-snes_qn_powell_downhill","Powell descent tolerance", "SNESQN", qn->powell_downhill, &qn->powell_downhill, NULL);
583: PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", monflg, &monflg, NULL);
584: PetscOptionsBool("-snes_qn_single_reduction", "Aggregate reductions", "SNESQN", qn->singlereduction, &qn->singlereduction, NULL);
585: PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);
586: if (flg) SNESQNSetScaleType(snes,stype);
588: PetscOptionsEnum("-snes_qn_restart_type","Restart type","SNESQNSetRestartType",SNESQNRestartTypes,(PetscEnum)rtype,(PetscEnum*)&rtype,&flg);
589: if (flg) SNESQNSetRestartType(snes,rtype);
591: PetscOptionsEnum("-snes_qn_type","Quasi-Newton update type","",SNESQNTypes,
592: (PetscEnum)qn->type,(PetscEnum*)&qn->type,NULL);
593: PetscOptionsTail();
594: if (!snes->linesearch) {
595: SNESGetLineSearch(snes, &linesearch);
596: SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);
597: }
598: if (monflg) {
599: qn->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
600: }
601: return(0);
602: }
607: /*@
608: SNESQNSetRestartType - Sets the restart type for SNESQN.
610: Logically Collective on SNES
612: Input Parameters:
613: + snes - the iterative context
614: - rtype - restart type
616: Options Database:
617: + -snes_qn_restart_type<powell,periodic,none> - set the restart type
618: - -snes_qn_restart[10] - sets the number of iterations before restart for periodic
620: Level: intermediate
622: SNESQNRestartTypes:
623: + SNES_QN_RESTART_NONE - never restart
624: . SNES_QN_RESTART_POWELL - restart based upon descent criteria
625: - SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations
627: Notes:
628: The default line search used is the L2 line search and it requires two additional function evaluations.
630: .keywords: SNES, SNESQN, restart, type, set SNESLineSearch
631: @*/
632: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
633: {
638: PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));
639: return(0);
640: }
644: /*@
645: SNESQNSetScaleType - Sets the scaling type for the inner inverse jacobian in SNESQN.
647: Logically Collective on SNES
649: Input Parameters:
650: + snes - the iterative context
651: - stype - scale type
653: Options Database:
654: . -snes_qn_scale_type<shanno,none,linesearch,jacobian>
656: Level: intermediate
658: SNESQNSelectTypes:
659: + SNES_QN_SCALE_NONE - don't scale the problem
660: . SNES_QN_SCALE_SHANNO - use shanno scaling
661: . SNES_QN_SCALE_LINESEARCH - scale based upon line search lambda
662: - SNES_QN_SCALE_JACOBIAN - scale by inverting a previously computed Jacobian.
664: .keywords: SNES, SNESQN, scaling, type, set SNESLineSearch
665: @*/
667: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
668: {
673: PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));
674: return(0);
675: }
679: PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
680: {
681: SNES_QN *qn = (SNES_QN*)snes->data;
684: qn->scale_type = stype;
685: return(0);
686: }
690: PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
691: {
692: SNES_QN *qn = (SNES_QN*)snes->data;
695: qn->restart_type = rtype;
696: return(0);
697: }
699: /* -------------------------------------------------------------------------- */
700: /*MC
701: SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.
703: Options Database:
705: + -snes_qn_m - Number of past states saved for the L-Broyden methods.
706: . -snes_qn_powell_angle - Angle condition for restart.
707: . -snes_qn_powell_descent - Descent condition for restart.
708: . -snes_linesearch_type <cp, l2, basic> - Type of line search.
709: - -snes_qn_monitor - Monitors the quasi-newton jacobian.
711: Notes: This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of F(x) = b using
712: previous change in F(x) and x to form the approximate inverse Jacobian using a series of multiplicative rank-one
713: updates.
715: When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of
716: these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this
717: iteration as the current iteration's values when constructing the approximate jacobian. The second, composed,
718: perturbs the problem the jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner.
720: References:
722: L-Broyden Methods: a generalization of the L-BFGS method to the limited memory Broyden family, M. B. Reed,
723: International Journal of Computer Mathematics, vol. 86, 2009.
725: Experiments with Quasi-Newton Methods in Solving Stiff ODE Systems, Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker
726: SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985.
729: Level: beginner
731: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR
733: M*/
736: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
737: {
739: SNES_QN *qn;
742: snes->ops->setup = SNESSetUp_QN;
743: snes->ops->solve = SNESSolve_QN;
744: snes->ops->destroy = SNESDestroy_QN;
745: snes->ops->setfromoptions = SNESSetFromOptions_QN;
746: snes->ops->view = 0;
747: snes->ops->reset = SNESReset_QN;
749: snes->usespc = PETSC_TRUE;
750: snes->usesksp = PETSC_FALSE;
752: if (!snes->tolerancesset) {
753: snes->max_funcs = 30000;
754: snes->max_its = 10000;
755: }
757: PetscNewLog(snes,SNES_QN,&qn);
758: snes->data = (void*) qn;
759: qn->m = 10;
760: qn->scaling = 1.0;
761: qn->U = NULL;
762: qn->V = NULL;
763: qn->dXtdF = NULL;
764: qn->dFtdX = NULL;
765: qn->dXdFmat = NULL;
766: qn->monitor = NULL;
767: qn->singlereduction = PETSC_FALSE;
768: qn->powell_gamma = 0.9999;
769: qn->scale_type = SNES_QN_SCALE_SHANNO;
770: qn->restart_type = SNES_QN_RESTART_POWELL;
771: qn->type = SNES_QN_LBFGS;
773: PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);
774: PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);
775: return(0);
776: }