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++) {
52: H(i,j) = Q(i,j)-xi[i]-xi[j]+nu;
53: }
54: }
55: if (l == 1) {
56: /* simply set alpha[0] = beta[0] / H[0, 0] */
57: if (H(0,0) != 0.) beta[0] = beta[0]/H(0,0);
58: else beta[0] = 0.;
59: } else {
60: PetscBLASIntCast(l,&ngmres->m);
61: PetscBLASIntCast(l,&ngmres->n);
62: ngmres->info = 0;
63: ngmres->rcond = -1.;
64: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
65: #if defined(PETSC_USE_COMPLEX)
66: PetscStackCallBLAS("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));
67: #else
68: PetscStackCallBLAS("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));
69: #endif
70: PetscFPTrapPop();
73: }
74: for (i=0; i<l; i++) {
76: }
77: alph_total = 0.;
78: for (i = 0; i < l; i++) alph_total += beta[i];
80: VecCopy(XM,XA);
81: VecScale(XA,1.-alph_total);
82: VecMAXPY(XA,l,beta,Xdot);
83: /* check the validity of the step */
84: VecCopy(XA,Y);
85: VecAXPY(Y,-1.0,X);
86: SNESLineSearchPostCheck(snes->linesearch,X,Y,XA,&changed_y,&changed_w);
87: if (!ngmres->approxfunc) {
88: if (snes->npc && snes->npcside== PC_LEFT) {
89: SNESApplyNPC(snes,XA,NULL,FA);
90: } else {
91: SNESComputeFunction(snes,XA,FA);
92: }
93: } else {
94: VecCopy(FM,FA);
95: VecScale(FA,1.-alph_total);
96: VecMAXPY(FA,l,beta,Fdot);
97: }
98: return 0;
99: }
101: 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)
102: {
103: SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data;
104: PetscReal dcurnorm,dmin = -1.0;
105: Vec *Xdot = ngmres->Xdot;
106: PetscInt i;
108: if (xMnorm) {
109: VecNormBegin(XM,NORM_2,xMnorm);
110: }
111: if (fMnorm) {
112: VecNormBegin(FM,NORM_2,fMnorm);
113: }
114: if (yMnorm) {
115: VecCopy(X,D);
116: VecAXPY(D,-1.0,XM);
117: VecNormBegin(D,NORM_2,yMnorm);
118: }
119: if (xAnorm) {
120: VecNormBegin(XA,NORM_2,xAnorm);
121: }
122: if (fAnorm) {
123: VecNormBegin(FA,NORM_2,fAnorm);
124: }
125: if (yAnorm) {
126: VecCopy(X,D);
127: VecAXPY(D,-1.0,XA);
128: VecNormBegin(D,NORM_2,yAnorm);
129: }
130: if (dnorm) {
131: VecCopy(XA,D);
132: VecAXPY(D,-1.0,XM);
133: VecNormBegin(D,NORM_2,dnorm);
134: }
135: if (dminnorm) {
136: for (i=0; i<l; i++) {
137: VecCopy(Xdot[i],D);
138: VecAXPY(D,-1.0,XA);
139: VecNormBegin(D,NORM_2,&ngmres->xnorms[i]);
140: }
141: }
142: if (xMnorm) VecNormEnd(XM,NORM_2,xMnorm);
143: if (fMnorm) VecNormEnd(FM,NORM_2,fMnorm);
144: if (yMnorm) VecNormEnd(D,NORM_2,yMnorm);
145: if (xAnorm) VecNormEnd(XA,NORM_2,xAnorm);
146: if (fAnorm) VecNormEnd(FA,NORM_2,fAnorm);
147: if (yAnorm) VecNormEnd(D,NORM_2,yAnorm);
148: if (dnorm) VecNormEnd(D,NORM_2,dnorm);
149: if (dminnorm) {
150: for (i=0; i<l; i++) {
151: VecNormEnd(D,NORM_2,&ngmres->xnorms[i]);
152: dcurnorm = ngmres->xnorms[i];
153: if ((dcurnorm < dmin) || (dmin < 0.0)) dmin = dcurnorm;
154: }
155: *dminnorm = dmin;
156: }
157: return 0;
158: }
160: 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)
161: {
162: SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data;
163: SNESLineSearchReason lssucceed;
164: PetscBool selectA;
166: if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) {
167: /* X = X + \lambda(XA - X) */
168: if (ngmres->monitor) {
169: PetscViewerASCIIPrintf(ngmres->monitor,"||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);
170: }
171: VecCopy(FM,F);
172: VecCopy(XM,X);
173: VecCopy(XA,Y);
174: VecAYPX(Y,-1.0,X);
175: *fnorm = fMnorm;
176: SNESLineSearchApply(ngmres->additive_linesearch,X,F,fnorm,Y);
177: SNESLineSearchGetReason(ngmres->additive_linesearch,&lssucceed);
178: SNESLineSearchGetNorms(ngmres->additive_linesearch,xnorm,fnorm,ynorm);
179: if (lssucceed) {
180: if (++snes->numFailures >= snes->maxFailures) {
181: snes->reason = SNES_DIVERGED_LINE_SEARCH;
182: return 0;
183: }
184: }
185: if (ngmres->monitor) {
186: PetscViewerASCIIPrintf(ngmres->monitor,"Additive solution: ||F||_2 = %e\n",*fnorm);
187: }
188: } else if (ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) {
189: selectA = PETSC_TRUE;
190: /* Conditions for choosing the accelerated answer */
191: /* Criterion A -- the norm of the function isn't increased above the minimum by too much */
192: if (fAnorm >= ngmres->gammaA*fminnorm) selectA = PETSC_FALSE;
194: /* Criterion B -- the choice of x^A isn't too close to some other choice */
195: if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(*fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) {
196: } else selectA=PETSC_FALSE;
198: if (selectA) {
199: if (ngmres->monitor) {
200: PetscViewerASCIIPrintf(ngmres->monitor,"picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);
201: }
202: /* copy it over */
203: *xnorm = xAnorm;
204: *fnorm = fAnorm;
205: *ynorm = yAnorm;
206: VecCopy(FA,F);
207: VecCopy(XA,X);
208: } else {
209: if (ngmres->monitor) {
210: PetscViewerASCIIPrintf(ngmres->monitor,"picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);
211: }
212: *xnorm = xMnorm;
213: *fnorm = fMnorm;
214: *ynorm = yMnorm;
215: VecCopy(XM,Y);
216: VecAXPY(Y,-1.0,X);
217: VecCopy(FM,F);
218: VecCopy(XM,X);
219: }
220: } else { /* none */
221: *xnorm = xAnorm;
222: *fnorm = fAnorm;
223: *ynorm = yAnorm;
224: VecCopy(FA,F);
225: VecCopy(XA,X);
226: }
227: return 0;
228: }
230: PetscErrorCode SNESNGMRESSelectRestart_Private(SNES snes,PetscInt l,PetscReal fMnorm, PetscReal fAnorm,PetscReal dnorm,PetscReal fminnorm,PetscReal dminnorm,PetscBool *selectRestart)
231: {
232: SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
234: *selectRestart = PETSC_FALSE;
235: /* difference stagnation restart */
236: if ((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm)) && l > 0) {
237: if (ngmres->monitor) {
238: PetscViewerASCIIPrintf(ngmres->monitor,"difference restart: %e > %e\n",ngmres->epsilonB*dnorm,dminnorm);
239: }
240: *selectRestart = PETSC_TRUE;
241: }
242: /* residual stagnation restart */
243: if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) {
244: if (ngmres->monitor) {
245: PetscViewerASCIIPrintf(ngmres->monitor,"residual restart: %e > %e\n",PetscSqrtReal(fAnorm),ngmres->gammaC*PetscSqrtReal(fminnorm));
246: }
247: *selectRestart = PETSC_TRUE;
248: }
250: /* F_M stagnation restart */
251: if (ngmres->restart_fm_rise && fMnorm > snes->norm) {
252: if (ngmres->monitor) {
253: PetscViewerASCIIPrintf(ngmres->monitor,"F_M rise restart: %e > %e\n",fMnorm,snes->norm);
254: }
255: *selectRestart = PETSC_TRUE;
256: }
258: return 0;
259: }