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