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, Y, 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 Y = snes->vec_sol_update; 100 F = snes->vec_func; 101 rOld = snes->work[0]; 102 103 ierr = SNESGetPC(snes, &pc);CHKERRQ(ierr); 104 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 105 snes->iter = 0; 106 snes->norm = 0.; 107 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 108 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr); 109 if (snes->domainerror) { 110 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 111 PetscFunctionReturn(0); 112 } 113 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 114 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FP, "Infinite or not-a-number generated in norm"); 115 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 116 snes->norm = fnorm; 117 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 118 SNESLogConvHistory(snes, fnorm, 0); 119 ierr = SNESMonitor(snes, 0, fnorm);CHKERRQ(ierr); 120 121 /* set parameter for default relative tolerance convergence test */ 122 snes->ttol = fnorm*snes->rtol; 123 /* test convergence */ 124 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 125 if (snes->reason) PetscFunctionReturn(0); 126 127 #ifndef OLD_BARRY_CODE 128 ierr = VecCopy(F, rOld);CHKERRQ(ierr); 129 #else 130 /* Barry: What the heck is this part doing? */ 131 /* determine optimal scale factor -- slow code */ 132 ierr = VecDuplicate(P, &y);CHKERRQ(ierr); 133 ierr = VecDuplicate(P, &w);CHKERRQ(ierr); 134 ierr = MatMult(Amat, Pold, y);CHKERRQ(ierr); 135 /*ierr = KSP_PCApplyBAorAB(ksp,Pold,y,w);CHKERRQ(ierr); */ /* y = BAp */ 136 ierr = VecDotNorm2(Pold, y, &rdot, &abr);CHKERRQ(ierr); /* rdot = (p)^T(BAp); abr = (BAp)^T (BAp) */ 137 ierr = VecDestroy(&y);CHKERRQ(ierr); 138 ierr = VecDestroy(&w);CHKERRQ(ierr); 139 A0 = rdot/abr; 140 ierr = VecAXPY(X,A0,Pold);CHKERRQ(ierr); /* x <- x + scale p */ 141 #endif 142 143 /* Loop over batches of directions */ 144 /* Code from Barry I do not understand to solve the least-squares problem, Time to try again */ 145 for(k = 0; k < snes->max_its; k += ngmres->msize) { 146 /* Loop over updates for this batch */ 147 /* TODO: Incorporate the variant which use the analytic Jacobian */ 148 /* TODO: Incorporate selection criteria for u^A from paper (need to store residual and update norms) */ 149 /* TODO: Incorporate criteria for restarting from paper */ 150 for(i = 0; i < ngmres->msize && k+i < snes->max_its; ++i) { 151 /* Get a new approxiate solution u^k */ 152 ierr = SNESSolve(pc, PETSC_NULL, X);CHKERRQ(ierr); 153 /* Solve least squares problem using last msize iterates */ 154 ierr = SNESGetFunction(pc, &r, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); /* r = F(u^M) */ 155 for(j = 0; j < i; ++j) { /* r = \prod_j (I + v_j w^T_j) r */ 156 ierr = VecDot(W[j], r, &wdot);CHKERRQ(ierr); 157 ierr = VecAXPY(r, wdot, V[j]);CHKERRQ(ierr); 158 } 159 ierr = VecCopy(rOld, W[i]);CHKERRQ(ierr); /* w_i = r_{old} */ 160 ierr = VecAXPY(rOld, -1.0, r);CHKERRQ(ierr); /* v_i = r */ 161 ierr = VecDot(W[i], rOld, &wdot);CHKERRQ(ierr); /* ------------------- */ 162 ierr = VecCopy(r, V[i]);CHKERRQ(ierr); /* w^T_i (r_{old} - r) */ 163 ierr = VecScale(V[i], 1.0/wdot);CHKERRQ(ierr); 164 ierr = VecDot(W[i], r, &wdot);CHKERRQ(ierr); /* r = (I + v_i w^T_i) r */ 165 ierr = VecAXPY(r, wdot, V[i]);CHKERRQ(ierr); 166 ierr = VecCopy(r, rOld);CHKERRQ(ierr); /* r_{old} = r */ 167 168 ierr = VecAXPY(X, 1.0, r);CHKERRQ(ierr); /* u^A = u^M + r */ 169 /* Monitor and log history */ 170 ierr = SNESComputeFunction(snes, X, F);CHKERRQ(ierr); /* F(u^k) */ 171 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); 172 SNESLogConvHistory(snes, fnorm, 0); 173 ierr = SNESMonitor(snes, i+k+1, fnorm);CHKERRQ(ierr); 174 /* Check convergence */ 175 ierr = (*snes->ops->converged)(snes, i+k+1, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);CHKERRQ(ierr); 176 if (snes->reason) PetscFunctionReturn(0); 177 /* Select update between u^M and u^A */ 178 /* Decide whether to restart */ 179 } 180 } 181 snes->reason = SNES_DIVERGED_MAX_IT; 182 PetscFunctionReturn(0); 183 } 184 185 /*MC 186 SNESNGMRES - The Nonlinear Generalized Minimum Residual (NGMRES) method of Oosterlee and Washio. 187 188 Level: beginner 189 190 Notes: Supports only left preconditioning 191 192 "Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows", C. W. Oosterlee and T. Washio, 193 SIAM Journal on Scientific Computing, 21(5), 2000. 194 195 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 196 M*/ 197 EXTERN_C_BEGIN 198 #undef __FUNCT__ 199 #define __FUNCT__ "SNESCreate_NGMRES" 200 PetscErrorCode SNESCreate_NGMRES(SNES snes) 201 { 202 SNES_NGMRES *ngmres; 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 snes->ops->destroy = SNESDestroy_NGMRES; 207 snes->ops->setup = SNESSetUp_NGMRES; 208 snes->ops->setfromoptions = SNESSetFromOptions_NGMRES; 209 snes->ops->view = SNESView_NGMRES; 210 snes->ops->solve = SNESSolve_NGMRES; 211 snes->ops->reset = SNESReset_NGMRES; 212 213 ierr = PetscNewLog(snes, SNES_NGMRES, &ngmres);CHKERRQ(ierr); 214 snes->data = (void*) ngmres; 215 ngmres->msize = 30; 216 ngmres->csize = 0; 217 218 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 219 #if 0 220 if (ksp->pc_side != PC_LEFT) {ierr = PetscInfo(ksp,"WARNING! Setting PC_SIDE for NGMRES to left!\n");CHKERRQ(ierr);} 221 snes->pc_side = PC_LEFT; 222 #endif 223 PetscFunctionReturn(0); 224 } 225 EXTERN_C_END 226