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