Actual source code: ngmresfunc.c
1: #include <../src/snes/impls/ngmres/snesngmres.h>
2: #include <petscblaslapack.h>
4: PetscErrorCode SNESNGMRESUpdateSubspace_Private(SNES snes, PetscInt ivec, PetscInt l, Vec F, PetscReal fnorm, Vec X)
5: {
6: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
7: Vec *Fdot = ngmres->Fdot;
8: Vec *Xdot = ngmres->Xdot;
11: VecCopy(F, Fdot[ivec]);
12: VecCopy(X, Xdot[ivec]);
14: ngmres->fnorms[ivec] = fnorm;
15: return 0;
16: }
18: PetscErrorCode SNESNGMRESFormCombinedSolution_Private(SNES snes, PetscInt ivec, PetscInt l, Vec XM, Vec FM, PetscReal fMnorm, Vec X, Vec XA, Vec FA)
19: {
20: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
21: PetscInt i, j;
22: Vec *Fdot = ngmres->Fdot;
23: Vec *Xdot = ngmres->Xdot;
24: PetscScalar *beta = ngmres->beta;
25: PetscScalar *xi = ngmres->xi;
26: PetscScalar alph_total = 0.;
27: PetscReal nu;
28: Vec Y = snes->work[2];
29: PetscBool changed_y, changed_w;
31: nu = fMnorm * fMnorm;
33: /* construct the right hand side and xi factors */
34: if (l > 0) {
35: VecMDotBegin(FM, l, Fdot, xi);
36: VecMDotBegin(Fdot[ivec], l, Fdot, beta);
37: VecMDotEnd(FM, l, Fdot, xi);
38: VecMDotEnd(Fdot[ivec], l, Fdot, beta);
39: for (i = 0; i < l; i++) {
40: Q(i, ivec) = beta[i];
41: Q(ivec, i) = beta[i];
42: }
43: } else {
44: Q(0, 0) = ngmres->fnorms[ivec] * ngmres->fnorms[ivec];
45: }
47: for (i = 0; i < l; i++) beta[i] = nu - xi[i];
49: /* construct h */
50: for (j = 0; j < l; j++) {
51: for (i = 0; i < l; i++) H(i, j) = Q(i, j) - xi[i] - xi[j] + nu;
52: }
53: if (l == 1) {
54: /* simply set alpha[0] = beta[0] / H[0, 0] */
55: if (H(0, 0) != 0.) beta[0] = beta[0] / H(0, 0);
56: else beta[0] = 0.;
57: } else {
58: PetscBLASIntCast(l, &ngmres->m);
59: PetscBLASIntCast(l, &ngmres->n);
60: ngmres->info = 0;
61: ngmres->rcond = -1.;
62: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
63: #if defined(PETSC_USE_COMPLEX)
64: PetscCallBLAS("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, &ngmres->info));
65: #else
66: PetscCallBLAS("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->info));
67: #endif
68: PetscFPTrapPop();
71: }
73: alph_total = 0.;
74: for (i = 0; i < l; i++) alph_total += beta[i];
76: VecCopy(XM, XA);
77: VecScale(XA, 1. - alph_total);
78: VecMAXPY(XA, l, beta, Xdot);
79: /* check the validity of the step */
80: VecCopy(XA, Y);
81: VecAXPY(Y, -1.0, X);
82: SNESLineSearchPostCheck(snes->linesearch, X, Y, XA, &changed_y, &changed_w);
83: if (!ngmres->approxfunc) {
84: if (snes->npc && snes->npcside == PC_LEFT) {
85: SNESApplyNPC(snes, XA, NULL, FA);
86: } else {
87: SNESComputeFunction(snes, XA, FA);
88: }
89: } else {
90: VecCopy(FM, FA);
91: VecScale(FA, 1. - alph_total);
92: VecMAXPY(FA, l, beta, Fdot);
93: }
94: return 0;
95: }
97: 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)
98: {
99: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
100: PetscReal dcurnorm, dmin = -1.0;
101: Vec *Xdot = ngmres->Xdot;
102: PetscInt i;
104: if (xMnorm) VecNormBegin(XM, NORM_2, xMnorm);
105: if (fMnorm) VecNormBegin(FM, NORM_2, fMnorm);
106: if (yMnorm) {
107: VecCopy(X, D);
108: VecAXPY(D, -1.0, XM);
109: VecNormBegin(D, NORM_2, yMnorm);
110: }
111: if (xAnorm) VecNormBegin(XA, NORM_2, xAnorm);
112: if (fAnorm) VecNormBegin(FA, NORM_2, fAnorm);
113: if (yAnorm) {
114: VecCopy(X, D);
115: VecAXPY(D, -1.0, XA);
116: VecNormBegin(D, NORM_2, yAnorm);
117: }
118: if (dnorm) {
119: VecCopy(XA, D);
120: VecAXPY(D, -1.0, XM);
121: VecNormBegin(D, NORM_2, dnorm);
122: }
123: if (dminnorm) {
124: for (i = 0; i < l; i++) {
125: VecCopy(Xdot[i], D);
126: VecAXPY(D, -1.0, XA);
127: VecNormBegin(D, NORM_2, &ngmres->xnorms[i]);
128: }
129: }
130: if (xMnorm) VecNormEnd(XM, NORM_2, xMnorm);
131: if (fMnorm) VecNormEnd(FM, NORM_2, fMnorm);
132: if (yMnorm) VecNormEnd(D, NORM_2, yMnorm);
133: if (xAnorm) VecNormEnd(XA, NORM_2, xAnorm);
134: if (fAnorm) VecNormEnd(FA, NORM_2, fAnorm);
135: if (yAnorm) VecNormEnd(D, NORM_2, yAnorm);
136: if (dnorm) VecNormEnd(D, NORM_2, dnorm);
137: if (dminnorm) {
138: for (i = 0; i < l; i++) {
139: VecNormEnd(D, NORM_2, &ngmres->xnorms[i]);
140: dcurnorm = ngmres->xnorms[i];
141: if ((dcurnorm < dmin) || (dmin < 0.0)) dmin = dcurnorm;
142: }
143: *dminnorm = dmin;
144: }
145: return 0;
146: }
148: PetscErrorCode SNESNGMRESSelect_Private(SNES snes, PetscInt k_restart, Vec XM, Vec FM, PetscReal xMnorm, PetscReal fMnorm, PetscReal yMnorm, Vec XA, Vec FA, PetscReal xAnorm, PetscReal fAnorm, PetscReal yAnorm, PetscReal dnorm, PetscReal fminnorm, PetscReal dminnorm, Vec X, Vec F, Vec Y, PetscReal *xnorm, PetscReal *fnorm, PetscReal *ynorm)
149: {
150: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
151: SNESLineSearchReason lssucceed;
152: PetscBool selectA;
154: if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) {
155: /* X = X + \lambda(XA - X) */
156: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "||F_A||_2 = %e, ||F_M||_2 = %e\n", (double)fAnorm, (double)fMnorm);
157: VecCopy(FM, F);
158: VecCopy(XM, X);
159: VecCopy(XA, Y);
160: VecAYPX(Y, -1.0, X);
161: *fnorm = fMnorm;
162: SNESLineSearchApply(ngmres->additive_linesearch, X, F, fnorm, Y);
163: SNESLineSearchGetReason(ngmres->additive_linesearch, &lssucceed);
164: SNESLineSearchGetNorms(ngmres->additive_linesearch, xnorm, fnorm, ynorm);
165: if (lssucceed) {
166: if (++snes->numFailures >= snes->maxFailures) {
167: snes->reason = SNES_DIVERGED_LINE_SEARCH;
168: return 0;
169: }
170: }
171: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "Additive solution: ||F||_2 = %e\n", (double)*fnorm);
172: } else if (ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) {
173: selectA = PETSC_TRUE;
174: /* Conditions for choosing the accelerated answer */
175: /* Criterion A -- the norm of the function isn't increased above the minimum by too much */
176: if (fAnorm >= ngmres->gammaA * fminnorm) selectA = PETSC_FALSE;
178: /* Criterion B -- the choice of x^A isn't too close to some other choice */
179: if (ngmres->epsilonB * dnorm < dminnorm || PetscSqrtReal(*fnorm) < ngmres->deltaB * PetscSqrtReal(fminnorm)) {
180: } else selectA = PETSC_FALSE;
182: if (selectA) {
183: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n", (double)fAnorm, (double)fMnorm);
184: /* copy it over */
185: *xnorm = xAnorm;
186: *fnorm = fAnorm;
187: *ynorm = yAnorm;
188: VecCopy(FA, F);
189: VecCopy(XA, X);
190: } else {
191: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n", (double)fAnorm, (double)fMnorm);
192: *xnorm = xMnorm;
193: *fnorm = fMnorm;
194: *ynorm = yMnorm;
195: VecCopy(XM, Y);
196: VecAXPY(Y, -1.0, X);
197: VecCopy(FM, F);
198: VecCopy(XM, X);
199: }
200: } else { /* none */
201: *xnorm = xAnorm;
202: *fnorm = fAnorm;
203: *ynorm = yAnorm;
204: VecCopy(FA, F);
205: VecCopy(XA, X);
206: }
207: return 0;
208: }
210: PetscErrorCode SNESNGMRESSelectRestart_Private(SNES snes, PetscInt l, PetscReal fMnorm, PetscReal fAnorm, PetscReal dnorm, PetscReal fminnorm, PetscReal dminnorm, PetscBool *selectRestart)
211: {
212: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
214: *selectRestart = PETSC_FALSE;
215: /* difference stagnation restart */
216: if ((ngmres->epsilonB * dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB * PetscSqrtReal(fminnorm)) && l > 0) {
217: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", (double)(ngmres->epsilonB * dnorm), (double)dminnorm);
218: *selectRestart = PETSC_TRUE;
219: }
220: /* residual stagnation restart */
221: if (PetscSqrtReal(fAnorm) > ngmres->gammaC * PetscSqrtReal(fminnorm)) {
222: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", (double)PetscSqrtReal(fAnorm), (double)(ngmres->gammaC * PetscSqrtReal(fminnorm)));
223: *selectRestart = PETSC_TRUE;
224: }
226: /* F_M stagnation restart */
227: if (ngmres->restart_fm_rise && fMnorm > snes->norm) {
228: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "F_M rise restart: %e > %e\n", (double)fMnorm, (double)snes->norm);
229: *selectRestart = PETSC_TRUE;
230: }
232: return 0;
233: }