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