1 /* Defines the basic SNES object */ 2 #include <private/snesimpl.h> 3 4 /* Private structure for the Anderson mixing method aka nonlinear Krylov */ 5 typedef struct { 6 Vec *v, *w; 7 PetscReal *f2; /* 2-norms of function (residual) at each stage */ 8 PetscInt msize; /* maximum size of space */ 9 PetscInt csize; /* current size of space */ 10 } SNES_NGMRES; 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "SNESReset_NGMRES" 14 PetscErrorCode SNESReset_NGMRES(SNES snes) 15 { 16 SNES_NGMRES *ngmres = (SNES_NGMRES*) snes->data; 17 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 ierr = VecDestroyVecs(ngmres->msize, &ngmres->v);CHKERRQ(ierr); 21 ierr = VecDestroyVecs(ngmres->msize, &ngmres->w);CHKERRQ(ierr); 22 PetscFunctionReturn(0); 23 } 24 25 #undef __FUNCT__ 26 #define __FUNCT__ "SNESDestroy_NGMRES" 27 PetscErrorCode SNESDestroy_NGMRES(SNES snes) 28 { 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 ierr = SNESReset_NGMRES(snes);CHKERRQ(ierr); 33 if (snes->work) {ierr = VecDestroyVecs(snes->nwork, &snes->work);CHKERRQ(ierr);} 34 PetscFunctionReturn(0); 35 } 36 37 #undef __FUNCT__ 38 #define __FUNCT__ "SNESSetUp_NGMRES" 39 PetscErrorCode SNESSetUp_NGMRES(SNES snes) 40 { 41 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 42 PetscErrorCode ierr; 43 44 PetscFunctionBegin; 45 #if 0 46 if (snes->pc_side != PC_LEFT) {SETERRQ(((PetscObject) snes)->comm, PETSC_ERR_SUP, "Only left preconditioning allowed for SNESNGMRES");} 47 #endif 48 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->v);CHKERRQ(ierr); 49 ierr = VecDuplicateVecs(snes->vec_sol, ngmres->msize, &ngmres->w);CHKERRQ(ierr); 50 ierr = SNESDefaultGetWork(snes, 1);CHKERRQ(ierr); 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNCT__ 55 #define __FUNCT__ "SNESSetFromOptions_NGMRES" 56 PetscErrorCode SNESSetFromOptions_NGMRES(SNES snes) 57 { 58 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 59 PetscErrorCode ierr; 60 61 PetscFunctionBegin; 62 ierr = PetscOptionsHead("SNES NGMRES options");CHKERRQ(ierr); 63 ierr = PetscOptionsInt("-snes_gmres_restart", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, PETSC_NULL);CHKERRQ(ierr); 64 ierr = PetscOptionsTail();CHKERRQ(ierr); 65 PetscFunctionReturn(0); 66 } 67 68 #undef __FUNCT__ 69 #define __FUNCT__ "SNESView_NGMRES" 70 PetscErrorCode SNESView_NGMRES(SNES snes, PetscViewer viewer) 71 { 72 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 73 PetscBool iascii; 74 PetscErrorCode ierr; 75 76 PetscFunctionBegin; 77 ierr = PetscTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 78 if (iascii) { 79 ierr = PetscViewerASCIIPrintf(viewer, " Size of space %d\n", ngmres->msize);CHKERRQ(ierr); 80 } 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "SNESSolve_NGMRES" 86 PetscErrorCode SNESSolve_NGMRES(SNES snes) 87 { 88 SNES_NGMRES *ngmres = (SNES_NGMRES *) snes->data; 89 SNES pc; 90 Vec X, F, r, rOld, *V = ngmres->v, *W = ngmres->w; 91 PetscScalar wdot; 92 PetscReal fnorm; 93 PetscInt i, j, k; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 snes->reason = SNES_CONVERGED_ITERATING; 98 X = snes->vec_sol; 99 F = snes->vec_func; 100 rOld = snes->work[0]; 101 102 ierr = SNESGetPC(snes, &pc);CHKERRQ(ierr); 103 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 104 snes->iter = 0; 105 snes->norm = 0.; 106 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 107 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr); 108 if (snes->domainerror) { 109 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 110 PetscFunctionReturn(0); 111 } 112 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 113 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in norm"); 114 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 115 snes->norm = fnorm; 116 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 117 SNESLogConvHistory(snes, fnorm, 0); 118 ierr = SNESMonitor(snes, 0, fnorm);CHKERRQ(ierr); 119 120 /* set parameter for default relative tolerance convergence test */ 121 snes->ttol = fnorm*snes->rtol; 122 /* test convergence */ 123 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 124 if (snes->reason) PetscFunctionReturn(0); 125 126 #ifndef OLD_BARRY_CODE 127 ierr = VecCopy(F, rOld);CHKERRQ(ierr); 128 #else 129 /* Barry: What the heck is this part doing? */ 130 /* determine optimal scale factor -- slow code */ 131 ierr = VecDuplicate(P, &y);CHKERRQ(ierr); 132 ierr = VecDuplicate(P, &w);CHKERRQ(ierr); 133 ierr = MatMult(Amat, Pold, y);CHKERRQ(ierr); 134 /*ierr = KSP_PCApplyBAorAB(ksp,Pold,y,w);CHKERRQ(ierr); */ /* y = BAp */ 135 ierr = VecDotNorm2(Pold, y, &rdot, &abr);CHKERRQ(ierr); /* rdot = (p)^T(BAp); abr = (BAp)^T (BAp) */ 136 ierr = VecDestroy(&y);CHKERRQ(ierr); 137 ierr = VecDestroy(&w);CHKERRQ(ierr); 138 A0 = rdot/abr; 139 ierr = VecAXPY(X,A0,Pold);CHKERRQ(ierr); /* x <- x + scale p */ 140 #endif 141 142 /* Loop over batches of directions */ 143 /* Code from Barry I do not understand to solve the least-squares problem, Time to try again */ 144 for(k = 0; k < snes->max_its; k += ngmres->msize) { 145 /* Loop over updates for this batch */ 146 /* TODO: Incorporate the variant which use the analytic Jacobian */ 147 /* TODO: Incorporate selection criteria for u^A from paper (need to store residual and update norms) */ 148 /* TODO: Incorporate criteria for restarting from paper */ 149 for(i = 0; i < ngmres->msize && k+i < snes->max_its; ++i) { 150 /* Get a new approxiate solution u^k */ 151 ierr = SNESSolve(pc, PETSC_NULL, X);CHKERRQ(ierr); 152 /* Solve least squares problem using last msize iterates */ 153 ierr = SNESGetFunction(pc, &r, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); /* r = F(u^M) */ 154 for(j = 0; j < i; ++j) { /* r = \prod_j (I + v_j w^T_j) r */ 155 ierr = VecDot(W[j], r, &wdot);CHKERRQ(ierr); 156 ierr = VecAXPY(r, wdot, V[j]);CHKERRQ(ierr); 157 } 158 ierr = VecCopy(rOld, W[i]);CHKERRQ(ierr); /* w_i = r_{old} */ 159 ierr = VecAXPY(rOld, -1.0, r);CHKERRQ(ierr); /* v_i = r */ 160 ierr = VecDot(W[i], rOld, &wdot);CHKERRQ(ierr); /* ------------------- */ 161 ierr = VecCopy(r, V[i]);CHKERRQ(ierr); /* w^T_i (r_{old} - r) */ 162 ierr = VecScale(V[i], 1.0/wdot);CHKERRQ(ierr); 163 ierr = VecDot(W[i], r, &wdot);CHKERRQ(ierr); /* r = (I + v_i w^T_i) r */ 164 ierr = VecAXPY(r, wdot, V[i]);CHKERRQ(ierr); 165 ierr = VecCopy(r, rOld);CHKERRQ(ierr); /* r_{old} = r */ 166 167 ierr = VecAXPY(X, 1.0, r);CHKERRQ(ierr); /* u^A = u^M + r */ 168 /* Monitor and log history */ 169 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr); /* F(u^k) */ 170 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 171 SNESLogConvHistory(snes, fnorm, 0); 172 ierr = SNESMonitor(snes, i+k+1, fnorm);CHKERRQ(ierr); 173 /* Check convergence */ 174 ierr = (*snes->ops->converged)(snes, i+k+1, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);CHKERRQ(ierr); 175 if (snes->reason) PetscFunctionReturn(0); 176 /* Select update between u^M and u^A */ 177 /* Decide whether to restart */ 178 } 179 } 180 snes->reason = SNES_DIVERGED_MAX_IT; 181 PetscFunctionReturn(0); 182 } 183 184 /*MC 185 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 186 187 Level: beginner 188 189 Notes: Supports only left preconditioning 190 191 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 192 SIAM Journal on Scientific Computing, 21(5), 2000. 193 194 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 195 M*/ 196 EXTERN_C_BEGIN 197 #undef __FUNCT__ 198 #define __FUNCT__ "SNESCreate_NGMRES" 199 PetscErrorCode SNESCreate_NGMRES(SNES snes) 200 { 201 SNES_NGMRES *ngmres; 202 PetscErrorCode ierr; 203 204 PetscFunctionBegin; 205 snes->ops->destroy = SNESDestroy_NGMRES; 206 snes->ops->setup = SNESSetUp_NGMRES; 207 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 208 snes->ops->view = SNESView_NGMRES; 209 snes->ops->solve = SNESSolve_NGMRES; 210 snes->ops->reset = SNESReset_NGMRES; 211 212 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 213 snes->data = (void*) ngmres; 214 ngmres->msize = 30; 215 ngmres->csize = 0; 216 217 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 218 #if 0 219 if (ksp->pc_side != PC_LEFT) {ierr = PetscInfo(ksp,"WARNING! Setting PC_SIDE for NGMRES to left!\n");CHKERRQ(ierr);} 220 snes->pc_side = PC_LEFT; 221 #endif 222 PetscFunctionReturn(0); 223 } 224 EXTERN_C_END 225