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