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