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: }