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