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