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 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 137 PetscReal fnorm,gnorm,xnorm,ynorm; 138 Vec Y,X,F,G,W; 139 KSPConvergedReason kspreason; 140 PetscBool domainerror; 141 SNESLineSearch linesearch; 142 SNESConvergedReason reason; 143 PCSide npcside; 144 145 PetscFunctionBegin; 146 snes->numFailures = 0; 147 snes->numLinearSolveFailures = 0; 148 snes->reason = SNES_CONVERGED_ITERATING; 149 150 maxits = snes->max_its; /* maximum number of iterations */ 151 X = snes->vec_sol; /* solution vector */ 152 F = snes->vec_func; /* residual vector */ 153 Y = snes->vec_sol_update; /* newton step */ 154 G = snes->work[0]; 155 W = snes->work[1]; 156 157 ierr = PetscObjectAMSTakeAccess((PetscObject)snes);CHKERRQ(ierr); 158 snes->iter = 0; 159 snes->norm = 0.0; 160 ierr = PetscObjectAMSGrantAccess((PetscObject)snes);CHKERRQ(ierr); 161 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 162 163 /* compute the preconditioned function first in the case of left preconditioning with preconditioned function */ 164 if (snes->pc && snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) { 165 ierr = SNESApplyPC(snes,X,PETSC_NULL,PETSC_NULL,F);CHKERRQ(ierr); 166 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 167 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 168 snes->reason = SNES_DIVERGED_INNER; 169 PetscFunctionReturn(0); 170 } 171 172 ierr = VecNormBegin(F,NORM_2,&fnorm);CHKERRQ(ierr); 173 ierr = VecNormEnd(F,NORM_2,&fnorm);CHKERRQ(ierr); 174 } else { 175 if (!snes->vec_func_init_set) { 176 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 177 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 178 if (domainerror) { 179 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 180 PetscFunctionReturn(0); 181 } 182 } else snes->vec_func_init_set = PETSC_FALSE; 183 184 if (!snes->norm_init_set) { 185 ierr = VecNormBegin(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 186 ierr = VecNormEnd(F,NORM_2,&fnorm);CHKERRQ(ierr); 187 if (PetscIsInfOrNanReal(fnorm)) { 188 snes->reason = SNES_DIVERGED_FNORM_NAN; 189 PetscFunctionReturn(0); 190 } 191 } else { 192 fnorm = snes->norm_init; 193 snes->norm_init_set = PETSC_FALSE; 194 } 195 } 196 197 ierr = PetscObjectAMSTakeAccess((PetscObject)snes);CHKERRQ(ierr); 198 snes->norm = fnorm; 199 ierr = PetscObjectAMSGrantAccess((PetscObject)snes);CHKERRQ(ierr); 200 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 201 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 202 203 /* set parameter for default relative tolerance convergence test */ 204 snes->ttol = fnorm*snes->rtol; 205 /* test convergence */ 206 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 207 if (snes->reason) PetscFunctionReturn(0); 208 209 for (i=0; i<maxits; i++) { 210 211 /* Call general purpose update function */ 212 if (snes->ops->update) { 213 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 214 } 215 216 /* apply the nonlinear preconditioner */ 217 if (snes->pc) { 218 if (snes->pcside == PC_RIGHT) { 219 ierr = SNESSetInitialFunction(snes->pc, F);CHKERRQ(ierr); 220 ierr = SNESSetInitialFunctionNorm(snes->pc, fnorm);CHKERRQ(ierr); 221 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,snes->vec_rhs,0);CHKERRQ(ierr); 222 ierr = SNESSolve(snes->pc, snes->vec_rhs, X);CHKERRQ(ierr); 223 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,snes->vec_rhs,0);CHKERRQ(ierr); 224 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 225 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 226 snes->reason = SNES_DIVERGED_INNER; 227 PetscFunctionReturn(0); 228 } 229 ierr = SNESGetPCFunction(snes,F,&fnorm);CHKERRQ(ierr); 230 } else if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 231 ierr = SNESApplyPC(snes,X,F,&fnorm,F);CHKERRQ(ierr); 232 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 233 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 234 snes->reason = SNES_DIVERGED_INNER; 235 PetscFunctionReturn(0); 236 } 237 } 238 } 239 240 /* Solve J Y = F, where J is Jacobian matrix */ 241 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 242 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 243 ierr = KSPSolve(snes->ksp,F,Y);CHKERRQ(ierr); 244 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 245 if (kspreason < 0) { 246 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 247 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 248 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 249 break; 250 } 251 } 252 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 253 snes->linear_its += lits; 254 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 255 256 if (PetscLogPrintInfo) { 257 ierr = SNESNEWTONLSCheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 258 } 259 260 /* Compute a (scaled) negative update in the line search routine: 261 X <- X - lambda*Y 262 and evaluate F = function(X) (depends on the line search). 263 */ 264 gnorm = fnorm; 265 ierr = SNESLineSearchApply(linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 266 ierr = SNESLineSearchGetSuccess(linesearch, &lssucceed);CHKERRQ(ierr); 267 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 268 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); 269 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 270 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 271 if (domainerror) { 272 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 273 PetscFunctionReturn(0); 274 } 275 if (!lssucceed) { 276 if (snes->stol*xnorm > ynorm) { 277 snes->reason = SNES_CONVERGED_SNORM_RELATIVE; 278 PetscFunctionReturn(0); 279 } 280 if (++snes->numFailures >= snes->maxFailures) { 281 PetscBool ismin; 282 snes->reason = SNES_DIVERGED_LINE_SEARCH; 283 ierr = SNESNEWTONLSCheckLocalMin_Private(snes,snes->jacobian,F,W,fnorm,&ismin);CHKERRQ(ierr); 284 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 285 break; 286 } 287 } 288 /* Monitor convergence */ 289 ierr = PetscObjectAMSTakeAccess((PetscObject)snes);CHKERRQ(ierr); 290 snes->iter = i+1; 291 snes->norm = fnorm; 292 ierr = PetscObjectAMSGrantAccess((PetscObject)snes);CHKERRQ(ierr); 293 ierr = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr); 294 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 295 /* Test for convergence */ 296 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 297 if (snes->reason) break; 298 } 299 if (i == maxits) { 300 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 301 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 302 } 303 PetscFunctionReturn(0); 304 } 305 /* -------------------------------------------------------------------------- */ 306 /* 307 SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use 308 of the SNESNEWTONLS nonlinear solver. 309 310 Input Parameter: 311 . snes - the SNES context 312 . x - the solution vector 313 314 Application Interface Routine: SNESSetUp() 315 316 Notes: 317 For basic use of the SNES solvers, the user need not explicitly call 318 SNESSetUp(), since these actions will automatically occur during 319 the call to SNESSolve(). 320 */ 321 #undef __FUNCT__ 322 #define __FUNCT__ "SNESSetUp_NEWTONLS" 323 PetscErrorCode SNESSetUp_NEWTONLS(SNES snes) 324 { 325 PetscErrorCode ierr; 326 327 PetscFunctionBegin; 328 ierr = SNESSetWorkVecs(snes,2);CHKERRQ(ierr); 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(SNES snes) 400 { 401 PetscErrorCode ierr; 402 SNESLineSearch linesearch; 403 404 PetscFunctionBegin; 405 ierr = PetscOptionsHead("SNESNEWTONLS options");CHKERRQ(ierr); 406 ierr = PetscOptionsTail();CHKERRQ(ierr); 407 /* set the default line search type */ 408 if (!snes->linesearch) { 409 ierr = SNESGetLineSearch(snes, &linesearch);CHKERRQ(ierr); 410 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);CHKERRQ(ierr); 411 } 412 PetscFunctionReturn(0); 413 } 414 415 /* -------------------------------------------------------------------------- */ 416 /*MC 417 SNESNEWTONLS - Newton based nonlinear solver that uses a line search 418 419 Options Database: 420 + -snes_linesearch_type <bt> - bt,basic. Select line search type 421 . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt 422 . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch 423 . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 424 . -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) 425 . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate 426 . -snes_linesearch_monitor - print information about progress of line searches 427 - -snes_linesearch_damping - damping factor used for basic line search 428 429 Notes: This is the default nonlinear solver in SNES 430 431 Level: beginner 432 433 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONTR, SNESQN, SNESLineSearchSetType(), SNESLineSearchSetOrder() 434 SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck() SNESLineSearchSetComputeNorms() 435 436 M*/ 437 #undef __FUNCT__ 438 #define __FUNCT__ "SNESCreate_NEWTONLS" 439 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONLS(SNES snes) 440 { 441 PetscErrorCode ierr; 442 SNES_NEWTONLS *neP; 443 444 PetscFunctionBegin; 445 snes->ops->setup = SNESSetUp_NEWTONLS; 446 snes->ops->solve = SNESSolve_NEWTONLS; 447 snes->ops->destroy = SNESDestroy_NEWTONLS; 448 snes->ops->setfromoptions = SNESSetFromOptions_NEWTONLS; 449 snes->ops->view = SNESView_NEWTONLS; 450 snes->ops->reset = SNESReset_NEWTONLS; 451 452 snes->pcside = PC_RIGHT; 453 snes->usesksp = PETSC_TRUE; 454 snes->usespc = PETSC_TRUE; 455 ierr = PetscNewLog(snes,SNES_NEWTONLS,&neP);CHKERRQ(ierr); 456 snes->data = (void*)neP; 457 PetscFunctionReturn(0); 458 } 459