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