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 snes->numFailures = 0; 154 snes->numLinearSolveFailures = 0; 155 snes->reason = SNES_CONVERGED_ITERATING; 156 157 maxits = snes->max_its; /* maximum number of iterations */ 158 X = snes->vec_sol; /* solution vector */ 159 F = snes->vec_func; /* residual vector */ 160 Y = snes->vec_sol_update; /* newton step */ 161 162 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 163 snes->iter = 0; 164 snes->norm = 0.0; 165 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 166 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 167 168 /* compute the preconditioned function first in the case of left preconditioning with preconditioned function */ 169 if (snes->pc && snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 170 ierr = SNESApplyNPC(snes,X,NULL,F);CHKERRQ(ierr); 171 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 172 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 173 snes->reason = SNES_DIVERGED_INNER; 174 PetscFunctionReturn(0); 175 } 176 177 ierr = VecNormBegin(F,NORM_2,&fnorm);CHKERRQ(ierr); 178 ierr = VecNormEnd(F,NORM_2,&fnorm);CHKERRQ(ierr); 179 } else { 180 if (!snes->vec_func_init_set) { 181 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 182 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 183 if (domainerror) { 184 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 185 PetscFunctionReturn(0); 186 } 187 } else snes->vec_func_init_set = PETSC_FALSE; 188 } 189 190 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 191 if (PetscIsInfOrNanReal(fnorm)) { 192 snes->reason = SNES_DIVERGED_FNORM_NAN; 193 PetscFunctionReturn(0); 194 } 195 196 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 197 snes->norm = fnorm; 198 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 199 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 200 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 201 202 /* test convergence */ 203 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 204 if (snes->reason) PetscFunctionReturn(0); 205 206 for (i=0; i<maxits; i++) { 207 208 /* Call general purpose update function */ 209 if (snes->ops->update) { 210 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 211 } 212 213 /* apply the nonlinear preconditioner */ 214 if (snes->pc) { 215 if (snes->pcside == PC_RIGHT) { 216 ierr = SNESSetInitialFunction(snes->pc, F);CHKERRQ(ierr); 217 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,snes->vec_rhs,0);CHKERRQ(ierr); 218 ierr = SNESSolve(snes->pc, snes->vec_rhs, X);CHKERRQ(ierr); 219 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,snes->vec_rhs,0);CHKERRQ(ierr); 220 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 221 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 222 snes->reason = SNES_DIVERGED_INNER; 223 PetscFunctionReturn(0); 224 } 225 ierr = SNESGetNPCFunction(snes,F,&fnorm);CHKERRQ(ierr); 226 } else if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 227 ierr = SNESApplyNPC(snes,X,F,F);CHKERRQ(ierr); 228 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 229 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 230 snes->reason = SNES_DIVERGED_INNER; 231 PetscFunctionReturn(0); 232 } 233 } 234 } 235 236 /* Solve J Y = F, where J is Jacobian matrix */ 237 ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 238 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 239 ierr = KSPSolve(snes->ksp,F,Y);CHKERRQ(ierr); 240 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 241 if (kspreason < 0) { 242 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 243 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 244 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 245 break; 246 } 247 } 248 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 249 snes->linear_its += lits; 250 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 251 252 if (PetscLogPrintInfo) { 253 ierr = SNESNEWTONLSCheckResidual_Private(snes,snes->jacobian,F,Y);CHKERRQ(ierr); 254 } 255 256 /* Compute a (scaled) negative update in the line search routine: 257 X <- X - lambda*Y 258 and evaluate F = function(X) (depends on the line search). 259 */ 260 gnorm = fnorm; 261 ierr = SNESLineSearchApply(linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 262 ierr = SNESLineSearchGetSuccess(linesearch, &lssucceed);CHKERRQ(ierr); 263 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 264 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); 265 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 266 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 267 if (domainerror) { 268 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 269 PetscFunctionReturn(0); 270 } 271 if (!lssucceed) { 272 if (snes->stol*xnorm > ynorm) { 273 snes->reason = SNES_CONVERGED_SNORM_RELATIVE; 274 PetscFunctionReturn(0); 275 } 276 if (++snes->numFailures >= snes->maxFailures) { 277 PetscBool ismin; 278 snes->reason = SNES_DIVERGED_LINE_SEARCH; 279 ierr = SNESNEWTONLSCheckLocalMin_Private(snes,snes->jacobian,F,fnorm,&ismin);CHKERRQ(ierr); 280 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 281 break; 282 } 283 } 284 /* Monitor convergence */ 285 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 286 snes->iter = i+1; 287 snes->norm = fnorm; 288 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 289 ierr = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr); 290 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 291 /* Test for convergence */ 292 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 293 if (snes->reason) break; 294 } 295 if (i == maxits) { 296 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 297 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 298 } 299 PetscFunctionReturn(0); 300 } 301 /* -------------------------------------------------------------------------- */ 302 /* 303 SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use 304 of the SNESNEWTONLS nonlinear solver. 305 306 Input Parameter: 307 . snes - the SNES context 308 . x - the solution vector 309 310 Application Interface Routine: SNESSetUp() 311 312 Notes: 313 For basic use of the SNES solvers, the user need not explicitly call 314 SNESSetUp(), since these actions will automatically occur during 315 the call to SNESSolve(). 316 */ 317 #undef __FUNCT__ 318 #define __FUNCT__ "SNESSetUp_NEWTONLS" 319 PetscErrorCode SNESSetUp_NEWTONLS(SNES snes) 320 { 321 PetscErrorCode ierr; 322 323 PetscFunctionBegin; 324 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 325 if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED; 326 PetscFunctionReturn(0); 327 } 328 /* -------------------------------------------------------------------------- */ 329 330 #undef __FUNCT__ 331 #define __FUNCT__ "SNESReset_NEWTONLS" 332 PetscErrorCode SNESReset_NEWTONLS(SNES snes) 333 { 334 PetscFunctionBegin; 335 PetscFunctionReturn(0); 336 } 337 338 /* 339 SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created 340 with SNESCreate_NEWTONLS(). 341 342 Input Parameter: 343 . snes - the SNES context 344 345 Application Interface Routine: SNESDestroy() 346 */ 347 #undef __FUNCT__ 348 #define __FUNCT__ "SNESDestroy_NEWTONLS" 349 PetscErrorCode SNESDestroy_NEWTONLS(SNES snes) 350 { 351 PetscErrorCode ierr; 352 353 PetscFunctionBegin; 354 ierr = SNESReset_NEWTONLS(snes);CHKERRQ(ierr); 355 ierr = PetscFree(snes->data);CHKERRQ(ierr); 356 PetscFunctionReturn(0); 357 } 358 /* -------------------------------------------------------------------------- */ 359 360 /* 361 SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure. 362 363 Input Parameters: 364 . SNES - the SNES context 365 . viewer - visualization context 366 367 Application Interface Routine: SNESView() 368 */ 369 #undef __FUNCT__ 370 #define __FUNCT__ "SNESView_NEWTONLS" 371 static PetscErrorCode SNESView_NEWTONLS(SNES snes,PetscViewer viewer) 372 { 373 PetscErrorCode ierr; 374 PetscBool iascii; 375 376 PetscFunctionBegin; 377 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 378 if (iascii) { 379 } 380 PetscFunctionReturn(0); 381 } 382 383 /* -------------------------------------------------------------------------- */ 384 /* 385 SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method. 386 387 Input Parameter: 388 . snes - the SNES context 389 390 Application Interface Routine: SNESSetFromOptions() 391 */ 392 #undef __FUNCT__ 393 #define __FUNCT__ "SNESSetFromOptions_NEWTONLS" 394 static PetscErrorCode SNESSetFromOptions_NEWTONLS(PetscOptions *PetscOptionsObject,SNES snes) 395 { 396 PetscErrorCode ierr; 397 SNESLineSearch linesearch; 398 399 PetscFunctionBegin; 400 if (!snes->linesearch) { 401 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 402 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);CHKERRQ(ierr); 403 } 404 PetscFunctionReturn(0); 405 } 406 407 /* -------------------------------------------------------------------------- */ 408 /*MC 409 SNESNEWTONLS - Newton based nonlinear solver that uses a line search 410 411 Options Database: 412 + -snes_linesearch_type <bt> - bt,basic. Select line search type 413 . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt 414 . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch 415 . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 416 . -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) 417 . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate 418 . -snes_linesearch_monitor - print information about progress of line searches 419 - -snes_linesearch_damping - damping factor used for basic line search 420 421 Notes: This is the default nonlinear solver in SNES 422 423 Level: beginner 424 425 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONTR, SNESQN, SNESLineSearchSetType(), SNESLineSearchSetOrder() 426 SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck() SNESLineSearchSetComputeNorms() 427 428 M*/ 429 #undef __FUNCT__ 430 #define __FUNCT__ "SNESCreate_NEWTONLS" 431 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes) 432 { 433 PetscErrorCode ierr; 434 SNES_NEWTONLS *neP; 435 436 PetscFunctionBegin; 437 snes->ops->setup = SNESSetUp_NEWTONLS; 438 snes->ops->solve = SNESSolve_NEWTONLS; 439 snes->ops->destroy = SNESDestroy_NEWTONLS; 440 snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS; 441 snes->ops->view = SNESView_NEWTONLS; 442 snes->ops->reset = SNESReset_NEWTONLS; 443 444 snes->pcside = PC_RIGHT; 445 snes->usesksp = PETSC_TRUE; 446 snes->usespc = PETSC_TRUE; 447 ierr = PetscNewLog(snes,&neP);CHKERRQ(ierr); 448 snes->data = (void*)neP; 449 PetscFunctionReturn(0); 450 } 451