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