Actual source code: qn.c
petsc-3.3-p7 2013-05-11
1: #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/
3: #define H(i,j) qn->dXdFmat[i*qn->m + j]
5: const char *SNESQNCompositionTypes[] = {"SEQUENTIAL","COMPOSED","SNESQNCompositionType","SNES_QN_",0};
6: const char *SNESQNScaleTypes[] = {"NONE","SHANNO","LINESEARCH","JACOBIAN","SNESQNScaleType","SNES_QN_SCALING_",0};
7: const char *SNESQNRestartTypes[] = {"NONE","POWELL","PERIODIC","SNESQNRestartType","SNES_QN_RESTART_",0};
9: typedef struct {
10: Vec *dX; /* The change in X */
11: Vec *dF; /* The change in F */
12: PetscInt m; /* The number of kept previous steps */
13: PetscScalar *alpha, *beta;
14: PetscScalar *dXtdF, *dFtdX, *YtdX;
15: PetscBool singlereduction; /* Aggregated reduction implementation */
16: PetscScalar *dXdFmat; /* A matrix of values for dX_i dot dF_j */
17: PetscViewer monitor;
18: PetscReal powell_gamma; /* Powell angle restart condition */
19: PetscReal powell_downhill; /* Powell descent restart condition */
20: PetscReal scaling; /* scaling of H0 */
21: PetscInt restart_periodic; /* the maximum iterations between restart */
23: SNESQNCompositionType composition_type; /* determine if the composition is done sequentially or as a composition */
24: SNESQNScaleType scale_type; /* determine if the composition is done sequentially or as a composition */
25: SNESQNRestartType restart_type; /* determine the frequency and type of restart conditions */
26: } SNES_QN;
30: PetscErrorCode SNESQNApplyJinv_Private(SNES snes, PetscInt it, Vec D, Vec Y) {
34: SNES_QN *qn = (SNES_QN*)snes->data;
36: Vec Yin = snes->work[3];
38: Vec *dX = qn->dX;
39: Vec *dF = qn->dF;
41: PetscScalar *alpha = qn->alpha;
42: PetscScalar *beta = qn->beta;
43: PetscScalar *dXtdF = qn->dXtdF;
44: PetscScalar *YtdX = qn->YtdX;
46: /* ksp thing for jacobian scaling */
47: KSPConvergedReason kspreason;
48: MatStructure flg = DIFFERENT_NONZERO_PATTERN;
50: PetscInt k, i, j, g, lits;
51: PetscInt m = qn->m;
52: PetscScalar t;
53: PetscInt l = m;
55: Mat jac, jac_pre;
59: VecCopy(D, Y);
61: if (it < m) l = it;
63: if (qn->singlereduction) {
64: VecMDot(Y, l, qn->dX, YtdX);
65: }
66: /* outward recursion starting at iteration k's update and working back */
67: for (i = 0; i < l; i++) {
68: k = (it - i - 1) % l;
69: if (qn->singlereduction) {
70: /* construct t = dX[k] dot Y as Y_0 dot dX[k] + sum(-alpha[j]dX[k]dF[j]) */
71: t = YtdX[k];
72: for (j = 0; j < i; j++) {
73: g = (it - j - 1) % l;
74: t += -alpha[g]*H(g, k);
75: }
76: alpha[k] = t / H(k, k);
77: } else {
78: VecDot(dX[k], Y, &t);
79: alpha[k] = t / dXtdF[k];
80: }
81: if (qn->monitor) {
82: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
83: PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d alpha: %14.12e\n", it, k, PetscRealPart(alpha[k]));
84: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
85: }
86: VecAXPY(Y, -alpha[k], dF[k]);
87: }
89: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
90: SNESGetJacobian(snes, &jac, &jac_pre, PETSC_NULL, PETSC_NULL);
91: KSPSetOperators(snes->ksp,jac,jac_pre,flg);
92: SNES_KSPSolve(snes,snes->ksp,Y,Yin);
93: KSPGetConvergedReason(snes->ksp,&kspreason);
94: if (kspreason < 0) {
95: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
96: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
97: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
98: return(0);
99: }
100: }
101: KSPGetIterationNumber(snes->ksp,&lits);
102: snes->linear_its += lits;
103: VecCopy(Yin, Y);
104: } else {
105: VecScale(Y, qn->scaling);
106: }
107: if (qn->singlereduction) {
108: VecMDot(Y, l, qn->dF, YtdX);
109: }
110: /* inward recursion starting at the first update and working forward */
111: for (i = 0; i < l; i++) {
112: k = (it + i - l) % l;
113: if (qn->singlereduction) {
114: t = YtdX[k];
115: for (j = 0; j < i; j++) {
116: g = (it + j - l) % l;
117: t += (alpha[g] - beta[g])*H(k, g);
118: }
119: beta[k] = t / H(k, k);
120: } else {
121: VecDot(dF[k], Y, &t);
122: beta[k] = t / dXtdF[k];
123: }
124: VecAXPY(Y, (alpha[k] - beta[k]), dX[k]);
125: if (qn->monitor) {
126: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
127: PetscViewerASCIIPrintf(qn->monitor, "it: %d k: %d alpha - beta: %14.12e\n", it, k, PetscRealPart(alpha[k] - beta[k]));
128: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
129: }
130: }
131: return(0);
132: }
136: static PetscErrorCode SNESSolve_QN(SNES snes)
137: {
140: SNES_QN *qn = (SNES_QN*) snes->data;
142: Vec X, Xold;
143: Vec F, B;
144: Vec Y, FPC, D, Dold;
145: SNESConvergedReason reason;
146: PetscInt i, i_r, k, l, j;
148: PetscReal fnorm, xnorm, ynorm, gnorm;
149: PetscInt m = qn->m;
150: PetscBool lssucceed,powell,periodic;
152: Vec *dX = qn->dX;
153: Vec *dF = qn->dF;
154: PetscScalar *dXtdF = qn->dXtdF;
155: PetscScalar *dFtdX = qn->dFtdX;
156: PetscScalar DolddotD, DolddotDold, DdotD, YdotD, a;
158: MatStructure flg = DIFFERENT_NONZERO_PATTERN;
160: /* basically just a regular newton's method except for the application of the jacobian */
163: X = snes->vec_sol; /* solution vector */
164: F = snes->vec_func; /* residual vector */
165: Y = snes->vec_sol_update; /* search direction generated by J^-1D*/
166: B = snes->vec_rhs;
167: Xold = snes->work[0];
169: /* directions generated by the preconditioned problem with F_pre = F or x - M(x, b) */
170: D = snes->work[1];
171: Dold = snes->work[2];
173: snes->reason = SNES_CONVERGED_ITERATING;
175: PetscObjectTakeAccess(snes);
176: snes->iter = 0;
177: snes->norm = 0.;
178: PetscObjectGrantAccess(snes);
179: if (!snes->vec_func_init_set){
180: SNESComputeFunction(snes,X,F);
181: if (snes->domainerror) {
182: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
183: return(0);
184: }
185: } else {
186: snes->vec_func_init_set = PETSC_FALSE;
187: }
189: if (!snes->norm_init_set) {
190: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
191: if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
192: } else {
193: fnorm = snes->norm_init;
194: snes->norm_init_set = PETSC_FALSE;
195: }
197: PetscObjectTakeAccess(snes);
198: snes->norm = fnorm;
199: PetscObjectGrantAccess(snes);
200: SNESLogConvHistory(snes,fnorm,0);
201: SNESMonitor(snes,0,fnorm);
203: /* set parameter for default relative tolerance convergence test */
204: snes->ttol = fnorm*snes->rtol;
205: /* test convergence */
206: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
207: if (snes->reason) return(0);
209: /* composed solve -- either sequential or composed */
210: if (snes->pc) {
211: if (qn->composition_type == SNES_QN_SEQUENTIAL) {
212: SNESSetInitialFunction(snes->pc, F);
213: SNESSetInitialFunctionNorm(snes->pc, fnorm);
214: SNESSolve(snes->pc, B, X);
215: SNESGetConvergedReason(snes->pc,&reason);
216: if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) {
217: snes->reason = SNES_DIVERGED_INNER;
218: return(0);
219: }
220: SNESGetFunction(snes->pc, &FPC, PETSC_NULL, PETSC_NULL);
221: VecCopy(FPC, F);
222: SNESGetFunctionNorm(snes->pc, &fnorm);
223: VecCopy(F, Y);
224: } else {
225: VecCopy(X, Y);
226: SNESSetInitialFunction(snes->pc, F);
227: SNESSetInitialFunctionNorm(snes->pc, fnorm);
228: SNESSolve(snes->pc, B, Y);
229: SNESGetConvergedReason(snes->pc,&reason);
230: if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) {
231: snes->reason = SNES_DIVERGED_INNER;
232: return(0);
233: }
234: VecAYPX(Y,-1.0,X);
235: }
236: } else {
237: VecCopy(F, Y);
238: }
239: VecCopy(Y, D);
241: /* scale the initial update */
242: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
243: SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
244: }
246: for(i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
247: SNESQNApplyJinv_Private(snes, i_r, D, Y);
248: /* line search for lambda */
249: ynorm = 1; gnorm = fnorm;
250: VecCopy(D, Dold);
251: VecCopy(X, Xold);
252: SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
253: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
254: if (snes->domainerror) {
255: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
256: return(0);
257: }
258: SNESLineSearchGetSuccess(snes->linesearch, &lssucceed);
259: if (!lssucceed) {
260: if (++snes->numFailures >= snes->maxFailures) {
261: snes->reason = SNES_DIVERGED_LINE_SEARCH;
262: break;
263: }
264: }
265: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
266: if (qn->scale_type == SNES_QN_SCALE_LINESEARCH) {
267: SNESLineSearchGetLambda(snes->linesearch, &qn->scaling);
268: }
270: /* convergence monitoring */
271: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);
273: SNESSetIterationNumber(snes, i+1);
274: SNESSetFunctionNorm(snes, fnorm);
276: SNESLogConvHistory(snes,snes->norm,snes->iter);
277: SNESMonitor(snes,snes->iter,snes->norm);
278: /* set parameter for default relative tolerance convergence test */
279: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
280: if (snes->reason) return(0);
283: if (snes->pc) {
284: if (qn->composition_type == SNES_QN_SEQUENTIAL) {
285: SNESSetInitialFunction(snes->pc, F);
286: SNESSetInitialFunctionNorm(snes->pc, fnorm);
287: SNESSolve(snes->pc, B, X);
288: SNESGetConvergedReason(snes->pc,&reason);
289: if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) {
290: snes->reason = SNES_DIVERGED_INNER;
291: return(0);
292: }
293: SNESGetFunction(snes->pc, &FPC, PETSC_NULL, PETSC_NULL);
294: VecCopy(FPC, F);
295: SNESGetFunctionNorm(snes->pc, &fnorm);
296: VecCopy(F, D);
297: } else {
298: VecCopy(X, D);
299: SNESSetInitialFunction(snes->pc, F);
300: SNESSetInitialFunctionNorm(snes->pc, fnorm);
301: SNESSolve(snes->pc, B, D);
302: SNESGetConvergedReason(snes->pc,&reason);
303: if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) {
304: snes->reason = SNES_DIVERGED_INNER;
305: return(0);
306: }
307: VecAYPX(D,-1.0,X);
308: }
309: } else {
310: VecCopy(F, D);
311: }
313: powell = PETSC_FALSE;
314: if (qn->restart_type == SNES_QN_RESTART_POWELL) {
315: /* check restart by Powell's Criterion: |F^T H_0 Fold| > 0.2 * |Fold^T H_0 Fold| */
316: VecDotBegin(Dold, Dold, &DolddotDold);
317: VecDotBegin(Dold, D, &DolddotD);
318: VecDotBegin(D, D, &DdotD);
319: VecDotBegin(Y, D, &YdotD);
320: VecDotEnd(Dold, Dold, &DolddotDold);
321: VecDotEnd(Dold, D, &DolddotD);
322: VecDotEnd(D, D, &DdotD);
323: VecDotEnd(Y, D, &YdotD);
324: if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma*PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE;
325: }
326: periodic = PETSC_FALSE;
327: if (qn->restart_type != SNES_QN_RESTART_NONE) {
328: if ((i_r > qn->restart_periodic - 1 && qn->restart_periodic > 0)) periodic = PETSC_TRUE;
329: }
331: /* restart if either powell or periodic restart is satisfied. */
332: if (powell || periodic) {
333: if (qn->monitor) {
334: PetscViewerASCIIAddTab(qn->monitor,((PetscObject)snes)->tablevel+2);
335: PetscViewerASCIIPrintf(qn->monitor, "restart! |%14.12e| > %4.2f*|%14.12e| or i_r = %d\n", PetscRealPart(DolddotD), qn->powell_gamma, PetscRealPart(DolddotDold), i_r);
336: PetscViewerASCIISubtractTab(qn->monitor,((PetscObject)snes)->tablevel+2);
337: }
338: i_r = -1;
339: /* general purpose update */
340: if (snes->ops->update) {
341: (*snes->ops->update)(snes, snes->iter);
342: }
343: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
344: SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
345: }
346: } else {
347: /* set the differences */
348: k = i_r % m;
349: l = m;
350: if (i_r + 1 < m) l = i_r + 1;
351: VecCopy(D, dF[k]);
352: VecAXPY(dF[k], -1.0, Dold);
353: VecCopy(X, dX[k]);
354: VecAXPY(dX[k], -1.0, Xold);
355: if (qn->singlereduction) {
356: VecMDot(dF[k], l, dX, dXtdF);
357: VecMDot(dX[k], l, dF, dFtdX);
358: for (j = 0; j < l; j++) {
359: H(k, j) = dFtdX[j];
360: H(j, k) = dXtdF[j];
361: }
362: /* copy back over to make the computation of alpha and beta easier */
363: for (j = 0; j < l; j++) {
364: dXtdF[j] = H(j, j);
365: }
366: } else {
367: VecDot(dX[k], dF[k], &dXtdF[k]);
368: }
369: /* set scaling to be shanno scaling */
370: if (qn->scale_type == SNES_QN_SCALE_SHANNO) {
371: VecDot(dF[k], dF[k], &a);
372: qn->scaling = PetscRealPart(dXtdF[k]) / PetscRealPart(a);
373: }
374: /* general purpose update */
375: if (snes->ops->update) {
376: (*snes->ops->update)(snes, snes->iter);
377: }
378: }
379: }
380: if (i == snes->max_its) {
381: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", snes->max_its);
382: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
383: }
384: return(0);
385: }
390: static PetscErrorCode SNESSetUp_QN(SNES snes)
391: {
392: SNES_QN *qn = (SNES_QN*)snes->data;
396: VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dX);
397: VecDuplicateVecs(snes->vec_sol, qn->m, &qn->dF);
398: PetscMalloc3(qn->m, PetscScalar, &qn->alpha,
399: qn->m, PetscScalar, &qn->beta,
400: qn->m, PetscScalar, &qn->dXtdF);
402: if (qn->singlereduction) {
403: PetscMalloc3(qn->m*qn->m, PetscScalar, &qn->dXdFmat,
404: qn->m, PetscScalar, &qn->dFtdX,
405: qn->m, PetscScalar, &qn->YtdX);
406: }
407: SNESDefaultGetWork(snes,4);
409: /* set up the line search */
410: if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
411: SNESSetUpMatrices(snes);
412: }
413: return(0);
414: }
418: static PetscErrorCode SNESReset_QN(SNES snes)
419: {
421: SNES_QN *qn;
423: if (snes->data) {
424: qn = (SNES_QN*)snes->data;
425: if (qn->dX) {
426: VecDestroyVecs(qn->m, &qn->dX);
427: }
428: if (qn->dF) {
429: VecDestroyVecs(qn->m, &qn->dF);
430: }
431: if (qn->singlereduction) {
432: PetscFree3(qn->dXdFmat, qn->dFtdX, qn->YtdX);
433: }
434: PetscFree3(qn->alpha, qn->beta, qn->dXtdF);
435: }
436: return(0);
437: }
441: static PetscErrorCode SNESDestroy_QN(SNES snes)
442: {
445: SNESReset_QN(snes);
446: PetscFree(snes->data);
447: PetscObjectComposeFunctionDynamic((PetscObject)snes,"","",PETSC_NULL);
448: return(0);
449: }
453: static PetscErrorCode SNESSetFromOptions_QN(SNES snes)
454: {
457: SNES_QN *qn;
458: PetscBool monflg = PETSC_FALSE;
459: SNESLineSearch linesearch;
462: qn = (SNES_QN*)snes->data;
464: PetscOptionsHead("SNES QN options");
465: PetscOptionsInt("-snes_qn_m","Number of past states saved for L-BFGS methods","SNESQN",qn->m,&qn->m,PETSC_NULL);
466: PetscOptionsInt("-snes_qn_restart","Maximum number of iterations between restarts","SNESQN",qn->restart_periodic,&qn->restart_periodic, PETSC_NULL);
467: PetscOptionsReal("-snes_qn_powell_gamma","Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, PETSC_NULL);
468: PetscOptionsReal("-snes_qn_powell_downhill","Powell descent tolerance", "SNESQN", qn->powell_downhill, &qn->powell_downhill, PETSC_NULL);
469: PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", monflg, &monflg, PETSC_NULL);
470: PetscOptionsBool("-snes_qn_single_reduction", "Aggregate reductions", "SNESQN", qn->singlereduction, &qn->singlereduction, PETSC_NULL);
471: PetscOptionsEnum("-snes_qn_scale_type","Scaling type","SNESQNSetScaleType",SNESQNScaleTypes,(PetscEnum)qn->scale_type,(PetscEnum*)&qn->scale_type,PETSC_NULL);
472: PetscOptionsEnum("-snes_qn_composition_type","Composition type","SNESQNSetCompositionType",SNESQNCompositionTypes,
473: (PetscEnum)qn->composition_type,(PetscEnum*)&qn->composition_type,PETSC_NULL);
474: PetscOptionsEnum("-snes_qn_restart_type","Restart type","SNESQNSetRestartType",SNESQNRestartTypes,
475: (PetscEnum)qn->restart_type,(PetscEnum*)&qn->restart_type,PETSC_NULL);
476: PetscOptionsTail();
477: if (!snes->linesearch) {
478: SNESGetSNESLineSearch(snes, &linesearch);
479: if (!snes->pc || qn->composition_type == SNES_QN_SEQUENTIAL) {
480: SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);
481: } else {
482: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
483: }
484: }
485: if (monflg) {
486: qn->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);
487: }
488: return(0);
489: }
494: /*@
495: SNESQNSetRestartType - Sets the restart type for SNESQN.
497: Logically Collective on SNES
499: Input Parameters:
500: + snes - the iterative context
501: - rtype - restart type
503: Options Database:
504: + -snes_qn_restart_type<powell,periodic,none> - set the restart type
505: - -snes_qn_restart[30] - sets the number of iterations before restart for periodic
507: Level: intermediate
509: SNESQNRestartTypes:
510: + SNES_QN_RESTART_NONE - never restart
511: . SNES_QN_RESTART_POWELL - restart based upon descent criteria
512: - SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations
514: Notes:
515: The default line search used is the L2 line search and it requires two additional function evaluations.
517: .keywords: SNES, SNESQN, restart, type, set SNESLineSearch
518: @*/
519: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype) {
523: PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));
524: return(0);
525: }
529: /*@
530: SNESQNSetScaleType - Sets the scaling type for the inner inverse jacobian in SNESQN.
532: Logically Collective on SNES
534: Input Parameters:
535: + snes - the iterative context
536: - stype - scale type
538: Options Database:
539: . -snes_qn_scale_type<shanno,none,linesearch,jacobian>
541: Level: intermediate
543: SNESQNSelectTypes:
544: + SNES_QN_SCALE_NONE - don't scale the problem
545: . SNES_QN_SCALE_SHANNO - use shanno scaling
546: . SNES_QN_SCALE_LINESEARCH - scale based upon line search lambda
547: - SNES_QN_SCALE_JACOBIAN - scale by inverting a previously computed Jacobian.
549: .keywords: SNES, SNESQN, scaling, type, set SNESLineSearch
550: @*/
552: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype) {
556: PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));
557: return(0);
558: }
562: /*@
563: SNESQNSetCompositionType - Sets the composition type
565: Logically Collective on SNES
567: Input Parameters:
568: + snes - the iterative context
569: - stype - composition type
571: Options Database:
572: . -snes_qn_composition_type<sequential, composed>
574: Level: intermediate
576: SNESQNSelectTypes:
577: + SNES_QN_COMPOSITION_SEQUENTIAL - Solve the system with X = PC(X) and D = F(PC(X))
578: - SNES_QN_COMPOSITION_COMPOSED - solve the system with X = X and D = PC(X) - X
580: .keywords: SNES, SNESQN, scaling, type, set SNESLineSearch
581: @*/
583: PetscErrorCode SNESQNSetCompositionType(SNES snes, SNESQNCompositionType ctype) {
587: PetscTryMethod(snes,"SNESQNSetCompositionType_C",(SNES,SNESQNCompositionType),(snes,ctype));
588: return(0);
589: }
591: EXTERN_C_BEGIN
594: PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype) {
595: SNES_QN *qn = (SNES_QN *)snes->data;
597: qn->scale_type = stype;
598: return(0);
599: }
603: PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype) {
604: SNES_QN *qn = (SNES_QN *)snes->data;
606: qn->restart_type = rtype;
607: return(0);
608: }
613: PetscErrorCode SNESQNSetCompositionType_QN(SNES snes, SNESQNCompositionType ctype) {
614: SNES_QN *qn = (SNES_QN *)snes->data;
616: qn->composition_type = ctype;
617: return(0);
618: }
619: EXTERN_C_END
621: /* -------------------------------------------------------------------------- */
622: /*MC
623: SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.
625: Options Database:
627: + -snes_qn_m - Number of past states saved for the L-Broyden methods.
628: . -snes_qn_powell_angle - Angle condition for restart.
629: . -snes_qn_powell_descent - Descent condition for restart.
630: . -snes_qn_composition <sequential, composed>- Type of composition.
631: . -snes_linesearch_type <cp, l2, basic> - Type of line search.
632: - -snes_qn_monitor - Monitors the quasi-newton jacobian.
634: Notes: This implements the L-BFGS algorithm for the solution of F(x) = b using previous change in F(x) and x to
635: form the approximate inverse Jacobian using a series of multiplicative rank-one updates. This will eventually be
636: generalized to implement several limited-memory Broyden methods.
638: When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of
639: these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this
640: iteration as the current iteration's values when constructing the approximate jacobian. The second, composed,
641: perturbs the problem the jacobian represents to be P(x, b) - x = 0, where P(x, b) is the preconditioner.
643: References:
645: L-Broyden Methods: a generalization of the L-BFGS method to the limited memory Broyden family, M. B. Reed,
646: International Journal of Computer Mathematics, vol. 86, 2009.
649: Level: beginner
651: .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR
653: M*/
654: EXTERN_C_BEGIN
657: PetscErrorCode SNESCreate_QN(SNES snes)
658: {
661: SNES_QN *qn;
664: snes->ops->setup = SNESSetUp_QN;
665: snes->ops->solve = SNESSolve_QN;
666: snes->ops->destroy = SNESDestroy_QN;
667: snes->ops->setfromoptions = SNESSetFromOptions_QN;
668: snes->ops->view = 0;
669: snes->ops->reset = SNESReset_QN;
671: snes->usespc = PETSC_TRUE;
672: snes->usesksp = PETSC_FALSE;
674: if (!snes->tolerancesset) {
675: snes->max_funcs = 30000;
676: snes->max_its = 10000;
677: }
679: PetscNewLog(snes,SNES_QN,&qn);
680: snes->data = (void *) qn;
681: qn->m = 10;
682: qn->scaling = 1.0;
683: qn->dX = PETSC_NULL;
684: qn->dF = PETSC_NULL;
685: qn->dXtdF = PETSC_NULL;
686: qn->dFtdX = PETSC_NULL;
687: qn->dXdFmat = PETSC_NULL;
688: qn->monitor = PETSC_NULL;
689: qn->singlereduction = PETSC_FALSE;
690: qn->powell_gamma = 0.9;
691: qn->powell_downhill = 0.2;
692: qn->composition_type= SNES_QN_SEQUENTIAL;
693: qn->scale_type = SNES_QN_SCALE_SHANNO;
694: qn->restart_type = SNES_QN_RESTART_POWELL;
695: qn->restart_periodic= -1;
697: return(0);
698: }
700: EXTERN_C_END