Actual source code: ngmresfunc.c
petsc-3.6.4 2016-04-12
1: #include <../src/snes/impls/ngmres/snesngmres.h> /*I "petscsnes.h" I*/
2: #include <petscblaslapack.h>
6: PetscErrorCode SNESNGMRESUpdateSubspace_Private(SNES snes,PetscInt ivec,PetscInt l,Vec F,PetscReal fnorm,Vec X)
7: {
8: SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data;
9: Vec *Fdot = ngmres->Fdot;
10: Vec *Xdot = ngmres->Xdot;
14: if (ivec > l) SETERRQ2(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Cannot update vector %d with space size %d!",ivec,l);
15: VecCopy(F,Fdot[ivec]);
16: VecCopy(X,Xdot[ivec]);
18: ngmres->fnorms[ivec] = fnorm;
19: return(0);
20: }
24: PetscErrorCode SNESNGMRESFormCombinedSolution_Private(SNES snes,PetscInt ivec,PetscInt l,Vec XM,Vec FM,PetscReal fMnorm,Vec X,Vec XA,Vec FA)
25: {
26: SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data;
27: PetscInt i,j;
28: Vec *Fdot = ngmres->Fdot;
29: Vec *Xdot = ngmres->Xdot;
30: PetscScalar *beta = ngmres->beta;
31: PetscScalar *xi = ngmres->xi;
32: PetscScalar alph_total = 0.;
34: PetscReal nu;
35: Vec Y = snes->work[2];
36: PetscBool changed_y,changed_w;
39: nu = fMnorm*fMnorm;
41: /* construct the right hand side and xi factors */
42: if (l > 0) {
43: VecMDotBegin(FM,l,Fdot,xi);
44: VecMDotBegin(Fdot[ivec],l,Fdot,beta);
45: VecMDotEnd(FM,l,Fdot,xi);
46: VecMDotEnd(Fdot[ivec],l,Fdot,beta);
47: for (i = 0; i < l; i++) {
48: Q(i,ivec) = beta[i];
49: Q(ivec,i) = beta[i];
50: }
51: } else {
52: Q(0,0) = ngmres->fnorms[ivec]*ngmres->fnorms[ivec];
53: }
55: for (i = 0; i < l; i++) beta[i] = nu - xi[i];
57: /* construct h */
58: for (j = 0; j < l; j++) {
59: for (i = 0; i < l; i++) {
60: H(i,j) = Q(i,j)-xi[i]-xi[j]+nu;
61: }
62: }
63: if (l == 1) {
64: /* simply set alpha[0] = beta[0] / H[0, 0] */
65: if (H(0,0) != 0.) beta[0] = beta[0]/H(0,0);
66: else beta[0] = 0.;
67: } else {
68: #if defined(PETSC_MISSING_LAPACK_GELSS)
69: SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"NGMRES with LS requires the LAPACK GELSS routine.");
70: #else
71: PetscBLASIntCast(l,&ngmres->m);
72: PetscBLASIntCast(l,&ngmres->n);
73: ngmres->info = 0;
74: ngmres->rcond = -1.;
75: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
76: #if defined(PETSC_USE_COMPLEX)
77: 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));
78: #else
79: 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));
80: #endif
81: PetscFPTrapPop();
82: if (ngmres->info < 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"Bad argument to GELSS");
83: if (ngmres->info > 0) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD failed to converge");
84: #endif
85: }
86: for (i=0; i<l; i++) {
87: if (PetscIsInfOrNanScalar(beta[i])) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_LIB,"SVD generated inconsistent output");
88: }
89: alph_total = 0.;
90: for (i = 0; i < l; i++) alph_total += beta[i];
92: VecCopy(XM,XA);
93: VecScale(XA,1.-alph_total);
94: VecMAXPY(XA,l,beta,Xdot);
95: /* check the validity of the step */
96: VecCopy(XA,Y);
97: VecAXPY(Y,-1.0,X);
98: SNESLineSearchPostCheck(snes->linesearch,X,Y,XA,&changed_y,&changed_w);
99: if (!ngmres->approxfunc) {
100: if (snes->pc && snes->pcside == PC_LEFT) {
101: SNESApplyNPC(snes,XA,NULL,FA);
102: } else {
103: ierr =SNESComputeFunction(snes,XA,FA);
104: }
105: }
106: else {
107: VecCopy(FM,FA);
108: VecScale(FA,1.-alph_total);
109: VecMAXPY(FA,l,beta,Fdot);
110: }
111: return(0);
112: }
116: 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)
117: {
119: SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data;
120: PetscReal dcurnorm,dmin = -1.0;
121: Vec *Xdot = ngmres->Xdot;
122: PetscInt i;
125: if (xMnorm) {
126: VecNormBegin(XM,NORM_2,xMnorm);
127: }
128: if (fMnorm) {
129: VecNormBegin(FM,NORM_2,fMnorm);
130: }
131: if (yMnorm) {
132: VecCopy(X,D);
133: VecAXPY(D,-1.0,XM);
134: VecNormBegin(D,NORM_2,yMnorm);
135: }
136: if (xAnorm) {
137: VecNormBegin(XA,NORM_2,xAnorm);
138: }
139: if (fAnorm) {
140: VecNormBegin(FA,NORM_2,fAnorm);
141: }
142: if (yAnorm) {
143: VecCopy(X,D);
144: VecAXPY(D,-1.0,XA);
145: VecNormBegin(D,NORM_2,yAnorm);
146: }
147: if (dnorm) {
148: VecCopy(XA,D);
149: VecAXPY(D,-1.0,XM);
150: VecNormBegin(D,NORM_2,dnorm);
151: }
152: if (dminnorm) {
153: for (i=0; i<l; i++) {
154: VecCopy(Xdot[i],D);
155: VecAXPY(D,-1.0,XA);
156: VecNormBegin(D,NORM_2,&ngmres->xnorms[i]);
157: }
158: }
159: if (xMnorm) {VecNormEnd(XM,NORM_2,xMnorm);}
160: if (fMnorm) {VecNormEnd(FM,NORM_2,fMnorm);}
161: if (yMnorm) {VecNormEnd(D,NORM_2,yMnorm);}
162: if (xAnorm) {VecNormEnd(XA,NORM_2,xAnorm);}
163: if (fAnorm) {VecNormEnd(FA,NORM_2,fAnorm);}
164: if (yAnorm) {VecNormEnd(D,NORM_2,yAnorm);}
165: if (dnorm) {VecNormEnd(D,NORM_2,dnorm);}
166: if (dminnorm) {
167: for (i=0; i<l; i++) {
168: VecNormEnd(D,NORM_2,&ngmres->xnorms[i]);
169: dcurnorm = ngmres->xnorms[i];
170: if ((dcurnorm < dmin) || (dmin < 0.0)) dmin = dcurnorm;
171: }
172: *dminnorm = dmin;
173: }
174: return(0);
175: }
179: 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)
180: {
181: SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data;
182: PetscErrorCode ierr;
183: SNESLineSearchReason lssucceed;
184: PetscBool selectA;
187: if (ngmres->select_type == SNES_NGMRES_SELECT_LINESEARCH) {
188: /* X = X + \lambda(XA - X) */
189: if (ngmres->monitor) {
190: PetscViewerASCIIPrintf(ngmres->monitor,"||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);
191: }
192: VecCopy(FM,F);
193: VecCopy(XM,X);
194: VecCopy(XA,Y);
195: VecAYPX(Y,-1.0,X);
196: *fnorm = fMnorm;
197: SNESLineSearchApply(ngmres->additive_linesearch,X,F,fnorm,Y);
198: SNESLineSearchGetReason(ngmres->additive_linesearch,&lssucceed);
199: SNESLineSearchGetNorms(ngmres->additive_linesearch,xnorm,fnorm,ynorm);
200: if (lssucceed) {
201: if (++snes->numFailures >= snes->maxFailures) {
202: snes->reason = SNES_DIVERGED_LINE_SEARCH;
203: return(0);
204: }
205: }
206: if (ngmres->monitor) {
207: PetscViewerASCIIPrintf(ngmres->monitor,"Additive solution: ||F||_2 = %e\n",*fnorm);
208: }
209: } else if (ngmres->select_type == SNES_NGMRES_SELECT_DIFFERENCE) {
210: selectA = PETSC_TRUE;
211: /* Conditions for choosing the accelerated answer */
212: /* Criterion A -- the norm of the function isn't increased above the minimum by too much */
213: if (fAnorm >= ngmres->gammaA*fminnorm) selectA = PETSC_FALSE;
215: /* Criterion B -- the choice of x^A isn't too close to some other choice */
216: if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(*fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) {
217: } else selectA=PETSC_FALSE;
219: if (selectA) {
220: if (ngmres->monitor) {
221: PetscViewerASCIIPrintf(ngmres->monitor,"picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);
222: }
223: /* copy it over */
224: *xnorm = xAnorm;
225: *fnorm = fAnorm;
226: *ynorm = yAnorm;
227: VecCopy(FA,F);
228: VecCopy(XA,X);
229: } else {
230: if (ngmres->monitor) {
231: PetscViewerASCIIPrintf(ngmres->monitor,"picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n",fAnorm,fMnorm);
232: }
233: *xnorm = xMnorm;
234: *fnorm = fMnorm;
235: *ynorm = yMnorm;
236: VecCopy(XM,Y);
237: VecAXPY(Y,-1.0,X);
238: VecCopy(FM,F);
239: VecCopy(XM,X);
240: }
241: } else { /* none */
242: *xnorm = xAnorm;
243: *fnorm = fAnorm;
244: *ynorm = yAnorm;
245: VecCopy(FA,F);
246: VecCopy(XA,X);
247: }
248: return(0);
249: }
253: PetscErrorCode SNESNGMRESSelectRestart_Private(SNES snes,PetscInt l,PetscReal fAnorm,PetscReal dnorm,PetscReal fminnorm,PetscReal dminnorm,PetscBool *selectRestart)
254: {
255: SNES_NGMRES *ngmres = (SNES_NGMRES*)snes->data;
259: *selectRestart = PETSC_FALSE;
260: /* difference stagnation restart */
261: if ((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm)) && l > 0) {
262: if (ngmres->monitor) {
263: PetscViewerASCIIPrintf(ngmres->monitor,"difference restart: %e > %e\n",ngmres->epsilonB*dnorm,dminnorm);
264: }
265: *selectRestart = PETSC_TRUE;
266: }
267: /* residual stagnation restart */
268: if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) {
269: if (ngmres->monitor) {
270: PetscViewerASCIIPrintf(ngmres->monitor,"residual restart: %e > %e\n",PetscSqrtReal(fAnorm),ngmres->gammaC*PetscSqrtReal(fminnorm));
271: }
272: *selectRestart = PETSC_TRUE;
273: }
274: return(0);
275: }