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