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 184 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 185 if (PetscIsInfOrNanReal(fnorm)) { 186 snes->reason = SNES_DIVERGED_FNORM_NAN; 187 PetscFunctionReturn(0); 188 } 189 190 ierr = PetscObjectAMSTakeAccess((PetscObject)snes);CHKERRQ(ierr); 191 snes->norm = fnorm; 192 ierr = PetscObjectAMSGrantAccess((PetscObject)snes);CHKERRQ(ierr); 193 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 194 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 195 196 /* set parameter for default relative tolerance convergence test */ 197 snes->ttol = fnorm*snes->rtol; 198 /* test convergence */ 199 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 200 if (snes->reason) PetscFunctionReturn(0); 201 202 for (i=0; i<maxits; i++) { 203 204 /* Call general purpose update function */ 205 if (snes->ops->update) { 206 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 207 } 208 209 /* apply the nonlinear preconditioner */ 210 if (snes->pc) { 211 if (snes->pcside == PC_RIGHT) { 212 ierr = SNESSetInitialFunction(snes->pc, F);CHKERRQ(ierr); 213 ierr = PetscLogEventBegin(SNES_NPCSolve,snes->pc,X,snes->vec_rhs,0);CHKERRQ(ierr); 214 ierr = SNESSolve(snes->pc, snes->vec_rhs, X);CHKERRQ(ierr); 215 ierr = PetscLogEventEnd(SNES_NPCSolve,snes->pc,X,snes->vec_rhs,0);CHKERRQ(ierr); 216 ierr = SNESGetConvergedReason(snes->pc,&reason);CHKERRQ(ierr); 217 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) { 218 snes->reason = SNES_DIVERGED_INNER; 219 PetscFunctionReturn(0); 220 } 221 ierr = SNESGetPCFunction(snes,F,&fnorm);CHKERRQ(ierr); 222 } else if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) { 223 ierr = SNESApplyPC(snes,X,F,&fnorm,F);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 } 230 } 231 232 /* Solve J Y = F, where J is Jacobian matrix */ 233 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 234 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 235 ierr = KSPSolve(snes->ksp,F,Y);CHKERRQ(ierr); 236 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 237 if (kspreason < 0) { 238 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 239 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 240 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 241 break; 242 } 243 } 244 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 245 snes->linear_its += lits; 246 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 247 248 if (PetscLogPrintInfo) { 249 ierr = SNESNEWTONLSCheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 250 } 251 252 /* Compute a (scaled) negative update in the line search routine: 253 X <- X - lambda*Y 254 and evaluate F = function(X) (depends on the line search). 255 */ 256 gnorm = fnorm; 257 ierr = SNESLineSearchApply(linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 258 ierr = SNESLineSearchGetSuccess(linesearch, &lssucceed);CHKERRQ(ierr); 259 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 260 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); 261 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 262 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 263 if (domainerror) { 264 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 265 PetscFunctionReturn(0); 266 } 267 if (!lssucceed) { 268 if (snes->stol*xnorm > ynorm) { 269 snes->reason = SNES_CONVERGED_SNORM_RELATIVE; 270 PetscFunctionReturn(0); 271 } 272 if (++snes->numFailures >= snes->maxFailures) { 273 PetscBool ismin; 274 snes->reason = SNES_DIVERGED_LINE_SEARCH; 275 ierr = SNESNEWTONLSCheckLocalMin_Private(snes,snes->jacobian,F,W,fnorm,&ismin);CHKERRQ(ierr); 276 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 277 break; 278 } 279 } 280 /* Monitor convergence */ 281 ierr = PetscObjectAMSTakeAccess((PetscObject)snes);CHKERRQ(ierr); 282 snes->iter = i+1; 283 snes->norm = fnorm; 284 ierr = PetscObjectAMSGrantAccess((PetscObject)snes);CHKERRQ(ierr); 285 ierr = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr); 286 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 287 /* Test for convergence */ 288 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 289 if (snes->reason) break; 290 } 291 if (i == maxits) { 292 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 293 if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 294 } 295 PetscFunctionReturn(0); 296 } 297 /* -------------------------------------------------------------------------- */ 298 /* 299 SNESSetUp_NEWTONLS - Sets up the internal data structures for the later use 300 of the SNESNEWTONLS nonlinear solver. 301 302 Input Parameter: 303 . snes - the SNES context 304 . x - the solution vector 305 306 Application Interface Routine: SNESSetUp() 307 308 Notes: 309 For basic use of the SNES solvers, the user need not explicitly call 310 SNESSetUp(), since these actions will automatically occur during 311 the call to SNESSolve(). 312 */ 313 #undef __FUNCT__ 314 #define __FUNCT__ "SNESSetUp_NEWTONLS" 315 PetscErrorCode SNESSetUp_NEWTONLS(SNES snes) 316 { 317 PetscErrorCode ierr; 318 319 PetscFunctionBegin; 320 ierr = SNESSetWorkVecs(snes,2);CHKERRQ(ierr); 321 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 322 if (snes->pcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED; 323 PetscFunctionReturn(0); 324 } 325 /* -------------------------------------------------------------------------- */ 326 327 #undef __FUNCT__ 328 #define __FUNCT__ "SNESReset_NEWTONLS" 329 PetscErrorCode SNESReset_NEWTONLS(SNES snes) 330 { 331 PetscFunctionBegin; 332 PetscFunctionReturn(0); 333 } 334 335 /* 336 SNESDestroy_NEWTONLS - Destroys the private SNES_NEWTONLS context that was created 337 with SNESCreate_NEWTONLS(). 338 339 Input Parameter: 340 . snes - the SNES context 341 342 Application Interface Routine: SNESDestroy() 343 */ 344 #undef __FUNCT__ 345 #define __FUNCT__ "SNESDestroy_NEWTONLS" 346 PetscErrorCode SNESDestroy_NEWTONLS(SNES snes) 347 { 348 PetscErrorCode ierr; 349 350 PetscFunctionBegin; 351 ierr = SNESReset_NEWTONLS(snes);CHKERRQ(ierr); 352 ierr = PetscFree(snes->data);CHKERRQ(ierr); 353 PetscFunctionReturn(0); 354 } 355 /* -------------------------------------------------------------------------- */ 356 357 /* 358 SNESView_NEWTONLS - Prints info from the SNESNEWTONLS data structure. 359 360 Input Parameters: 361 . SNES - the SNES context 362 . viewer - visualization context 363 364 Application Interface Routine: SNESView() 365 */ 366 #undef __FUNCT__ 367 #define __FUNCT__ "SNESView_NEWTONLS" 368 static PetscErrorCode SNESView_NEWTONLS(SNES snes,PetscViewer viewer) 369 { 370 PetscErrorCode ierr; 371 PetscBool iascii; 372 373 PetscFunctionBegin; 374 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 375 if (iascii) { 376 } 377 PetscFunctionReturn(0); 378 } 379 380 /* -------------------------------------------------------------------------- */ 381 /* 382 SNESSetFromOptions_NEWTONLS - Sets various parameters for the SNESNEWTONLS method. 383 384 Input Parameter: 385 . snes - the SNES context 386 387 Application Interface Routine: SNESSetFromOptions() 388 */ 389 #undef __FUNCT__ 390 #define __FUNCT__ "SNESSetFromOptions_NEWTONLS" 391 static PetscErrorCode SNESSetFromOptions_NEWTONLS(SNES snes) 392 { 393 PetscErrorCode ierr; 394 SNESLineSearch linesearch; 395 396 PetscFunctionBegin; 397 ierr = PetscOptionsHead("SNESNEWTONLS options");CHKERRQ(ierr); 398 ierr = PetscOptionsTail();CHKERRQ(ierr); 399 /* set the default line search type */ 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,SNES_NEWTONLS,&neP);CHKERRQ(ierr); 448 snes->data = (void*)neP; 449 PetscFunctionReturn(0); 450 } 451