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