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