Actual source code: ngmresfunc.c
1: #include <../src/snes/impls/ngmres/snesngmres.h>
2: #include <petscblaslapack.h>
4: PetscErrorCode SNESNGMRESGetAdditiveLineSearch_Private(SNES snes, SNESLineSearch *linesearch)
5: {
6: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
8: PetscFunctionBegin;
9: if (!ngmres->additive_linesearch) {
10: const char *optionsprefix;
11: PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
12: PetscCall(SNESLineSearchCreate(PetscObjectComm((PetscObject)snes), &ngmres->additive_linesearch));
13: PetscCall(SNESLineSearchSetSNES(ngmres->additive_linesearch, snes));
14: PetscCall(SNESLineSearchSetType(ngmres->additive_linesearch, SNESLINESEARCHSECANT));
15: PetscCall(SNESLineSearchAppendOptionsPrefix(ngmres->additive_linesearch, "snes_ngmres_additive_"));
16: PetscCall(SNESLineSearchAppendOptionsPrefix(ngmres->additive_linesearch, optionsprefix));
17: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ngmres->additive_linesearch, (PetscObject)snes, 1));
18: }
19: *linesearch = ngmres->additive_linesearch;
20: PetscFunctionReturn(PETSC_SUCCESS);
21: }
23: PetscErrorCode SNESNGMRESUpdateSubspace_Private(SNES snes, PetscInt ivec, PetscInt l, Vec F, PetscReal fnorm, Vec X)
24: {
25: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
26: Vec *Fdot = ngmres->Fdot;
27: Vec *Xdot = ngmres->Xdot;
29: PetscFunctionBegin;
30: PetscCheck(ivec < ngmres->msize, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "Cannot update vector %" PetscInt_FMT " with history size %" PetscInt_FMT "!", ivec, ngmres->msize);
31: PetscCall(VecCopy(F, Fdot[ivec]));
32: PetscCall(VecCopy(X, Xdot[ivec]));
34: ngmres->fnorms[ivec] = fnorm;
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: /*
39: Set the window to the l most recently stored iterates of the circular history, whose newest entry is ivec.
40: The window lists the history indices in increasing order, so it is 0, ..., l - 1 until the history wraps around.
41: */
42: static PetscErrorCode SNESNGMRESSetWindow_Private(SNES snes, PetscInt ivec, PetscInt l)
43: {
44: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
45: PetscInt msize = ngmres->msize, n = 0;
47: PetscFunctionBegin;
48: for (PetscInt s = 0; s < msize; s++) {
49: if ((ivec - s + msize) % msize < l) {
50: ngmres->window[n] = s;
51: ngmres->Fwindow[n] = ngmres->Fdot[s];
52: ngmres->Xwindow[n] = ngmres->Xdot[s];
53: n++;
54: }
55: }
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: PetscErrorCode SNESNGMRESFormCombinedSolution_Private(SNES snes, PetscInt ivec, PetscInt l, Vec XM, Vec FM, PetscReal fMnorm, Vec X, Vec XA, Vec FA)
60: {
61: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
62: Vec *Fwindow = ngmres->Fwindow;
63: Vec *Xwindow = ngmres->Xwindow;
64: PetscInt *window = ngmres->window;
65: PetscScalar *beta = ngmres->beta;
66: PetscScalar *xi = ngmres->xi;
67: PetscScalar alph_total = 0.;
68: PetscReal nu, obj;
69: Vec Y = snes->vec_sol_update;
70: PetscBool changed_y, changed_w;
72: PetscFunctionBegin;
73: nu = fMnorm * fMnorm;
74: PetscCall(SNESNGMRESSetWindow_Private(snes, ivec, l));
76: /* construct the right-hand side and xi factors */
77: if (l > 0) {
78: PetscCall(VecMDotBegin(FM, l, Fwindow, xi));
79: PetscCall(VecMDotBegin(ngmres->Fdot[ivec], l, Fwindow, beta));
80: PetscCall(VecMDotEnd(FM, l, Fwindow, xi));
81: PetscCall(VecMDotEnd(ngmres->Fdot[ivec], l, Fwindow, beta));
82: for (PetscInt i = 0; i < l; i++) {
83: Q(window[i], ivec) = beta[i];
84: Q(ivec, window[i]) = beta[i];
85: }
86: } else {
87: Q(0, 0) = ngmres->fnorms[ivec] * ngmres->fnorms[ivec];
88: }
90: for (PetscInt i = 0; i < l; i++) beta[i] = nu - xi[i];
92: /* construct h */
93: for (PetscInt j = 0; j < l; j++) {
94: for (PetscInt i = 0; i < l; i++) H(i, j) = Q(window[i], window[j]) - xi[i] - xi[j] + nu;
95: }
96: if (l == 1) {
97: /* simply set alpha[0] = beta[0] / H[0, 0] */
98: if (H(0, 0) != 0.) beta[0] = beta[0] / H(0, 0);
99: else beta[0] = 0.;
100: } else {
101: PetscCall(PetscBLASIntCast(l, &ngmres->m));
102: PetscCall(PetscBLASIntCast(l, &ngmres->n));
103: ngmres->rcond = -1.;
104: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
105: #if PetscDefined(USE_COMPLEX)
106: PetscCallLAPACKInfo("LAPACKgelss", LAPACKgelss_(&ngmres->m, &ngmres->n, &ngmres->nrhs, ngmres->h, &ngmres->lda, ngmres->beta, &ngmres->ldb, ngmres->s, &ngmres->rcond, &ngmres->rank, ngmres->work, &ngmres->lwork, ngmres->rwork, &info));
107: #else
108: PetscCallLAPACKInfo("LAPACKgelss", LAPACKgelss_(&ngmres->m, &ngmres->n, &ngmres->nrhs, ngmres->h, &ngmres->lda, ngmres->beta, &ngmres->ldb, ngmres->s, &ngmres->rcond, &ngmres->rank, ngmres->work, &ngmres->lwork, &info));
109: #endif
110: PetscCall(PetscFPTrapPop());
111: }
112: for (PetscInt i = 0; i < l; i++) PetscCheck(!PetscIsInfOrNanScalar(beta[i]), PetscObjectComm((PetscObject)snes), PETSC_ERR_LIB, "SVD generated inconsistent output");
113: alph_total = 0.;
114: for (PetscInt i = 0; i < l; i++) alph_total += beta[i];
116: /*
117: The right-hand side nu - xi lies in the range of h, so the least-squares solution beta satisfies the normal equations,
118: and the squared objective at the minimizer is nu - sum_i beta_i (nu - xi_i).
119: */
120: obj = nu;
121: for (PetscInt i = 0; i < l; i++) obj -= PetscRealPart(beta[i] * (nu - xi[i]));
122: // Floating-point roundoff can make the computed squared norm slightly negative, so clamp it before taking the square root.
123: obj = PetscSqrtReal(PetscMax(obj, 0.0));
124: ngmres->lsnorm = obj;
125: ngmres->gamma = fMnorm > 0.0 ? obj / fMnorm : 0.0;
127: PetscCall(VecAXPBY(XA, 1.0 - alph_total, 0.0, XM));
128: PetscCall(VecMAXPY(XA, l, beta, Xwindow));
129: /* check the validity of the step */
130: PetscCall(VecWAXPY(Y, -1.0, X, XA));
131: PetscCall(SNESLineSearchPostCheck(snes->linesearch, X, Y, XA, &changed_y, &changed_w));
132: if (!ngmres->approxfunc) {
133: if (snes->npc && snes->npcside == PC_LEFT) PetscCall(SNESApplyNPC(snes, XA, NULL, FA));
134: else PetscCall(SNESComputeFunction(snes, XA, FA));
135: } else {
136: PetscCall(VecAXPBY(FA, 1.0 - alph_total, 0.0, FM));
137: PetscCall(VecMAXPY(FA, l, beta, Fwindow));
138: }
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: PetscErrorCode SNESNGMRESNorms_Private(SNES snes, PetscInt l, Vec X, Vec F, Vec XM, Vec FM, Vec XA, Vec FA, Vec D, PetscReal *dnorm, PetscReal *dminnorm, PetscReal *xMnorm, PetscReal *fMnorm, PetscReal *yMnorm, PetscReal *xAnorm, PetscReal *fAnorm, PetscReal *yAnorm)
143: {
144: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
145: PetscReal dcurnorm, dmin = -1.0;
146: Vec *Xwindow = ngmres->Xwindow;
148: PetscFunctionBegin;
149: if (xMnorm) PetscCall(VecNormBegin(XM, NORM_2, xMnorm));
150: if (fMnorm) PetscCall(VecNormBegin(FM, NORM_2, fMnorm));
151: if (yMnorm) {
152: PetscCall(VecWAXPY(D, -1.0, XM, X));
153: PetscCall(VecNormBegin(D, NORM_2, yMnorm));
154: }
155: if (xAnorm) PetscCall(VecNormBegin(XA, NORM_2, xAnorm));
156: if (fAnorm) PetscCall(VecNormBegin(FA, NORM_2, fAnorm));
157: if (yAnorm) {
158: PetscCall(VecWAXPY(D, -1.0, XA, X));
159: PetscCall(VecNormBegin(D, NORM_2, yAnorm));
160: }
161: if (dnorm) {
162: PetscCall(VecWAXPY(D, -1.0, XM, XA));
163: PetscCall(VecNormBegin(D, NORM_2, dnorm));
164: }
165: if (dminnorm) {
166: for (PetscInt i = 0; i < l; i++) {
167: PetscCall(VecWAXPY(D, -1.0, XA, Xwindow[i]));
168: PetscCall(VecNormBegin(D, NORM_2, &ngmres->xnorms[i]));
169: }
170: }
171: if (xMnorm) PetscCall(VecNormEnd(XM, NORM_2, xMnorm));
172: if (fMnorm) PetscCall(VecNormEnd(FM, NORM_2, fMnorm));
173: if (yMnorm) PetscCall(VecNormEnd(D, NORM_2, yMnorm));
174: if (xAnorm) PetscCall(VecNormEnd(XA, NORM_2, xAnorm));
175: if (fAnorm) PetscCall(VecNormEnd(FA, NORM_2, fAnorm));
176: if (yAnorm) PetscCall(VecNormEnd(D, NORM_2, yAnorm));
177: if (dnorm) PetscCall(VecNormEnd(D, NORM_2, dnorm));
178: if (dminnorm) {
179: for (PetscInt i = 0; i < l; i++) {
180: PetscCall(VecNormEnd(D, NORM_2, &ngmres->xnorms[i]));
181: dcurnorm = ngmres->xnorms[i];
182: if (dcurnorm < dmin || dmin < 0.0) dmin = dcurnorm;
183: }
184: *dminnorm = dmin;
185: }
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: PetscErrorCode SNESNGMRESSelect_Private(SNES snes, PetscInt k_restart, Vec XM, Vec FM, PetscReal xMnorm, PetscReal fMnorm, PetscReal yMnorm, PetscReal objM, Vec XA, Vec FA, PetscReal xAnorm, PetscReal fAnorm, PetscReal yAnorm, PetscReal objA, PetscReal dnorm, PetscReal objmin, PetscReal dminnorm, Vec X, Vec F, Vec Y, PetscReal *xnorm, PetscReal *fnorm, PetscReal *ynorm)
190: {
191: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
192: SNESLineSearchReason lsreason;
193: PetscBool selectA;
195: PetscFunctionBegin;
196: if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) {
197: /* X = X + \lambda(XA - X) */
198: if (ngmres->monitor) PetscCall(PetscViewerASCIIPrintf(ngmres->monitor, "obj(X_A) = %e, ||F_A||_2 = %e, obj(X_M) = %e, ||F_M||_2 = %e\n", (double)objA, (double)fAnorm, (double)objM, (double)fMnorm));
199: /* Test if is XA - XM is a descent direction: we want < F(XM), XA - XM > not positive
200: If positive, GMRES will be restarted see https://epubs.siam.org/doi/pdf/10.1137/110835530 */
201: PetscCall(VecCopy(FM, F));
202: PetscCall(VecCopy(XM, X));
203: PetscCall(VecWAXPY(Y, -1.0, XA, X)); /* minus sign since linesearch expects to find Xnew = X - lambda * Y */
204: PetscCall(VecDotRealPart(FM, Y, &ngmres->descent_ls_test)); /* this is actually < F(XM), XM - XA > */
205: *fnorm = fMnorm;
206: if (ngmres->descent_ls_test < 0) { /* XA - XM is not a descent direction, select XM */
207: *xnorm = xMnorm;
208: *fnorm = fMnorm;
209: *ynorm = yMnorm;
210: PetscCall(VecWAXPY(Y, -1.0, X, XM));
211: PetscCall(VecCopy(FM, F));
212: PetscCall(VecCopy(XM, X));
213: } else {
214: PetscCall(SNESNGMRESGetAdditiveLineSearch_Private(snes, &ngmres->additive_linesearch));
215: PetscCall(SNESLineSearchApply(ngmres->additive_linesearch, X, F, fnorm, Y));
216: PetscCall(SNESLineSearchGetReason(ngmres->additive_linesearch, &lsreason));
217: if (lsreason == SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN) {
218: PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNES solver has not converged");
219: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
222: if (lsreason == SNES_LINESEARCH_FAILED_NANORINF) {
223: PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNES solver has not converged");
224: snes->reason = SNES_DIVERGED_FUNCTION_NANORINF;
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
227: if (lsreason && ++snes->numFailures >= snes->maxFailures) {
228: PetscCheck(!snes->errorifnotconverged, PetscObjectComm((PetscObject)snes), PETSC_ERR_NOT_CONVERGED, "SNES solver has not converged");
229: snes->reason = SNES_DIVERGED_LINE_SEARCH;
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
232: PetscCall(SNESLineSearchGetNorms(ngmres->additive_linesearch, xnorm, fnorm, ynorm));
233: }
234: if (ngmres->monitor) {
235: PetscReal objT = *fnorm;
236: SNESObjectiveFn *objective;
238: PetscCall(SNESGetObjective(snes, &objective, NULL));
239: if (objective) PetscCall(SNESComputeObjective(snes, X, &objT));
240: PetscCall(PetscViewerASCIIPrintf(ngmres->monitor, "Additive solution: objective = %e\n", (double)objT));
241: }
242: } else if (ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) {
243: /* Conditions for choosing the accelerated answer:
244: Criterion A -- the objective function isn't increased above the minimum by too much
245: Criterion B -- the choice of x^A isn't too close to some other choice
246: */
247: selectA = (PetscBool)(/* A */ (objA < ngmres->gammaA * objmin) && /* B */ (ngmres->epsilonB * dnorm < dminnorm || objA < ngmres->deltaB * objmin));
249: if (selectA) {
250: if (ngmres->monitor) PetscCall(PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, obj(X_A) = %e, ||F_A||_2 = %e, obj(X_M) = %e, ||F_M||_2 = %e\n", (double)objA, (double)fAnorm, (double)objM, (double)fMnorm));
251: /* copy it over */
252: *xnorm = xAnorm;
253: *fnorm = fAnorm;
254: *ynorm = yAnorm;
255: PetscCall(VecCopy(FA, F));
256: PetscCall(VecCopy(XA, X));
257: } else {
258: if (ngmres->monitor) PetscCall(PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, obj(X_A) = %e, ||F_A||_2 = %e, obj(X_M) = %e, ||F_M||_2 = %e\n", (double)objA, (double)fAnorm, (double)objM, (double)fMnorm));
259: *xnorm = xMnorm;
260: *fnorm = fMnorm;
261: *ynorm = yMnorm;
262: PetscCall(VecWAXPY(Y, -1.0, X, XM));
263: PetscCall(VecCopy(FM, F));
264: PetscCall(VecCopy(XM, X));
265: }
266: } else { /* none */
267: *xnorm = xAnorm;
268: *fnorm = fAnorm;
269: *ynorm = yAnorm;
270: PetscCall(VecCopy(FA, F));
271: PetscCall(VecCopy(XA, X));
272: }
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: PetscErrorCode SNESNGMRESSelectRestart_Private(SNES snes, PetscInt l, PetscReal obj, PetscReal objM, PetscReal objA, PetscReal dnorm, PetscReal objmin, PetscReal dminnorm, PetscBool *selectRestart)
277: {
278: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
280: PetscFunctionBegin;
281: *selectRestart = PETSC_FALSE;
282: if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) {
283: if (ngmres->descent_ls_test < 0) { /* XA - XM is not a descent direction */
284: if (ngmres->monitor) PetscCall(PetscViewerASCIIPrintf(ngmres->monitor, "ascent restart: %e > 0\n", (double)-ngmres->descent_ls_test));
285: *selectRestart = PETSC_TRUE;
286: }
287: } else if (ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) {
288: /* difference stagnation restart */
289: if (ngmres->epsilonB * dnorm > dminnorm && objA > ngmres->deltaB * objmin && l > 0) {
290: if (ngmres->monitor) PetscCall(PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", (double)(ngmres->epsilonB * dnorm), (double)dminnorm));
291: *selectRestart = PETSC_TRUE;
292: }
293: /* residual stagnation restart */
294: if (objA > ngmres->gammaC * objmin) {
295: if (ngmres->monitor) PetscCall(PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", (double)objA, (double)(ngmres->gammaC * objmin)));
296: *selectRestart = PETSC_TRUE;
297: }
299: /* F_M stagnation restart */
300: if (ngmres->restart_fm_rise && objM > obj) {
301: if (ngmres->monitor) PetscCall(PetscViewerASCIIPrintf(ngmres->monitor, "F_M rise restart: %e > %e\n", (double)objM, (double)obj));
302: *selectRestart = PETSC_TRUE;
303: }
304: }
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }