xref: /petsc/src/snes/impls/ngmres/snesngmres.c (revision 8cc86e312ab34c05aeaabff4b10649f02e511a2a)
1a312c225SMatthew G Knepley /* Defines the basic SNES object */
219653cdaSPeter Brune #include <../src/snes/impls/ngmres/snesngmres.h>
319653cdaSPeter Brune #include <petscblaslapack.h>
4a312c225SMatthew G Knepley 
5087dfb9eSxuemin 
6087dfb9eSxuemin 
7087dfb9eSxuemin 
8a312c225SMatthew G Knepley #undef __FUNCT__
9a312c225SMatthew G Knepley #define __FUNCT__ "SNESReset_NGMRES"
10a312c225SMatthew G Knepley PetscErrorCode SNESReset_NGMRES(SNES snes)
11a312c225SMatthew G Knepley {
12a312c225SMatthew G Knepley   SNES_NGMRES   *ngmres = (SNES_NGMRES*) snes->data;
13a312c225SMatthew G Knepley   PetscErrorCode ierr;
14a312c225SMatthew G Knepley 
15a312c225SMatthew G Knepley   PetscFunctionBegin;
16f109b39eSPeter Brune   ierr = VecDestroyVecs(ngmres->msize, &ngmres->Fdot);CHKERRQ(ierr);
17f109b39eSPeter Brune   ierr = VecDestroyVecs(ngmres->msize, &ngmres->Xdot);CHKERRQ(ierr);
18732c72c5SBarry Smith   if (snes->work) {ierr = VecDestroyVecs(snes->nwork, &snes->work);CHKERRQ(ierr);}
19a312c225SMatthew G Knepley   PetscFunctionReturn(0);
20a312c225SMatthew G Knepley }
21a312c225SMatthew G Knepley 
22a312c225SMatthew G Knepley #undef __FUNCT__
23a312c225SMatthew G Knepley #define __FUNCT__ "SNESDestroy_NGMRES"
24a312c225SMatthew G Knepley PetscErrorCode SNESDestroy_NGMRES(SNES snes)
25a312c225SMatthew G Knepley {
26a312c225SMatthew G Knepley   PetscErrorCode ierr;
27a312c225SMatthew G Knepley 
28a312c225SMatthew G Knepley   PetscFunctionBegin;
29a312c225SMatthew G Knepley   ierr = SNESReset_NGMRES(snes);CHKERRQ(ierr);
3019653cdaSPeter Brune   if (snes->data) {
3119653cdaSPeter Brune     SNES_NGMRES * ngmres = (SNES_NGMRES *)snes->data;
32f109b39eSPeter Brune     ierr = PetscFree5(ngmres->h, ngmres->beta, ngmres->xi, ngmres->fnorms, ngmres->q);CHKERRQ(ierr);
3319653cdaSPeter Brune     ierr = PetscFree(ngmres->s);CHKERRQ(ierr);
3419653cdaSPeter Brune #if PETSC_USE_COMPLEX
3519653cdaSPeter Brune     ierr = PetscFree(ngmres->rwork);
3619653cdaSPeter Brune #endif
3719653cdaSPeter Brune     ierr = PetscFree(ngmres->work);
3819653cdaSPeter Brune   }
3919653cdaSPeter Brune   ierr = PetscFree(snes->data);
40a312c225SMatthew G Knepley   PetscFunctionReturn(0);
41a312c225SMatthew G Knepley }
42a312c225SMatthew G Knepley 
43a312c225SMatthew G Knepley #undef __FUNCT__
44a312c225SMatthew G Knepley #define __FUNCT__ "SNESSetUp_NGMRES"
45a312c225SMatthew G Knepley PetscErrorCode SNESSetUp_NGMRES(SNES snes)
46a312c225SMatthew G Knepley {
47a312c225SMatthew G Knepley   SNES_NGMRES    *ngmres = (SNES_NGMRES *) snes->data;
4819653cdaSPeter Brune   PetscInt       msize,hsize;
49a312c225SMatthew G Knepley   PetscErrorCode ierr;
50a312c225SMatthew G Knepley 
51a312c225SMatthew G Knepley   PetscFunctionBegin;
52087dfb9eSxuemin   msize         = ngmres->msize;  /* restart size */
5319653cdaSPeter Brune   hsize         = msize * msize;
54087dfb9eSxuemin 
55087dfb9eSxuemin 
5698b3e84cSPeter Brune   /* explicit least squares minimization solve */
5719653cdaSPeter Brune   ierr = PetscMalloc5(hsize,PetscScalar,&ngmres->h,
5819653cdaSPeter Brune                       msize,PetscScalar,&ngmres->beta,
5919653cdaSPeter Brune                       msize,PetscScalar,&ngmres->xi,
60f109b39eSPeter Brune                       msize,PetscReal,  &ngmres->fnorms,
6119653cdaSPeter Brune                       hsize,PetscScalar,&ngmres->q);CHKERRQ(ierr);
6219653cdaSPeter Brune   ngmres->nrhs = 1;
6319653cdaSPeter Brune   ngmres->lda = msize;
6419653cdaSPeter Brune   ngmres->ldb = msize;
6519653cdaSPeter Brune   ierr = PetscMalloc(msize*sizeof(PetscScalar),&ngmres->s);CHKERRQ(ierr);
6619653cdaSPeter Brune   ierr = PetscMemzero(ngmres->h,   hsize*sizeof(PetscScalar));CHKERRQ(ierr);
6719653cdaSPeter Brune   ierr = PetscMemzero(ngmres->q,   hsize*sizeof(PetscScalar));CHKERRQ(ierr);
6819653cdaSPeter Brune   ierr = PetscMemzero(ngmres->xi,  msize*sizeof(PetscScalar));CHKERRQ(ierr);
6919653cdaSPeter Brune   ierr = PetscMemzero(ngmres->beta,msize*sizeof(PetscScalar));CHKERRQ(ierr);
7019653cdaSPeter Brune   ngmres->lwork = 12*msize;
7119653cdaSPeter Brune #if PETSC_USE_COMPLEX
7219653cdaSPeter Brune   ierr = PetscMalloc(sizeof(PetscReal)*ngmres->lwork,&ngmres->rwork);
7319653cdaSPeter Brune #endif
7419653cdaSPeter Brune   ierr = PetscMalloc(sizeof(PetscScalar)*ngmres->lwork,&ngmres->work);
75211b2d39SPeter Brune 
76f109b39eSPeter Brune   ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->Xdot);CHKERRQ(ierr);
77f109b39eSPeter Brune   ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->Fdot);CHKERRQ(ierr);
78f109b39eSPeter Brune   ierr = SNESDefaultGetWork(snes, 5);CHKERRQ(ierr);
79a312c225SMatthew G Knepley   PetscFunctionReturn(0);
80a312c225SMatthew G Knepley }
81a312c225SMatthew G Knepley 
82a312c225SMatthew G Knepley #undef __FUNCT__
83a312c225SMatthew G Knepley #define __FUNCT__ "SNESSetFromOptions_NGMRES"
84a312c225SMatthew G Knepley PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes)
85a312c225SMatthew G Knepley {
86a312c225SMatthew G Knepley   SNES_NGMRES   *ngmres = (SNES_NGMRES *) snes->data;
87a312c225SMatthew G Knepley   PetscErrorCode ierr;
88dfbf837cSBarry Smith   PetscBool      debug;
89a312c225SMatthew G Knepley 
90a312c225SMatthew G Knepley   PetscFunctionBegin;
91a312c225SMatthew G Knepley   ierr = PetscOptionsHead("SNES NGMRES options");CHKERRQ(ierr);
9219653cdaSPeter Brune   ierr = PetscOptionsInt("-snes_ngmres_m",         "Number of directions",               "SNES", ngmres->msize,  &ngmres->msize, PETSC_NULL);CHKERRQ(ierr);
9328ed4a04SPeter Brune   ierr = PetscOptionsInt("-snes_ngmres_restart_it","Tolerance iterations before restart","SNES", ngmres->restart_it,  &ngmres->restart_it, PETSC_NULL);CHKERRQ(ierr);
94dfbf837cSBarry Smith   ierr = PetscOptionsBool("-snes_ngmres_monitor",  "Monitor actions of NGMRES",          "SNES", ngmres->monitor ? PETSC_TRUE: PETSC_FALSE, &debug, PETSC_NULL);CHKERRQ(ierr);
95dfbf837cSBarry Smith   if (debug) {
96dfbf837cSBarry Smith     ngmres->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);CHKERRQ(ierr);
97dfbf837cSBarry Smith   }
986a7cf640SPeter Brune   ierr = PetscOptionsReal("-snes_ngmres_gammaA",   "Residual selection constant",   "SNES", ngmres->gammaA, &ngmres->gammaA, PETSC_NULL);CHKERRQ(ierr);
996a7cf640SPeter Brune   ierr = PetscOptionsReal("-snes_ngmres_gammaC", "  Residual restart constant",     "SNES", ngmres->gammaC, &ngmres->gammaC, PETSC_NULL);CHKERRQ(ierr);
1006a7cf640SPeter Brune   ierr = PetscOptionsReal("-snes_ngmres_epsilonB", "Difference selection constant", "SNES", ngmres->epsilonB, &ngmres->epsilonB, PETSC_NULL);CHKERRQ(ierr);
1016a7cf640SPeter Brune   ierr = PetscOptionsReal("-snes_ngmres_deltaB",   "Difference residual selection constant", "SNES", ngmres->deltaB, &ngmres->deltaB, PETSC_NULL);CHKERRQ(ierr);
102a312c225SMatthew G Knepley   ierr = PetscOptionsTail();CHKERRQ(ierr);
1036a7cf640SPeter Brune   if ((ngmres->gammaA > ngmres->gammaC) && (ngmres->gammaC > 2.)) ngmres->gammaC = ngmres->gammaA;
104a312c225SMatthew G Knepley   PetscFunctionReturn(0);
105a312c225SMatthew G Knepley }
106a312c225SMatthew G Knepley 
107a312c225SMatthew G Knepley #undef __FUNCT__
108a312c225SMatthew G Knepley #define __FUNCT__ "SNESView_NGMRES"
109a312c225SMatthew G Knepley PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer)
110a312c225SMatthew G Knepley {
111a312c225SMatthew G Knepley   SNES_NGMRES   *ngmres = (SNES_NGMRES *) snes->data;
112a312c225SMatthew G Knepley   PetscBool      iascii;
113f109b39eSPeter Brune   const char     *cstr;
114a312c225SMatthew G Knepley   PetscErrorCode ierr;
115a312c225SMatthew G Knepley 
116a312c225SMatthew G Knepley   PetscFunctionBegin;
117a312c225SMatthew G Knepley   ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
118a312c225SMatthew G Knepley   if (iascii) {
119f109b39eSPeter Brune     cstr = SNESLineSearchTypeName(snes->ls_type);
120f109b39eSPeter Brune     ierr = PetscViewerASCIIPrintf(viewer, "  line search variant: %s\n",cstr);CHKERRQ(ierr);
121f109b39eSPeter Brune     ierr = PetscViewerASCIIPrintf(viewer, "  Number of stored past updates: %d\n", ngmres->msize);CHKERRQ(ierr);
122f109b39eSPeter Brune     ierr = PetscViewerASCIIPrintf(viewer, "  Residual selection: gammaA=%1.0e, gammaC=%1.0e\n", ngmres->gammaA, ngmres->gammaC);CHKERRQ(ierr);
123f109b39eSPeter Brune     ierr = PetscViewerASCIIPrintf(viewer, "  Difference restart: epsilonB=%1.0e, deltaB=%1.0e\n", ngmres->epsilonB, ngmres->deltaB);CHKERRQ(ierr);
124a312c225SMatthew G Knepley   }
125a312c225SMatthew G Knepley   PetscFunctionReturn(0);
126a312c225SMatthew G Knepley }
127a312c225SMatthew G Knepley 
128f109b39eSPeter Brune EXTERN_C_BEGIN
129f109b39eSPeter Brune #undef __FUNCT__
130f109b39eSPeter Brune #define __FUNCT__ "SNESLineSearchSetType_NGMRES"
131f109b39eSPeter Brune PetscErrorCode  SNESLineSearchSetType_NGMRES(SNES snes, SNESLineSearchType type)
132f109b39eSPeter Brune {
133f109b39eSPeter Brune   PetscErrorCode ierr;
134f109b39eSPeter Brune   PetscFunctionBegin;
135f109b39eSPeter Brune 
136f109b39eSPeter Brune   switch (type) {
137f109b39eSPeter Brune   case SNES_LS_BASIC:
138f109b39eSPeter Brune     ierr = SNESLineSearchSet(snes,SNESLineSearchNo,PETSC_NULL);CHKERRQ(ierr);
139f109b39eSPeter Brune     break;
140f109b39eSPeter Brune   case SNES_LS_BASIC_NONORMS:
141f109b39eSPeter Brune     ierr = SNESLineSearchSet(snes,SNESLineSearchNoNorms,PETSC_NULL);CHKERRQ(ierr);
142f109b39eSPeter Brune     break;
143f109b39eSPeter Brune   case SNES_LS_QUADRATIC:
144f109b39eSPeter Brune     ierr = SNESLineSearchSet(snes,SNESLineSearchQuadraticSecant,PETSC_NULL);CHKERRQ(ierr);
145f109b39eSPeter Brune     break;
146f109b39eSPeter Brune   default:
147f109b39eSPeter Brune     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,"Unknown line search type.");
148f109b39eSPeter Brune     break;
149f109b39eSPeter Brune   }
150f109b39eSPeter Brune   snes->ls_type = type;
151f109b39eSPeter Brune   PetscFunctionReturn(0);
152f109b39eSPeter Brune }
153f109b39eSPeter Brune EXTERN_C_END
154f109b39eSPeter Brune 
15519653cdaSPeter Brune 
156a312c225SMatthew G Knepley #undef __FUNCT__
157a312c225SMatthew G Knepley #define __FUNCT__ "SNESSolve_NGMRES"
158211b2d39SPeter Brune 
159a312c225SMatthew G Knepley PetscErrorCode SNESSolve_NGMRES(SNES snes)
160a312c225SMatthew G Knepley {
161087dfb9eSxuemin   SNES_NGMRES        *ngmres = (SNES_NGMRES *) snes->data;
16298b3e84cSPeter Brune   /* present solution, residual, and preconditioned residual */
163f109b39eSPeter Brune   Vec                 X, F, B, D, G, W, Y;
164f109b39eSPeter Brune 
165f109b39eSPeter Brune   /* candidate linear combination answers */
166f109b39eSPeter Brune   Vec                 XA, FA;
16719653cdaSPeter Brune 
16898b3e84cSPeter Brune   /* previous iterations to construct the subspace */
169f109b39eSPeter Brune   Vec                 *Fdot = ngmres->Fdot;
170f109b39eSPeter Brune   Vec                 *Xdot = ngmres->Xdot;
17119653cdaSPeter Brune 
17298b3e84cSPeter Brune   /* coefficients and RHS to the minimization problem */
17319653cdaSPeter Brune   PetscScalar         *beta = ngmres->beta;
17419653cdaSPeter Brune   PetscScalar         *xi   = ngmres->xi;
175f109b39eSPeter Brune   PetscReal           fnorm, fAnorm, gnorm, ynorm, xnorm = 0.0;
17619653cdaSPeter Brune   PetscReal           nu;
17719653cdaSPeter Brune   PetscScalar         alph_total = 0.;
17819653cdaSPeter Brune   PetscScalar         qentry;
17928ed4a04SPeter Brune   PetscInt            i, j, k, k_restart, l, ivec, restart_count = 0;
18019653cdaSPeter Brune 
18198b3e84cSPeter Brune   /* solution selection data */
18219653cdaSPeter Brune   PetscBool           selectA, selectRestart;
183f109b39eSPeter Brune   PetscReal           dnorm, dminnorm, dcurnorm;
184f109b39eSPeter Brune   PetscReal           fminnorm;
18519653cdaSPeter Brune 
1861e633543SBarry Smith   SNESConvergedReason reason;
187f109b39eSPeter Brune   PetscBool           lssucceed;
188a312c225SMatthew G Knepley   PetscErrorCode      ierr;
189a312c225SMatthew G Knepley 
190a312c225SMatthew G Knepley   PetscFunctionBegin;
19198b3e84cSPeter Brune   /* variable initialization */
192a312c225SMatthew G Knepley   snes->reason  = SNES_CONVERGED_ITERATING;
193f109b39eSPeter Brune   X             = snes->vec_sol;
194f109b39eSPeter Brune   F             = snes->vec_func;
195f109b39eSPeter Brune   B             = snes->vec_rhs;
196f109b39eSPeter Brune   XA            = snes->vec_sol_update;
197f109b39eSPeter Brune   FA            = snes->work[0];
198f109b39eSPeter Brune   D             = snes->work[1];
199f109b39eSPeter Brune 
200f109b39eSPeter Brune   /* work for the line search */
201f109b39eSPeter Brune   Y             = snes->work[2];
202f109b39eSPeter Brune   G             = snes->work[3];
203f109b39eSPeter Brune   W             = snes->work[4];
204a312c225SMatthew G Knepley 
205a312c225SMatthew G Knepley   ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
206a312c225SMatthew G Knepley   snes->iter = 0;
207a312c225SMatthew G Knepley   snes->norm = 0.;
208a312c225SMatthew G Knepley   ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
20919653cdaSPeter Brune 
21098b3e84cSPeter Brune   /* initialization */
21119653cdaSPeter Brune 
21298b3e84cSPeter Brune   /* r = F(x) */
213f109b39eSPeter Brune   ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
214a312c225SMatthew G Knepley   if (snes->domainerror) {
215a312c225SMatthew G Knepley     snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
216a312c225SMatthew G Knepley     PetscFunctionReturn(0);
217a312c225SMatthew G Knepley   }
21819653cdaSPeter Brune 
21998b3e84cSPeter Brune   /* nu = (r, r) */
220f109b39eSPeter Brune   ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr);
221f109b39eSPeter Brune   fminnorm = fnorm;
222f109b39eSPeter Brune   nu = fnorm*fnorm;
223f109b39eSPeter Brune   if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in function evaluation");
22419653cdaSPeter Brune 
22598b3e84cSPeter Brune   /* q_{00} = nu  */
22619653cdaSPeter Brune   Q(0,0) = nu;
227f109b39eSPeter Brune   ngmres->fnorms[0] = fnorm;
228f109b39eSPeter Brune   /* Fdot[0] = F */
229f109b39eSPeter Brune   ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr);
230f109b39eSPeter Brune   ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr);
231087dfb9eSxuemin 
232a312c225SMatthew G Knepley   ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
233f109b39eSPeter Brune   snes->norm = fnorm;
234a312c225SMatthew G Knepley   ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
235f109b39eSPeter Brune   SNESLogConvHistory(snes, fnorm, 0);
236f109b39eSPeter Brune   ierr = SNESMonitor(snes, 0, fnorm);CHKERRQ(ierr);
237f109b39eSPeter Brune   ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
238a312c225SMatthew G Knepley   if (snes->reason) PetscFunctionReturn(0);
239a312c225SMatthew G Knepley 
24019653cdaSPeter Brune   k_restart = 1;
24119653cdaSPeter Brune   l = 1;
24219653cdaSPeter Brune   for (k=1; k<snes->max_its; k++) {
24398b3e84cSPeter Brune     /* select which vector of the stored subspace will be updated */
24498b3e84cSPeter Brune     ivec = k_restart % ngmres->msize; /* replace the last used part of the subspace */
24519653cdaSPeter Brune 
24698b3e84cSPeter Brune     /* Computation of x^M */
2478cc86e31SPeter Brune     if (snes->pc) {
2488cc86e31SPeter Brune       ierr = SNESSolve(snes->pc, B, X);CHKERRQ(ierr);
2498cc86e31SPeter Brune       ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr);
2508cc86e31SPeter Brune       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
2518cc86e31SPeter Brune         snes->reason = SNES_DIVERGED_INNER;
2528cc86e31SPeter Brune         PetscFunctionReturn(0);
2538cc86e31SPeter Brune       }
2548cc86e31SPeter Brune       ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
2558cc86e31SPeter Brune       ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr);
2568cc86e31SPeter Brune     } else if (snes->ops->computegs) {
2578cc86e31SPeter Brune       /* compute the update using the supplied Gauss-Seidel routine */
2588cc86e31SPeter Brune       ierr = SNESComputeGS(snes, B, X);CHKERRQ(ierr);
2598cc86e31SPeter Brune       ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr);
2608cc86e31SPeter Brune       ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr);
2618cc86e31SPeter Brune     } else {
262f109b39eSPeter Brune       /* no preconditioner -- just take gradient descent with line search */
263f109b39eSPeter Brune       ierr = VecCopy(F, Y);CHKERRQ(ierr);
264f109b39eSPeter Brune       ierr = VecScale(Y, -1.0);CHKERRQ(ierr);
265f109b39eSPeter Brune       ierr = (*snes->ops->linesearch)(snes,snes->lsP,X,F,Y,fnorm,xnorm,G,W,&ynorm,&gnorm,&lssucceed);CHKERRQ(ierr);
266f109b39eSPeter Brune       if (!lssucceed) {
267f109b39eSPeter Brune         if (++snes->numFailures >= snes->maxFailures) {
268f109b39eSPeter Brune           snes->reason = SNES_DIVERGED_LINE_SEARCH;
269f109b39eSPeter Brune           PetscFunctionReturn(0);
270f109b39eSPeter Brune         }
271f109b39eSPeter Brune       }
272f109b39eSPeter Brune       fnorm = gnorm;
273f109b39eSPeter Brune       ierr = VecCopy(G, F);CHKERRQ(ierr);
274f109b39eSPeter Brune       ierr = VecCopy(W, X);CHKERRQ(ierr);
2756634f59bSPeter Brune     }
2761e633543SBarry Smith 
27798b3e84cSPeter Brune     /* r = F(x) */
278f109b39eSPeter Brune     nu = fnorm*fnorm;
279f109b39eSPeter Brune     if (fminnorm > fnorm) fminnorm = fnorm;  /* the minimum norm is now of F^M */
28019653cdaSPeter Brune 
28198b3e84cSPeter Brune     /* construct the right hand side and xi factors */
28219653cdaSPeter Brune     for (i = 0; i < l; i++) {
283f109b39eSPeter Brune       ierr = VecDot(Fdot[i], F, &xi[i]);CHKERRQ(ierr);
28419653cdaSPeter Brune       beta[i] = nu - xi[i];
28519653cdaSPeter Brune     }
28619653cdaSPeter Brune 
28798b3e84cSPeter Brune     /* construct h */
28819653cdaSPeter Brune     for (j = 0; j < l; j++) {
28919653cdaSPeter Brune       for (i = 0; i < l; i++) {
29019653cdaSPeter Brune         H(i, j) = Q(i, j) - xi[i] - xi[j] + nu;
29119653cdaSPeter Brune       }
29219653cdaSPeter Brune     }
293f109b39eSPeter Brune 
294f109b39eSPeter Brune     if(l == 1) {
295f109b39eSPeter Brune       /* simply set alpha[0] = beta[0] / H[0, 0] */
296f109b39eSPeter Brune       beta[0] = beta[0] / H(0, 0);
297f109b39eSPeter Brune     } else {
29819653cdaSPeter Brune #ifdef PETSC_MISSING_LAPACK_GELSS
29919653cdaSPeter Brune     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "NGMRES with LS requires the LAPACK GELSS routine.");
3004a0c5b0cSMatthew G Knepley #else
30119653cdaSPeter Brune     ngmres->m = PetscBLASIntCast(l);
30219653cdaSPeter Brune     ngmres->n = PetscBLASIntCast(l);
30319653cdaSPeter Brune     ngmres->info = PetscBLASIntCast(0);
30419653cdaSPeter Brune     ngmres->rcond = -1.;
30519653cdaSPeter Brune #ifdef PETSC_USE_COMPLEX
30619653cdaSPeter Brune     LAPACKgelss_(&ngmres->m,
30719653cdaSPeter Brune                  &ngmres->n,
30819653cdaSPeter Brune                  &ngmres->nrhs,
30919653cdaSPeter Brune                  ngmres->h,
31019653cdaSPeter Brune                  &ngmres->lda,
31119653cdaSPeter Brune                  ngmres->beta,
31219653cdaSPeter Brune                  &ngmres->ldb,
31319653cdaSPeter Brune                  ngmres->s,
31419653cdaSPeter Brune                  &ngmres->rcond,
31519653cdaSPeter Brune                  &ngmres->rank,
31619653cdaSPeter Brune                  ngmres->work,
31719653cdaSPeter Brune                  &ngmres->lwork,
31819653cdaSPeter Brune                  ngmres->rwork,
31919653cdaSPeter Brune                  &ngmres->info);
32019653cdaSPeter Brune #else
32119653cdaSPeter Brune     LAPACKgelss_(&ngmres->m,
32219653cdaSPeter Brune                  &ngmres->n,
32319653cdaSPeter Brune                  &ngmres->nrhs,
32419653cdaSPeter Brune                  ngmres->h,
32519653cdaSPeter Brune                  &ngmres->lda,
32619653cdaSPeter Brune                  ngmres->beta,
32719653cdaSPeter Brune                  &ngmres->ldb,
32819653cdaSPeter Brune                  ngmres->s,
32919653cdaSPeter Brune                  &ngmres->rcond,
33019653cdaSPeter Brune                  &ngmres->rank,
33119653cdaSPeter Brune                  ngmres->work,
33219653cdaSPeter Brune                  &ngmres->lwork,
33319653cdaSPeter Brune                  &ngmres->info);
33419653cdaSPeter Brune #endif
33519653cdaSPeter Brune     if (ngmres->info < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELSS");
33619653cdaSPeter Brune     if (ngmres->info > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SVD failed to converge");
337a312c225SMatthew G Knepley #endif
338f109b39eSPeter Brune     }
339a312c225SMatthew G Knepley 
34019653cdaSPeter Brune     alph_total = 0.;
34119653cdaSPeter Brune     for (i = 0; i < l; i++) {
34219653cdaSPeter Brune       alph_total += beta[i];
343c0bbabecSxuemin     }
344f109b39eSPeter Brune 
345f109b39eSPeter Brune     ierr = VecCopy(X, XA);CHKERRQ(ierr);
346f109b39eSPeter Brune     ierr = VecScale(XA, 1. - alph_total);CHKERRQ(ierr);
347087dfb9eSxuemin 
34819653cdaSPeter Brune     for(i=0;i<l;i++){
349f109b39eSPeter Brune       ierr= VecAXPY(XA, beta[i], Xdot[i]);CHKERRQ(ierr);
35019653cdaSPeter Brune     }
351f109b39eSPeter Brune     ierr = SNESComputeFunction(snes, XA, FA);CHKERRQ(ierr);
352f109b39eSPeter Brune     ierr = VecNorm(FA, NORM_2, &fAnorm);CHKERRQ(ierr);
35319653cdaSPeter Brune 
35419653cdaSPeter Brune     selectA = PETSC_TRUE;
35598b3e84cSPeter Brune     /* Conditions for choosing the accelerated answer */
35619653cdaSPeter Brune 
35798b3e84cSPeter Brune     /* Criterion A -- the norm of the function isn't increased above the minimum by too much */
358f109b39eSPeter Brune     if (fAnorm >= ngmres->gammaA*fminnorm) {
35919653cdaSPeter Brune       selectA = PETSC_FALSE;
360a312c225SMatthew G Knepley     }
361087dfb9eSxuemin 
36298b3e84cSPeter Brune     /* Criterion B -- the choice of x^A isn't too close to some other choice */
363f109b39eSPeter Brune     ierr=VecCopy(XA, D);CHKERRQ(ierr);
364f109b39eSPeter Brune     ierr=VecAXPY(D, -1.0, X);CHKERRQ(ierr);
365f109b39eSPeter Brune     ierr=VecNorm(D, NORM_2, &dnorm);CHKERRQ(ierr);
366f109b39eSPeter Brune     dminnorm = -1.0;
36719653cdaSPeter Brune     for(i=0;i<l;i++) {
368f109b39eSPeter Brune       ierr=VecCopy(XA, D);CHKERRQ(ierr);
369f109b39eSPeter Brune       ierr=VecAXPY(D, -1.0, Xdot[i]);CHKERRQ(ierr);
370f109b39eSPeter Brune       ierr=VecNorm(D, NORM_2, &dcurnorm);CHKERRQ(ierr);
371f109b39eSPeter Brune       if((dcurnorm < dminnorm) || (dminnorm < 0.0)) dminnorm = dcurnorm;
37219653cdaSPeter Brune     }
373*cf22de31SPeter Brune     if (ngmres->epsilonB*dnorm<dminnorm || PetscSqrtReal(fnorm)<ngmres->deltaB*PetscSqrtReal(fminnorm)) {
37419653cdaSPeter Brune     } else {
37519653cdaSPeter Brune       selectA=PETSC_FALSE;
37619653cdaSPeter Brune     }
377211b2d39SPeter Brune 
378087dfb9eSxuemin 
37919653cdaSPeter Brune     if (selectA) {
380dfbf837cSBarry Smith       if (ngmres->monitor) {
381f109b39eSPeter Brune         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_A, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr);
382dfbf837cSBarry Smith       }
38398b3e84cSPeter Brune       /* copy it over */
384f109b39eSPeter Brune       fnorm = fAnorm;
385f109b39eSPeter Brune       nu = fnorm*fnorm;
386f109b39eSPeter Brune       ierr = VecCopy(FA, F);CHKERRQ(ierr);
387f109b39eSPeter Brune       ierr = VecCopy(XA, X);CHKERRQ(ierr);
38819653cdaSPeter Brune     } else {
389dfbf837cSBarry Smith       if (ngmres->monitor) {
390f109b39eSPeter Brune         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "picked X_M, ||F_A||_2 = %e, ||F_M||_2 = %e\n", fAnorm, fnorm);CHKERRQ(ierr);
391dfbf837cSBarry Smith       }
39219653cdaSPeter Brune     }
393211b2d39SPeter Brune 
39419653cdaSPeter Brune     selectRestart = PETSC_FALSE;
39519653cdaSPeter Brune 
39698b3e84cSPeter Brune     /* difference stagnation restart */
397*cf22de31SPeter Brune     if((ngmres->epsilonB*dnorm > dminnorm) && (PetscSqrtReal(fAnorm) > ngmres->deltaB*PetscSqrtReal(fminnorm))) {
398dfbf837cSBarry Smith       if (ngmres->monitor) {
399f109b39eSPeter Brune         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "difference restart: %e > %e\n", ngmres->epsilonB*dnorm, dminnorm);CHKERRQ(ierr);
400dfbf837cSBarry Smith       }
40119653cdaSPeter Brune       selectRestart = PETSC_TRUE;
40219653cdaSPeter Brune     }
40398b3e84cSPeter Brune     /* residual stagnation restart */
404*cf22de31SPeter Brune     if (PetscSqrtReal(fAnorm) > ngmres->gammaC*PetscSqrtReal(fminnorm)) {
405dfbf837cSBarry Smith       if (ngmres->monitor) {
406*cf22de31SPeter Brune         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "residual restart: %e > %e\n", PetscSqrtReal(fAnorm), ngmres->gammaC*PetscSqrtReal(fminnorm));CHKERRQ(ierr);
407dfbf837cSBarry Smith       }
40819653cdaSPeter Brune       selectRestart = PETSC_TRUE;
40919653cdaSPeter Brune     }
41019653cdaSPeter Brune 
41128ed4a04SPeter Brune     /* if the restart conditions persist for more than restart_it iterations, restart. */
41219653cdaSPeter Brune     if (selectRestart) {
41328ed4a04SPeter Brune       restart_count++;
41428ed4a04SPeter Brune     } else {
41528ed4a04SPeter Brune       restart_count = 0;
41628ed4a04SPeter Brune     }
41728ed4a04SPeter Brune 
41828ed4a04SPeter Brune     /* restart after restart conditions have persisted for a fixed number of iterations */
41928ed4a04SPeter Brune     if (restart_count >= ngmres->restart_it) {
420dfbf837cSBarry Smith       if (ngmres->monitor){
421dfbf837cSBarry Smith         ierr = PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %d\n", k_restart);CHKERRQ(ierr);
422dfbf837cSBarry Smith       }
42328ed4a04SPeter Brune       restart_count = 0;
42419653cdaSPeter Brune       k_restart = 1;
42519653cdaSPeter Brune       l = 1;
42698b3e84cSPeter Brune       /* q_{00} = nu */
427f109b39eSPeter Brune       ngmres->fnorms[0] = fnorm;
428f109b39eSPeter Brune       nu = fnorm*fnorm;
42919653cdaSPeter Brune       Q(0,0) = nu;
430f109b39eSPeter Brune       /* Fdot[0] = F */
431f109b39eSPeter Brune       ierr = VecCopy(X, Xdot[0]);CHKERRQ(ierr);
432f109b39eSPeter Brune       ierr = VecCopy(F, Fdot[0]);CHKERRQ(ierr);
43319653cdaSPeter Brune     } else {
43498b3e84cSPeter Brune       /* select the current size of the subspace */
4351e633543SBarry Smith       if (l < ngmres->msize) l++;
43619653cdaSPeter Brune       k_restart++;
43798b3e84cSPeter Brune       /* place the current entry in the list of previous entries */
438f109b39eSPeter Brune       ierr = VecCopy(F, Fdot[ivec]);CHKERRQ(ierr);
439f109b39eSPeter Brune       ierr = VecCopy(X, Xdot[ivec]);CHKERRQ(ierr);
440f109b39eSPeter Brune       ngmres->fnorms[ivec] = fnorm;
441f109b39eSPeter Brune       if (fminnorm > fnorm) fminnorm = fnorm;  /* the minimum norm is now of r^A */
44219653cdaSPeter Brune       for (i = 0; i < l; i++) {
443f109b39eSPeter Brune         ierr = VecDot(F, Fdot[i], &qentry);CHKERRQ(ierr);
44419653cdaSPeter Brune         Q(i, ivec) = qentry;
44519653cdaSPeter Brune         Q(ivec, i) = qentry;
44619653cdaSPeter Brune       }
44719653cdaSPeter Brune     }
44819653cdaSPeter Brune 
4498409ca45SMatthew G Knepley     ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr);
450087dfb9eSxuemin     snes->iter = k;
451f109b39eSPeter Brune     snes->norm = fnorm;
4528409ca45SMatthew G Knepley     ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr);
4538409ca45SMatthew G Knepley     SNESLogConvHistory(snes, snes->norm, snes->iter);
4548409ca45SMatthew G Knepley     ierr = SNESMonitor(snes, snes->iter, snes->norm);CHKERRQ(ierr);
4558409ca45SMatthew G Knepley 
456f109b39eSPeter Brune     ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr);
457087dfb9eSxuemin     if (snes->reason) PetscFunctionReturn(0);
458a312c225SMatthew G Knepley   }
459a312c225SMatthew G Knepley   snes->reason = SNES_DIVERGED_MAX_IT;
460a312c225SMatthew G Knepley   PetscFunctionReturn(0);
461a312c225SMatthew G Knepley }
462a312c225SMatthew G Knepley 
463dfbf837cSBarry Smith /*MC
464dfbf837cSBarry Smith   SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio.
465a312c225SMatthew G Knepley 
466dfbf837cSBarry Smith    Level: beginner
467dfbf837cSBarry Smith 
468dfbf837cSBarry Smith    "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio,
469dfbf837cSBarry Smith    SIAM Journal on Scientific Computing, 21(5), 2000.
470dfbf837cSBarry Smith 
471dfbf837cSBarry Smith    This is also the same as the algorithm called Anderson acceleration introduced in "D. G. Anderson. Iterative procedures for nonlinear integral equations.
472dfbf837cSBarry Smith    J. Assoc. Comput. Mach., 12:547–560, 1965."
473dfbf837cSBarry Smith 
474dfbf837cSBarry Smith .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types)
475dfbf837cSBarry Smith M*/
476a312c225SMatthew G Knepley 
477a312c225SMatthew G Knepley EXTERN_C_BEGIN
478a312c225SMatthew G Knepley #undef __FUNCT__
479a312c225SMatthew G Knepley #define __FUNCT__ "SNESCreate_NGMRES"
480a312c225SMatthew G Knepley PetscErrorCode SNESCreate_NGMRES(SNES snes)
481a312c225SMatthew G Knepley {
482a312c225SMatthew G Knepley   SNES_NGMRES   *ngmres;
483a312c225SMatthew G Knepley   PetscErrorCode ierr;
484a312c225SMatthew G Knepley 
485a312c225SMatthew G Knepley   PetscFunctionBegin;
486a312c225SMatthew G Knepley   snes->ops->destroy        = SNESDestroy_NGMRES;
487a312c225SMatthew G Knepley   snes->ops->setup          = SNESSetUp_NGMRES;
488a312c225SMatthew G Knepley   snes->ops->setfromoptions = SNESSetFromOptions_NGMRES;
489a312c225SMatthew G Knepley   snes->ops->view           = SNESView_NGMRES;
490a312c225SMatthew G Knepley   snes->ops->solve          = SNESSolve_NGMRES;
491a312c225SMatthew G Knepley   snes->ops->reset          = SNESReset_NGMRES;
492a312c225SMatthew G Knepley 
49342f4f86dSBarry Smith   snes->usespc          = PETSC_TRUE;
4942c155ee1SBarry Smith   snes->usesksp         = PETSC_FALSE;
4952c155ee1SBarry Smith 
496a312c225SMatthew G Knepley   ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr);
497a312c225SMatthew G Knepley   snes->data = (void*) ngmres;
49819653cdaSPeter Brune   ngmres->msize = 10;
49919653cdaSPeter Brune 
50028ed4a04SPeter Brune   ngmres->restart_it = 2;
501f109b39eSPeter Brune   ngmres->gammaA     = 2.0;
502f109b39eSPeter Brune   ngmres->gammaC     = 2.0;
503cac108bcSPeter Brune   ngmres->deltaB     = 0.9;
504cac108bcSPeter Brune   ngmres->epsilonB   = 0.1;
505f109b39eSPeter Brune   ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESLineSearchSetType_C","SNESLineSearchSetType_NGMRES",SNESLineSearchSetType_NGMRES);CHKERRQ(ierr);
506f109b39eSPeter Brune   ierr = SNESLineSearchSetType(snes, SNES_LS_QUADRATIC);CHKERRQ(ierr);
507a312c225SMatthew G Knepley   PetscFunctionReturn(0);
508a312c225SMatthew G Knepley }
509a312c225SMatthew G Knepley EXTERN_C_END
510