1 #include <../src/snes/impls/ncg/snesncgimpl.h> /*I "petscsnes.h" I*/ 2 const char *const SNESNCGTypes[] = {"FR", "PRP", "HS", "DY", "CD", "SNESNCGType", "SNES_NCG_", NULL}; 3 4 static PetscErrorCode SNESDestroy_NCG(SNES snes) 5 { 6 PetscFunctionBegin; 7 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNCGSetType_C", NULL)); 8 PetscCall(PetscFree(snes->data)); 9 PetscFunctionReturn(PETSC_SUCCESS); 10 } 11 12 static PetscErrorCode SNESSetUp_NCG(SNES snes) 13 { 14 PetscFunctionBegin; 15 PetscCall(SNESSetWorkVecs(snes, 2)); 16 PetscCheck(snes->npcside != PC_RIGHT, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNESNCG only supports left preconditioning"); 17 if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED; 18 PetscFunctionReturn(PETSC_SUCCESS); 19 } 20 21 static PetscErrorCode SNESLineSearchApply_NCGLinear(SNESLineSearch linesearch) 22 { 23 PetscScalar alpha, ptAp; 24 Vec X, Y, F, W; 25 SNES snes; 26 PetscReal *fnorm, *xnorm, *ynorm; 27 28 PetscFunctionBegin; 29 PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 30 X = linesearch->vec_sol; 31 W = linesearch->vec_sol_new; 32 F = linesearch->vec_func; 33 Y = linesearch->vec_update; 34 fnorm = &linesearch->fnorm; 35 xnorm = &linesearch->xnorm; 36 ynorm = &linesearch->ynorm; 37 38 if (!snes->jacobian) PetscCall(SNESSetUpMatrices(snes)); 39 40 /* 41 The exact step size for unpreconditioned linear CG is just: 42 alpha = (r, r) / (p, Ap) = (f, f) / (y, Jy) 43 */ 44 PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre)); 45 PetscCall(VecDot(F, F, &alpha)); 46 PetscCall(MatMult(snes->jacobian, Y, W)); 47 PetscCall(VecDot(Y, W, &ptAp)); 48 alpha = alpha / ptAp; 49 PetscCall(VecAXPY(X, -alpha, Y)); 50 PetscCall(SNESComputeFunction(snes, X, F)); 51 52 PetscCall(VecNorm(F, NORM_2, fnorm)); 53 PetscCall(VecNorm(X, NORM_2, xnorm)); 54 PetscCall(VecNorm(Y, NORM_2, ynorm)); 55 PetscFunctionReturn(PETSC_SUCCESS); 56 } 57 58 /*MC 59 SNESLINESEARCHNCGLINEAR - Special line search only for the nonlinear CG solver `SNESNCG` 60 61 This line search uses the length "as if" the problem is linear (that is what is computed by the linear CG method) using the Jacobian of the function. 62 alpha = (r, r) / (p, Ap) = (f, f) / (y, Jy) where r (f) is the current residual (function value), p (y) is the current search direction. 63 64 Level: advanced 65 66 Notes: 67 This requires a Jacobian-vector product but does not require the solution of a linear system with the Jacobian 68 69 This is a "odd-ball" line search, we don't know if it is in the literature or used in practice by anyone. 70 71 .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()` 72 M*/ 73 74 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NCGLinear(SNESLineSearch linesearch) 75 { 76 PetscFunctionBegin; 77 linesearch->ops->apply = SNESLineSearchApply_NCGLinear; 78 linesearch->ops->destroy = NULL; 79 linesearch->ops->setfromoptions = NULL; 80 linesearch->ops->reset = NULL; 81 linesearch->ops->view = NULL; 82 linesearch->ops->setup = NULL; 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 86 static PetscErrorCode SNESSetFromOptions_NCG(SNES snes, PetscOptionItems PetscOptionsObject) 87 { 88 SNES_NCG *ncg = (SNES_NCG *)snes->data; 89 PetscBool debug = PETSC_FALSE; 90 SNESNCGType ncgtype = ncg->type; 91 SNESLineSearch linesearch; 92 93 PetscFunctionBegin; 94 PetscOptionsHeadBegin(PetscOptionsObject, "SNES NCG options"); 95 PetscCall(PetscOptionsBool("-snes_ncg_monitor", "Monitor the beta values used in the NCG iterations", "SNES", ncg->monitor ? PETSC_TRUE : PETSC_FALSE, &debug, NULL)); 96 if (debug) ncg->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)); 97 PetscCall(PetscOptionsEnum("-snes_ncg_type", "NCG Beta type used", "SNESNCGSetType", SNESNCGTypes, (PetscEnum)ncg->type, (PetscEnum *)&ncgtype, NULL)); 98 PetscCall(SNESNCGSetType(snes, ncgtype)); 99 PetscOptionsHeadEnd(); 100 if (!snes->linesearch) { 101 PetscCall(SNESGetLineSearch(snes, &linesearch)); 102 if (!((PetscObject)linesearch)->type_name) { 103 if (!snes->npc) { 104 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHCP)); 105 } else { 106 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHL2)); 107 } 108 } 109 } 110 PetscFunctionReturn(PETSC_SUCCESS); 111 } 112 113 static PetscErrorCode SNESView_NCG(SNES snes, PetscViewer viewer) 114 { 115 SNES_NCG *ncg = (SNES_NCG *)snes->data; 116 PetscBool iascii; 117 118 PetscFunctionBegin; 119 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 120 if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " type: %s\n", SNESNCGTypes[ncg->type])); 121 PetscFunctionReturn(PETSC_SUCCESS); 122 } 123 124 /*@ 125 SNESNCGSetType - Sets the conjugate update type for nonlinear CG `SNESNCG`. 126 127 Logically Collective 128 129 Input Parameters: 130 + snes - the iterative context 131 - btype - update type, see `SNESNCGType` 132 133 Options Database Key: 134 . -snes_ncg_type <prp,fr,hs,dy,cd> - strategy for selecting algorithm for computing beta 135 136 Level: intermediate 137 138 Notes: 139 `SNES_NCG_PRP` is the default, and the only one that tolerates generalized search directions. 140 141 It is not clear what "generalized search directions" means, does it mean use with a nonlinear preconditioner, 142 that is using -npc_snes_type <type>, `SNESSetNPC()`, or `SNESGetNPC()`? 143 144 .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESNCGType`, `SNES_NCG_FR`, `SNES_NCG_PRP`, `SNES_NCG_HS`, `SNES_NCG_DY`, `SNES_NCG_CD` 145 @*/ 146 PetscErrorCode SNESNCGSetType(SNES snes, SNESNCGType btype) 147 { 148 PetscFunctionBegin; 149 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 150 PetscTryMethod(snes, "SNESNCGSetType_C", (SNES, SNESNCGType), (snes, btype)); 151 PetscFunctionReturn(PETSC_SUCCESS); 152 } 153 154 static PetscErrorCode SNESNCGSetType_NCG(SNES snes, SNESNCGType btype) 155 { 156 SNES_NCG *ncg = (SNES_NCG *)snes->data; 157 158 PetscFunctionBegin; 159 ncg->type = btype; 160 PetscFunctionReturn(PETSC_SUCCESS); 161 } 162 163 /* 164 SNESSolve_NCG - Solves a nonlinear system with the Nonlinear Conjugate Gradient method. 165 166 Input Parameter: 167 . snes - the `SNES` context 168 169 Output Parameter: 170 . outits - number of iterations until termination 171 172 Application Interface Routine: SNESSolve() 173 */ 174 static PetscErrorCode SNESSolve_NCG(SNES snes) 175 { 176 SNES_NCG *ncg = (SNES_NCG *)snes->data; 177 Vec X, dX, lX, F, dXold; 178 PetscReal fnorm, ynorm, xnorm, beta = 0.0; 179 PetscScalar dXdotdX, dXolddotdXold, dXdotdXold, lXdotdX, lXdotdXold; 180 PetscInt maxits, i; 181 SNESLineSearchReason lsresult = SNES_LINESEARCH_SUCCEEDED; 182 SNESLineSearch linesearch; 183 SNESConvergedReason reason; 184 185 PetscFunctionBegin; 186 PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name); 187 188 PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite)); 189 snes->reason = SNES_CONVERGED_ITERATING; 190 191 maxits = snes->max_its; /* maximum number of iterations */ 192 X = snes->vec_sol; /* X^n */ 193 dXold = snes->work[0]; /* The previous iterate of X */ 194 dX = snes->work[1]; /* the preconditioned direction */ 195 lX = snes->vec_sol_update; /* the conjugate direction */ 196 F = snes->vec_func; /* residual vector */ 197 198 PetscCall(SNESGetLineSearch(snes, &linesearch)); 199 200 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 201 snes->iter = 0; 202 snes->norm = 0.; 203 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 204 205 /* compute the initial function and preconditioned update dX */ 206 207 if (snes->npc && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 208 PetscCall(SNESApplyNPC(snes, X, NULL, dX)); 209 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 210 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 211 snes->reason = SNES_DIVERGED_INNER; 212 PetscFunctionReturn(PETSC_SUCCESS); 213 } 214 PetscCall(VecCopy(dX, F)); 215 PetscCall(VecNorm(F, NORM_2, &fnorm)); 216 } else { 217 if (!snes->vec_func_init_set) { 218 PetscCall(SNESComputeFunction(snes, X, F)); 219 } else snes->vec_func_init_set = PETSC_FALSE; 220 221 /* convergence test */ 222 PetscCall(VecNorm(F, NORM_2, &fnorm)); 223 SNESCheckFunctionNorm(snes, fnorm); 224 PetscCall(VecCopy(F, dX)); 225 } 226 if (snes->npc) { 227 if (snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 228 PetscCall(SNESApplyNPC(snes, X, F, dX)); 229 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 230 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 231 snes->reason = SNES_DIVERGED_INNER; 232 PetscFunctionReturn(PETSC_SUCCESS); 233 } 234 } 235 } 236 PetscCall(VecCopy(dX, lX)); 237 PetscCall(VecDot(dX, dX, &dXdotdX)); 238 239 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 240 snes->norm = fnorm; 241 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 242 PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0)); 243 244 /* test convergence */ 245 PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm)); 246 PetscCall(SNESMonitor(snes, 0, fnorm)); 247 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS); 248 249 /* Call general purpose update function */ 250 PetscTryTypeMethod(snes, update, snes->iter); 251 252 /* first update -- just use the (preconditioned) residual direction for the initial conjugate direction */ 253 254 for (i = 1; i < maxits + 1; i++) { 255 /* some update types require the old update direction or conjugate direction */ 256 if (ncg->type != SNES_NCG_FR) PetscCall(VecCopy(dX, dXold)); 257 PetscCall(SNESLineSearchApply(linesearch, X, F, &fnorm, lX)); 258 PetscCall(SNESLineSearchGetReason(linesearch, &lsresult)); 259 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm)); 260 if (lsresult) { 261 if (++snes->numFailures >= snes->maxFailures) { 262 snes->reason = SNES_DIVERGED_LINE_SEARCH; 263 PetscFunctionReturn(PETSC_SUCCESS); 264 } 265 } 266 if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 267 snes->reason = SNES_DIVERGED_FUNCTION_COUNT; 268 PetscFunctionReturn(PETSC_SUCCESS); 269 } 270 /* Monitor convergence */ 271 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes)); 272 snes->iter = i; 273 snes->norm = fnorm; 274 snes->xnorm = xnorm; 275 snes->ynorm = ynorm; 276 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes)); 277 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0)); 278 279 /* Test for convergence */ 280 PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm)); 281 PetscCall(SNESMonitor(snes, snes->iter, snes->norm)); 282 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS); 283 284 /* Call general purpose update function */ 285 PetscTryTypeMethod(snes, update, snes->iter); 286 if (snes->npc) { 287 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) { 288 PetscCall(SNESApplyNPC(snes, X, NULL, dX)); 289 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 290 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 291 snes->reason = SNES_DIVERGED_INNER; 292 PetscFunctionReturn(PETSC_SUCCESS); 293 } 294 PetscCall(VecCopy(dX, F)); 295 } else { 296 PetscCall(SNESApplyNPC(snes, X, F, dX)); 297 PetscCall(SNESGetConvergedReason(snes->npc, &reason)); 298 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 299 snes->reason = SNES_DIVERGED_INNER; 300 PetscFunctionReturn(PETSC_SUCCESS); 301 } 302 } 303 } else { 304 PetscCall(VecCopy(F, dX)); 305 } 306 307 /* compute the conjugate direction lX = dX + beta*lX with beta = ((dX, dX) / (dX_old, dX_old) (Fletcher-Reeves update)*/ 308 switch (ncg->type) { 309 case SNES_NCG_FR: /* Fletcher-Reeves */ 310 dXolddotdXold = dXdotdX; 311 PetscCall(VecDot(dX, dX, &dXdotdX)); 312 beta = PetscRealPart(dXdotdX / dXolddotdXold); 313 break; 314 case SNES_NCG_PRP: /* Polak-Ribiere-Poylak */ 315 dXolddotdXold = dXdotdX; 316 PetscCall(VecDotBegin(dX, dX, &dXdotdX)); 317 PetscCall(VecDotBegin(dXold, dX, &dXdotdXold)); 318 PetscCall(VecDotEnd(dX, dX, &dXdotdX)); 319 PetscCall(VecDotEnd(dXold, dX, &dXdotdXold)); 320 beta = PetscRealPart((dXdotdX - dXdotdXold) / dXolddotdXold); 321 if (beta < 0.0) beta = 0.0; /* restart */ 322 break; 323 case SNES_NCG_HS: /* Hestenes-Stiefel */ 324 PetscCall(VecDotBegin(dX, dX, &dXdotdX)); 325 PetscCall(VecDotBegin(dX, dXold, &dXdotdXold)); 326 PetscCall(VecDotBegin(lX, dX, &lXdotdX)); 327 PetscCall(VecDotBegin(lX, dXold, &lXdotdXold)); 328 PetscCall(VecDotEnd(dX, dX, &dXdotdX)); 329 PetscCall(VecDotEnd(dX, dXold, &dXdotdXold)); 330 PetscCall(VecDotEnd(lX, dX, &lXdotdX)); 331 PetscCall(VecDotEnd(lX, dXold, &lXdotdXold)); 332 beta = PetscRealPart((dXdotdX - dXdotdXold) / (lXdotdX - lXdotdXold)); 333 break; 334 case SNES_NCG_DY: /* Dai-Yuan */ 335 PetscCall(VecDotBegin(dX, dX, &dXdotdX)); 336 PetscCall(VecDotBegin(lX, dX, &lXdotdX)); 337 PetscCall(VecDotBegin(lX, dXold, &lXdotdXold)); 338 PetscCall(VecDotEnd(dX, dX, &dXdotdX)); 339 PetscCall(VecDotEnd(lX, dX, &lXdotdX)); 340 PetscCall(VecDotEnd(lX, dXold, &lXdotdXold)); 341 beta = PetscRealPart(dXdotdX / (lXdotdXold - lXdotdX)); 342 break; 343 case SNES_NCG_CD: /* Conjugate Descent */ 344 PetscCall(VecDotBegin(dX, dX, &dXdotdX)); 345 PetscCall(VecDotBegin(lX, dXold, &lXdotdXold)); 346 PetscCall(VecDotEnd(dX, dX, &dXdotdX)); 347 PetscCall(VecDotEnd(lX, dXold, &lXdotdXold)); 348 beta = PetscRealPart(dXdotdX / lXdotdXold); 349 break; 350 } 351 if (ncg->monitor) PetscCall(PetscViewerASCIIPrintf(ncg->monitor, "beta = %e\n", (double)beta)); 352 PetscCall(VecAYPX(lX, beta, dX)); 353 } 354 PetscFunctionReturn(PETSC_SUCCESS); 355 } 356 357 /*MC 358 SNESNCG - Nonlinear Conjugate-Gradient method for the solution of nonlinear systems {cite}`bruneknepleysmithtu15`. 359 360 Level: beginner 361 362 Options Database Keys: 363 + -snes_ncg_type <fr, prp, dy, hs, cd> - Choice of conjugate-gradient update parameter, default is `prp`. 364 . -snes_linesearch_type <cp,l2,basic> - Line search type. 365 - -snes_ncg_monitor - Print the beta values nonlinear Conjugate-Gradient used in the iteration, . 366 367 Notes: 368 This solves the nonlinear system of equations $ F(x) = 0 $ using the nonlinear generalization of the conjugate 369 gradient method. This may be used with a nonlinear preconditioner used to pick the new search directions, but otherwise 370 chooses the initial search direction as $ F(x) $ for the initial guess $x$. 371 372 Only supports left non-linear preconditioning. 373 374 Default line search is `SNESLINESEARCHCP`, unless a nonlinear preconditioner is used with `-npc_snes_type` <type>, `SNESSetNPC()`, or `SNESGetNPC()` then 375 `SNESLINESEARCHL2` is used. Also supports the special-purpose line search `SNESLINESEARCHNCGLINEAR` 376 377 .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNGMRES`, `SNESQN`, `SNESLINESEARCHNCGLINEAR`, `SNESNCGSetType()`, `SNESLineSearchSetType()` 378 M*/ 379 PETSC_EXTERN PetscErrorCode SNESCreate_NCG(SNES snes) 380 { 381 SNES_NCG *neP; 382 383 PetscFunctionBegin; 384 snes->ops->destroy = SNESDestroy_NCG; 385 snes->ops->setup = SNESSetUp_NCG; 386 snes->ops->setfromoptions = SNESSetFromOptions_NCG; 387 snes->ops->view = SNESView_NCG; 388 snes->ops->solve = SNESSolve_NCG; 389 390 snes->usesksp = PETSC_FALSE; 391 snes->usesnpc = PETSC_TRUE; 392 snes->npcside = PC_LEFT; 393 394 snes->alwayscomputesfinalresidual = PETSC_TRUE; 395 396 PetscCall(SNESParametersInitialize(snes)); 397 PetscObjectParameterSetDefault(snes, max_funcs, 30000); 398 PetscObjectParameterSetDefault(snes, max_its, 10000); 399 PetscObjectParameterSetDefault(snes, stol, 1e-20); 400 401 PetscCall(PetscNew(&neP)); 402 snes->data = (void *)neP; 403 neP->monitor = NULL; 404 neP->type = SNES_NCG_PRP; 405 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNCGSetType_C", SNESNCGSetType_NCG)); 406 PetscFunctionReturn(PETSC_SUCCESS); 407 } 408