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