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