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__ "SNESLSCheckLocalMin_Private" 12 PetscErrorCode SNESLSCheckLocalMin_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__ "SNESLSCheckResidual_Private" 50 PetscErrorCode SNESLSCheckResidual_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 _LS (e.g., SNESCreate_LS, SNESSolve_LS) 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_LS - 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_LS" 131 PetscErrorCode SNESSolve_LS(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; 139 KSPConvergedReason kspreason; 140 PetscBool domainerror; 141 SNESLineSearch linesearch; 142 143 PetscFunctionBegin; 144 snes->numFailures = 0; 145 snes->numLinearSolveFailures = 0; 146 snes->reason = SNES_CONVERGED_ITERATING; 147 148 maxits = snes->max_its; /* maximum number of iterations */ 149 X = snes->vec_sol; /* solution vector */ 150 F = snes->vec_func; /* residual vector */ 151 Y = snes->vec_sol_update; /* newton step */ 152 G = snes->work[0]; 153 W = snes->work[1]; 154 155 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 156 snes->iter = 0; 157 snes->norm = 0.0; 158 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 159 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 160 if (!snes->vec_func_init_set) { 161 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 162 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 163 if (domainerror) { 164 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 165 PetscFunctionReturn(0); 166 } 167 } else { 168 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 /* Solve J Y = F, where J is Jacobian matrix */ 198 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 199 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 200 ierr = SNES_KSPSolve(snes,snes->ksp,F,Y);CHKERRQ(ierr); 201 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 202 if (kspreason < 0) { 203 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 204 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 205 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 206 break; 207 } 208 } 209 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 210 snes->linear_its += lits; 211 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 212 213 if (PetscLogPrintInfo){ 214 ierr = SNESLSCheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 215 } 216 217 /* Compute a (scaled) negative update in the line search routine: 218 X <- X - lambda*Y 219 and evaluate F = function(X) (depends on the line search). 220 */ 221 gnorm = fnorm; 222 ierr = SNESLineSearchApply(linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 223 ierr = SNESLineSearchGetSuccess(linesearch, &lssucceed);CHKERRQ(ierr); 224 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 225 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); 226 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 227 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 228 if (domainerror) { 229 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 230 PetscFunctionReturn(0); 231 } 232 if (!lssucceed) { 233 if (snes->stol*xnorm > ynorm) { 234 snes->reason = SNES_CONVERGED_SNORM_RELATIVE; 235 PetscFunctionReturn(0); 236 } 237 if (++snes->numFailures >= snes->maxFailures) { 238 PetscBool ismin; 239 snes->reason = SNES_DIVERGED_LINE_SEARCH; 240 ierr = SNESLSCheckLocalMin_Private(snes,snes->jacobian,F,X,fnorm,&ismin);CHKERRQ(ierr); 241 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 242 break; 243 } 244 } 245 /* Monitor convergence */ 246 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 247 snes->iter = i+1; 248 snes->norm = fnorm; 249 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 250 SNESLogConvHistory(snes,snes->norm,lits); 251 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 252 /* Test for convergence, xnorm = || X || */ 253 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 254 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 255 if (snes->reason) break; 256 } 257 if (i == maxits) { 258 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 259 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 260 } 261 PetscFunctionReturn(0); 262 } 263 /* -------------------------------------------------------------------------- */ 264 /* 265 SNESSetUp_LS - Sets up the internal data structures for the later use 266 of the SNESLS nonlinear solver. 267 268 Input Parameter: 269 . snes - the SNES context 270 . x - the solution vector 271 272 Application Interface Routine: SNESSetUp() 273 274 Notes: 275 For basic use of the SNES solvers, the user need not explicitly call 276 SNESSetUp(), since these actions will automatically occur during 277 the call to SNESSolve(). 278 */ 279 #undef __FUNCT__ 280 #define __FUNCT__ "SNESSetUp_LS" 281 PetscErrorCode SNESSetUp_LS(SNES snes) 282 { 283 PetscErrorCode ierr; 284 285 PetscFunctionBegin; 286 ierr = SNESDefaultGetWork(snes,2);CHKERRQ(ierr); 287 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 288 289 PetscFunctionReturn(0); 290 } 291 /* -------------------------------------------------------------------------- */ 292 293 #undef __FUNCT__ 294 #define __FUNCT__ "SNESReset_LS" 295 PetscErrorCode SNESReset_LS(SNES snes) 296 { 297 PetscFunctionBegin; 298 PetscFunctionReturn(0); 299 } 300 301 /* 302 SNESDestroy_LS - Destroys the private SNES_LS context that was created 303 with SNESCreate_LS(). 304 305 Input Parameter: 306 . snes - the SNES context 307 308 Application Interface Routine: SNESDestroy() 309 */ 310 #undef __FUNCT__ 311 #define __FUNCT__ "SNESDestroy_LS" 312 PetscErrorCode SNESDestroy_LS(SNES snes) 313 { 314 PetscErrorCode ierr; 315 316 PetscFunctionBegin; 317 ierr = SNESReset_LS(snes);CHKERRQ(ierr); 318 ierr = PetscFree(snes->data);CHKERRQ(ierr); 319 PetscFunctionReturn(0); 320 } 321 /* -------------------------------------------------------------------------- */ 322 323 /* 324 SNESView_LS - Prints info from the SNESLS data structure. 325 326 Input Parameters: 327 . SNES - the SNES context 328 . viewer - visualization context 329 330 Application Interface Routine: SNESView() 331 */ 332 #undef __FUNCT__ 333 #define __FUNCT__ "SNESView_LS" 334 static PetscErrorCode SNESView_LS(SNES snes,PetscViewer viewer) 335 { 336 PetscErrorCode ierr; 337 PetscBool iascii; 338 339 PetscFunctionBegin; 340 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 341 if (iascii) { 342 } 343 PetscFunctionReturn(0); 344 } 345 346 /* -------------------------------------------------------------------------- */ 347 /* 348 SNESSetFromOptions_LS - Sets various parameters for the SNESLS method. 349 350 Input Parameter: 351 . snes - the SNES context 352 353 Application Interface Routine: SNESSetFromOptions() 354 */ 355 #undef __FUNCT__ 356 #define __FUNCT__ "SNESSetFromOptions_LS" 357 static PetscErrorCode SNESSetFromOptions_LS(SNES snes) 358 { 359 PetscErrorCode ierr; 360 SNESLineSearch linesearch; 361 362 PetscFunctionBegin; 363 ierr = PetscOptionsHead("SNESLS options");CHKERRQ(ierr); 364 ierr = PetscOptionsTail();CHKERRQ(ierr); 365 /* set the default line search type */ 366 if (!snes->linesearch) { 367 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 368 ierr = SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);CHKERRQ(ierr); 369 } 370 PetscFunctionReturn(0); 371 } 372 373 /* -------------------------------------------------------------------------- */ 374 /*MC 375 SNESLS - Newton based nonlinear solver that uses a line search 376 377 Options Database: 378 + -snes_linesearch_type <bt> - bt,basic. Select line search type 379 . -snes_linesearch_order <3> - 2, 3. Selects the order of the line search for bt 380 . -snes_linesearch_norms <true> - Turns on/off computation of the norms for basic linesearch 381 . -snes_linesearch_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 382 . -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) 383 . -snes_linesearch_minlambda <minlambda> - Sets the minimum lambda the line search will tolerate 384 . -snes_linesearch_monitor - print information about progress of line searches 385 - -snes_linesearch_damping - damping factor used for basic line search 386 387 Notes: This is the default nonlinear solver in SNES 388 389 Level: beginner 390 391 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESQN, SNESLineSearchSetType(), SNESLineSearchSetOrder() 392 SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck() SNESLineSearchSetComputeNorms() 393 394 M*/ 395 EXTERN_C_BEGIN 396 #undef __FUNCT__ 397 #define __FUNCT__ "SNESCreate_LS" 398 PetscErrorCode SNESCreate_LS(SNES snes) 399 { 400 PetscErrorCode ierr; 401 SNES_LS *neP; 402 403 PetscFunctionBegin; 404 snes->ops->setup = SNESSetUp_LS; 405 snes->ops->solve = SNESSolve_LS; 406 snes->ops->destroy = SNESDestroy_LS; 407 snes->ops->setfromoptions = SNESSetFromOptions_LS; 408 snes->ops->view = SNESView_LS; 409 snes->ops->reset = SNESReset_LS; 410 411 snes->usesksp = PETSC_TRUE; 412 snes->usespc = PETSC_FALSE; 413 ierr = PetscNewLog(snes,SNES_LS,&neP);CHKERRQ(ierr); 414 snes->data = (void*)neP; 415 PetscFunctionReturn(0); 416 } 417 EXTERN_C_END 418