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