1 2 #include <../src/snes/impls/ls/lsimpl.h> 3 4 /* 5 Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function, 6 || F(u) ||_2 but not a zero, F(u) = 0. In the case when one cannot compute J^T F we use the fact that 7 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More 8 for this trick. One assumes that the probability that W is in the null space of J is very, very small. 9 */ 10 static PetscErrorCode SNESNEWTONLSCheckLocalMin_Private(SNES snes,Mat A,Vec F,PetscReal fnorm,PetscBool *ismin) 11 { 12 PetscReal a1; 13 PetscErrorCode ierr; 14 PetscBool hastranspose; 15 Vec W; 16 17 PetscFunctionBegin; 18 *ismin = PETSC_FALSE; 19 ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr); 20 ierr = VecDuplicate(F,&W);CHKERRQ(ierr); 21 if (hastranspose) { 22 /* Compute || J^T F|| */ 23 ierr = MatMultTranspose(A,F,W);CHKERRQ(ierr); 24 ierr = VecNorm(W,NORM_2,&a1);CHKERRQ(ierr); 25 ierr = PetscInfo1(snes,"|| J^T F|| %14.12e near zero implies found a local minimum\n",(double)(a1/fnorm));CHKERRQ(ierr); 26 if (a1/fnorm < 1.e-4) *ismin = PETSC_TRUE; 27 } else { 28 Vec work; 29 PetscScalar result; 30 PetscReal wnorm; 31 32 ierr = VecSetRandom(W,NULL);CHKERRQ(ierr); 33 ierr = VecNorm(W,NORM_2,&wnorm);CHKERRQ(ierr); 34 ierr = VecDuplicate(W,&work);CHKERRQ(ierr); 35 ierr = MatMult(A,W,work);CHKERRQ(ierr); 36 ierr = VecDot(F,work,&result);CHKERRQ(ierr); 37 ierr = VecDestroy(&work);CHKERRQ(ierr); 38 a1 = PetscAbsScalar(result)/(fnorm*wnorm); 39 ierr = PetscInfo1(snes,"(F^T J random)/(|| F ||*||random|| %14.12e near zero implies found a local minimum\n",(double)a1);CHKERRQ(ierr); 40 if (a1 < 1.e-4) *ismin = PETSC_TRUE; 41 } 42 ierr = VecDestroy(&W);CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 46 /* 47 Checks if J^T(F - J*X) = 0 48 */ 49 static PetscErrorCode SNESNEWTONLSCheckResidual_Private(SNES snes,Mat A,Vec F,Vec X) 50 { 51 PetscReal a1,a2; 52 PetscErrorCode ierr; 53 PetscBool hastranspose; 54 55 PetscFunctionBegin; 56 ierr = MatHasOperation(A,MATOP_MULT_TRANSPOSE,&hastranspose);CHKERRQ(ierr); 57 if (hastranspose) { 58 Vec W1,W2; 59 60 ierr = VecDuplicate(F,&W1);CHKERRQ(ierr); 61 ierr = VecDuplicate(F,&W2);CHKERRQ(ierr); 62 ierr = MatMult(A,X,W1);CHKERRQ(ierr); 63 ierr = VecAXPY(W1,-1.0,F);CHKERRQ(ierr); 64 65 /* Compute || J^T W|| */ 66 ierr = MatMultTranspose(A,W1,W2);CHKERRQ(ierr); 67 ierr = VecNorm(W1,NORM_2,&a1);CHKERRQ(ierr); 68 ierr = VecNorm(W2,NORM_2,&a2);CHKERRQ(ierr); 69 if (a1 != 0.0) { 70 ierr = PetscInfo1(snes,"||J^T(F-Ax)||/||F-AX|| %14.12e near zero implies inconsistent rhs\n",(double)(a2/a1));CHKERRQ(ierr); 71 } 72 ierr = VecDestroy(&W1);CHKERRQ(ierr); 73 ierr = VecDestroy(&W2);CHKERRQ(ierr); 74 } 75 PetscFunctionReturn(0); 76 } 77 78 /* -------------------------------------------------------------------- 79 80 This file implements a truncated Newton method with a line search, 81 for solving a system of nonlinear equations, using the KSP, Vec, 82 and Mat interfaces for linear solvers, vectors, and matrices, 83 respectively. 84 85 The following basic routines are required for each nonlinear solver: 86 SNESCreate_XXX() - Creates a nonlinear solver context 87 SNESSetFromOptions_XXX() - Sets runtime options 88 SNESSolve_XXX() - Solves the nonlinear system 89 SNESDestroy_XXX() - Destroys the nonlinear solver context 90 The suffix "_XXX" denotes a particular implementation, in this case 91 we use _NEWTONLS (e.g., SNESCreate_NEWTONLS, SNESSolve_NEWTONLS) for solving 92 systems of nonlinear equations with a line search (LS) method. 93 These routines are actually called via the common user interface 94 routines SNESCreate(), SNESSetFromOptions(), SNESSolve(), and 95 SNESDestroy(), so the application code interface remains identical 96 for all nonlinear solvers. 97 98 Another key routine is: 99 SNESSetUp_XXX() - Prepares for the use of a nonlinear solver 100 by setting data structures and options. The interface routine SNESSetUp() 101 is not usually called directly by the user, but instead is called by 102 SNESSolve() if necessary. 103 104 Additional basic routines are: 105 SNESView_XXX() - Prints details of runtime options that 106 have actually been used. 107 These are called by application codes via the interface routines 108 SNESView(). 109 110 The various types of solvers (preconditioners, Krylov subspace methods, 111 nonlinear solvers, timesteppers) are all organized similarly, so the 112 above description applies to these categories also. 113 114 -------------------------------------------------------------------- */ 115 /* 116 SNESSolve_NEWTONLS - Solves a nonlinear system with a truncated Newton 117 method with a line search. 118 119 Input Parameters: 120 . snes - the SNES context 121 122 Output Parameter: 123 . outits - number of iterations until termination 124 125 Application Interface Routine: SNESSolve() 126 127 Notes: 128 This implements essentially a truncated Newton method with a 129 line search. By default a cubic backtracking line search 130 is employed, as described in the text "Numerical Methods for 131 Unconstrained Optimization and Nonlinear Equations" by Dennis 132 and Schnabel. 133 */ 134 PetscErrorCode SNESSolve_NEWTONLS(SNES snes) 135 { 136 PetscErrorCode ierr; 137 PetscInt maxits,i,lits; 138 SNESLineSearchReason lssucceed; 139 PetscReal fnorm,gnorm,xnorm,ynorm; 140 Vec Y,X,F; 141 SNESLineSearch linesearch; 142 SNESConvergedReason reason; 143 144 PetscFunctionBegin; 145 if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name); 146 147 snes->numFailures = 0; 148 snes->numLinearSolveFailures = 0; 149 snes->reason = SNES_CONVERGED_ITERATING; 150 151 maxits = snes->max_its; /* maximum number of iterations */ 152 X = snes->vec_sol; /* solution vector */ 153 F = snes->vec_func; /* residual vector */ 154 Y = snes->vec_sol_update; /* newton step */ 155 156 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 157 snes->iter = 0; 158 snes->norm = 0.0; 159 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 160 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 161 162 /* compute the preconditioned function first in the case of left preconditioning with preconditioned function */ 163 if (snes->npc && snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 164 ierr = SNESApplyNPC(snes,X,NULL,F);CHKERRQ(ierr); 165 ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr); 166 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 167 snes->reason = SNES_DIVERGED_INNER; 168 PetscFunctionReturn(0); 169 } 170 171 ierr = VecNormBegin(F,NORM_2,&fnorm);CHKERRQ(ierr); 172 ierr = VecNormEnd(F,NORM_2,&fnorm);CHKERRQ(ierr); 173 } else { 174 if (!snes->vec_func_init_set) { 175 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 176 } else snes->vec_func_init_set = PETSC_FALSE; 177 } 178 179 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 180 SNESCheckFunctionNorm(snes,fnorm); 181 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 182 snes->norm = fnorm; 183 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 184 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 185 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 186 187 /* test convergence */ 188 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 189 if (snes->reason) PetscFunctionReturn(0); 190 191 for (i=0; i<maxits; i++) { 192 193 /* Call general purpose update function */ 194 if (snes->ops->update) { 195 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 196 } 197 198 /* apply the nonlinear preconditioner */ 199 if (snes->npc) { 200 if (snes->npcside== PC_RIGHT) { 201 ierr = SNESSetInitialFunction(snes->npc, F);CHKERRQ(ierr); 202 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->npc,X,snes->vec_rhs,0);CHKERRQ(ierr); 203 ierr = SNESSolve(snes->npc, snes->vec_rhs, X);CHKERRQ(ierr); 204 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->npc,X,snes->vec_rhs,0);CHKERRQ(ierr); 205 ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr); 206 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 207 snes->reason = SNES_DIVERGED_INNER; 208 PetscFunctionReturn(0); 209 } 210 ierr = SNESGetNPCFunction(snes,F,&fnorm);CHKERRQ(ierr); 211 } else if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 212 ierr = SNESApplyNPC(snes,X,F,F);CHKERRQ(ierr); 213 ierr = SNESGetConvergedReason(snes->npc,&reason);CHKERRQ(ierr); 214 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 215 snes->reason = SNES_DIVERGED_INNER; 216 PetscFunctionReturn(0); 217 } 218 } 219 } 220 221 /* Solve J Y = F, where J is Jacobian matrix */ 222 ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 223 SNESCheckJacobianDomainerror(snes); 224 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 225 ierr = KSPSolve(snes->ksp,F,Y);CHKERRQ(ierr); 226 SNESCheckKSPSolve(snes); 227 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 228 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 229 230 if (PetscLogPrintInfo) { 231 ierr = SNESNEWTONLSCheckResidual_Private(snes,snes->jacobian,F,Y);CHKERRQ(ierr); 232 } 233 234 /* Compute a (scaled) negative update in the line search routine: 235 X <- X - lambda*Y 236 and evaluate F = function(X) (depends on the line search). 237 */ 238 gnorm = fnorm; 239 ierr = SNESLineSearchApply(linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 240 ierr = SNESLineSearchGetReason(linesearch, &lssucceed);CHKERRQ(ierr); 241 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 242 ierr = PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)gnorm,(double)fnorm,(double)ynorm,(int)lssucceed);CHKERRQ(ierr); 243 if (snes->reason) break; 244 SNESCheckFunctionNorm(snes,fnorm); 245 if (lssucceed) { 246 if (snes->stol*xnorm > ynorm) { 247 snes->reason = SNES_CONVERGED_SNORM_RELATIVE; 248 PetscFunctionReturn(0); 249 } 250 if (++snes->numFailures >= snes->maxFailures) { 251 PetscBool ismin; 252 snes->reason = SNES_DIVERGED_LINE_SEARCH; 253 ierr = SNESNEWTONLSCheckLocalMin_Private(snes,snes->jacobian,F,fnorm,&ismin);CHKERRQ(ierr); 254 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 255 if (snes->errorifnotconverged && snes->reason) { 256 PetscViewer monitor; 257 ierr = SNESLineSearchGetDefaultMonitor(linesearch,&monitor);CHKERRQ(ierr); 258 if (!monitor) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged due to %s. Suggest running with -snes_linesearch_monitor",SNESConvergedReasons[snes->reason]); 259 else SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged due %s.",SNESConvergedReasons[snes->reason]); 260 } 261 break; 262 } 263 } 264 /* Monitor convergence */ 265 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 266 snes->iter = i+1; 267 snes->norm = fnorm; 268 snes->ynorm = ynorm; 269 snes->xnorm = xnorm; 270 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 271 ierr = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr); 272 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 273 /* Test for convergence */ 274 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 275 if (snes->reason) break; 276 } 277 if (i == maxits) { 278 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 279 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 280 } 281 PetscFunctionReturn(0); 282 } 283 /* -------------------------------------------------------------------------- */ 284 /* 285 SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use 286 of the SNESNEWTONLS nonlinear solver. 287 288 Input Parameter: 289 . snes - the SNES context 290 . x - the solution vector 291 292 Application Interface Routine: SNESSetUp() 293 294 Notes: 295 For basic use of the SNES solvers, the user need not explicitly call 296 SNESSetUp(), since these actions will automatically occur during 297 the call to SNESSolve(). 298 */ 299 PetscErrorCode SNESSetUp_NEWTONLS(SNES snes) 300 { 301 PetscErrorCode ierr; 302 303 PetscFunctionBegin; 304 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 305 if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED; 306 PetscFunctionReturn(0); 307 } 308 /* -------------------------------------------------------------------------- */ 309 310 PetscErrorCode SNESReset_NEWTONLS(SNES snes) 311 { 312 PetscFunctionBegin; 313 PetscFunctionReturn(0); 314 } 315 316 /* 317 SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created 318 with SNESCreate_NEWTONLS(). 319 320 Input Parameter: 321 . snes - the SNES context 322 323 Application Interface Routine: SNESDestroy() 324 */ 325 PetscErrorCode SNESDestroy_NEWTONLS(SNES snes) 326 { 327 PetscErrorCode ierr; 328 329 PetscFunctionBegin; 330 ierr = SNESReset_NEWTONLS(snes);CHKERRQ(ierr); 331 ierr = PetscFree(snes->data);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 /* -------------------------------------------------------------------------- */ 335 336 /* 337 SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure. 338 339 Input Parameters: 340 . SNES - the SNES context 341 . viewer - visualization context 342 343 Application Interface Routine: SNESView() 344 */ 345 static PetscErrorCode SNESView_NEWTONLS(SNES snes,PetscViewer viewer) 346 { 347 PetscErrorCode ierr; 348 PetscBool iascii; 349 350 PetscFunctionBegin; 351 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 352 if (iascii) { 353 } 354 PetscFunctionReturn(0); 355 } 356 357 /* -------------------------------------------------------------------------- */ 358 /* 359 SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method. 360 361 Input Parameter: 362 . snes - the SNES context 363 364 Application Interface Routine: SNESSetFromOptions() 365 */ 366 static PetscErrorCode SNESSetFromOptions_NEWTONLS(PetscOptionItems *PetscOptionsObject,SNES snes) 367 { 368 PetscFunctionBegin; 369 PetscFunctionReturn(0); 370 } 371 372 /* -------------------------------------------------------------------------- */ 373 /*MC 374 SNESNEWTONLS - Newton based nonlinear solver that uses a line search 375 376 Options Database: 377 + -snes_linesearch_type <bt> - bt,basic. Select line search type 378 . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt 379 . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch (SNESLineSearchSetComputeNorms()) 380 . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 381 . -snes_linesearch_maxstep <maxstep> - Sets the maximum stepsize the line search will use (if the 2-norm(y) > maxstep then scale y to be y = (maxstep/2-norm(y)) *y) 382 . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate 383 . -snes_linesearch_monitor - print information about progress of line searches 384 - -snes_linesearch_damping - damping factor used for basic line search 385 386 Notes: 387 This is the default nonlinear solver in SNES 388 389 Level: beginner 390 391 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONTR, SNESQN, SNESLineSearchSetType(), SNESLineSearchSetOrder() 392 SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck() SNESLineSearchSetComputeNorms() 393 394 M*/ 395 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes) 396 { 397 PetscErrorCode ierr; 398 SNES_NEWTONLS *neP; 399 SNESLineSearch linesearch; 400 401 PetscFunctionBegin; 402 snes->ops->setup = SNESSetUp_NEWTONLS; 403 snes->ops->solve = SNESSolve_NEWTONLS; 404 snes->ops->destroy = SNESDestroy_NEWTONLS; 405 snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS; 406 snes->ops->view = SNESView_NEWTONLS; 407 snes->ops->reset = SNESReset_NEWTONLS; 408 409 snes->npcside = PC_RIGHT; 410 snes->usesksp = PETSC_TRUE; 411 snes->usesnpc = PETSC_TRUE; 412 413 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 414 if (!((PetscObject)linesearch)->type_name) { 415 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);CHKERRQ(ierr); 416 } 417 418 snes->alwayscomputesfinalresidual = PETSC_TRUE; 419 420 ierr = PetscNewLog(snes,&neP);CHKERRQ(ierr); 421 snes->data = (void*)neP; 422 PetscFunctionReturn(0); 423 } 424