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->work[0]; /* work vectors */ 152 G = snes->work[1]; 153 W = snes->work[2]; 154 155 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 156 snes->iter = 0; 157 snes->norm = 0.0; 158 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 159 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 160 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 161 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 162 if (domainerror) { 163 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 164 PetscFunctionReturn(0); 165 } 166 ierr = VecNormBegin(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 167 ierr = VecNormBegin(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm <- ||x|| */ 168 ierr = VecNormEnd(F,NORM_2,&fnorm);CHKERRQ(ierr); 169 ierr = VecNormEnd(X,NORM_2,&xnorm);CHKERRQ(ierr); 170 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"User provided compute function generated a Not-a-Number"); 171 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 172 snes->norm = fnorm; 173 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 174 SNESLogConvHistory(snes,fnorm,0); 175 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 176 177 /* set parameter for default relative tolerance convergence test */ 178 snes->ttol = fnorm*snes->rtol; 179 /* test convergence */ 180 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 181 if (snes->reason) PetscFunctionReturn(0); 182 183 for (i=0; i<maxits; i++) { 184 185 /* Call general purpose update function */ 186 if (snes->ops->update) { 187 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 188 } 189 190 /* Solve J Y = F, where J is Jacobian matrix */ 191 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 192 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 193 ierr = SNES_KSPSolve(snes,snes->ksp,F,Y);CHKERRQ(ierr); 194 ierr = KSPGetConvergedReason(snes->ksp,&kspreason);CHKERRQ(ierr); 195 if (kspreason < 0) { 196 if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) { 197 ierr = PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);CHKERRQ(ierr); 198 snes->reason = SNES_DIVERGED_LINEAR_SOLVE; 199 break; 200 } 201 } 202 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 203 snes->linear_its += lits; 204 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 205 206 if (PetscLogPrintInfo){ 207 ierr = SNESLSCheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 208 } 209 210 /* Compute a (scaled) negative update in the line search routine: 211 Y <- X - lambda*Y 212 and evaluate G = function(Y) (depends on the line search). 213 */ 214 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 215 gnorm = fnorm; 216 ierr = SNESLineSearchApply(linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 217 ierr = SNESLineSearchGetSuccess(linesearch, &lssucceed);CHKERRQ(ierr); 218 ierr = SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 219 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); 220 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 221 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 222 if (domainerror) { 223 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 224 PetscFunctionReturn(0); 225 } 226 if (!lssucceed) { 227 if (++snes->numFailures >= snes->maxFailures) { 228 PetscBool ismin; 229 snes->reason = SNES_DIVERGED_LINE_SEARCH; 230 ierr = SNESLSCheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 231 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 232 break; 233 } 234 } 235 /* Monitor convergence */ 236 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 237 snes->iter = i+1; 238 snes->norm = fnorm; 239 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 240 SNESLogConvHistory(snes,snes->norm,lits); 241 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 242 /* Test for convergence, xnorm = || X || */ 243 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 244 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 245 if (snes->reason) break; 246 } 247 if (i == maxits) { 248 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 249 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 250 } 251 PetscFunctionReturn(0); 252 } 253 /* -------------------------------------------------------------------------- */ 254 /* 255 SNESSetUp_LS - Sets up the internal data structures for the later use 256 of the SNESLS nonlinear solver. 257 258 Input Parameter: 259 . snes - the SNES context 260 . x - the solution vector 261 262 Application Interface Routine: SNESSetUp() 263 264 Notes: 265 For basic use of the SNES solvers, the user need not explicitly call 266 SNESSetUp(), since these actions will automatically occur during 267 the call to SNESSolve(). 268 */ 269 #undef __FUNCT__ 270 #define __FUNCT__ "SNESSetUp_LS" 271 PetscErrorCode SNESSetUp_LS(SNES snes) 272 { 273 PetscErrorCode ierr; 274 275 PetscFunctionBegin; 276 ierr = SNESDefaultGetWork(snes,3);CHKERRQ(ierr); 277 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 278 279 PetscFunctionReturn(0); 280 } 281 /* -------------------------------------------------------------------------- */ 282 283 #undef __FUNCT__ 284 #define __FUNCT__ "SNESReset_LS" 285 PetscErrorCode SNESReset_LS(SNES snes) 286 { 287 PetscFunctionBegin; 288 PetscFunctionReturn(0); 289 } 290 291 /* 292 SNESDestroy_LS - Destroys the private SNES_LS context that was created 293 with SNESCreate_LS(). 294 295 Input Parameter: 296 . snes - the SNES context 297 298 Application Interface Routine: SNESDestroy() 299 */ 300 #undef __FUNCT__ 301 #define __FUNCT__ "SNESDestroy_LS" 302 PetscErrorCode SNESDestroy_LS(SNES snes) 303 { 304 PetscErrorCode ierr; 305 306 PetscFunctionBegin; 307 ierr = SNESReset_LS(snes);CHKERRQ(ierr); 308 ierr = PetscFree(snes->data);CHKERRQ(ierr); 309 PetscFunctionReturn(0); 310 } 311 /* -------------------------------------------------------------------------- */ 312 313 /* 314 SNESView_LS - Prints info from the SNESLS data structure. 315 316 Input Parameters: 317 . SNES - the SNES context 318 . viewer - visualization context 319 320 Application Interface Routine: SNESView() 321 */ 322 #undef __FUNCT__ 323 #define __FUNCT__ "SNESView_LS" 324 static PetscErrorCode SNESView_LS(SNES snes,PetscViewer viewer) 325 { 326 PetscErrorCode ierr; 327 PetscBool iascii; 328 329 PetscFunctionBegin; 330 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 331 if (iascii) { 332 } 333 PetscFunctionReturn(0); 334 } 335 336 /* -------------------------------------------------------------------------- */ 337 /* 338 SNESSetFromOptions_LS - Sets various parameters for the SNESLS method. 339 340 Input Parameter: 341 . snes - the SNES context 342 343 Application Interface Routine: SNESSetFromOptions() 344 */ 345 #undef __FUNCT__ 346 #define __FUNCT__ "SNESSetFromOptions_LS" 347 static PetscErrorCode SNESSetFromOptions_LS(SNES snes) 348 { 349 PetscErrorCode ierr; 350 SNESLineSearch linesearch; 351 352 PetscFunctionBegin; 353 ierr = PetscOptionsHead("SNESLS options");CHKERRQ(ierr); 354 ierr = PetscOptionsTail();CHKERRQ(ierr); 355 /* set the default line search type */ 356 if (!snes->linesearch) { 357 ierr = SNESGetSNESLineSearch(snes, &linesearch);CHKERRQ(ierr); 358 ierr = SNESLineSearchSetType(linesearch, SNES_LINESEARCH_BT);CHKERRQ(ierr); 359 } 360 PetscFunctionReturn(0); 361 } 362 363 /* -------------------------------------------------------------------------- */ 364 /*MC 365 SNESLS - Newton based nonlinear solver that uses a line search 366 367 Options Database: 368 + -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search 369 . -snes_ls_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 370 . -snes_ls_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) 371 . -snes_ls_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 372 . -snes_ls_monitor - print information about progress of line searches 373 - -snes_ls_damping - damping factor used if -snes_ls is basic or basicnonorms 374 375 376 Notes: This is the default nonlinear solver in SNES 377 378 Level: beginner 379 380 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 381 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 382 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 383 384 M*/ 385 EXTERN_C_BEGIN 386 #undef __FUNCT__ 387 #define __FUNCT__ "SNESCreate_LS" 388 PetscErrorCode SNESCreate_LS(SNES snes) 389 { 390 PetscErrorCode ierr; 391 SNES_LS *neP; 392 393 PetscFunctionBegin; 394 snes->ops->setup = SNESSetUp_LS; 395 snes->ops->solve = SNESSolve_LS; 396 snes->ops->destroy = SNESDestroy_LS; 397 snes->ops->setfromoptions = SNESSetFromOptions_LS; 398 snes->ops->view = SNESView_LS; 399 snes->ops->reset = SNESReset_LS; 400 401 snes->usesksp = PETSC_TRUE; 402 snes->usespc = PETSC_FALSE; 403 ierr = PetscNewLog(snes,SNES_LS,&neP);CHKERRQ(ierr); 404 snes->data = (void*)neP; 405 PetscFunctionReturn(0); 406 } 407 EXTERN_C_END 408