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 SNES_LS *ls; 142 143 PetscFunctionBegin; 144 ls = (SNES_LS *)snes->data; 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->work[0]; /* work vectors */ 153 G = snes->work[1]; 154 W = snes->work[2]; 155 156 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 157 snes->iter = 0; 158 snes->norm = 0.0; 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 (snes->ops->precheckstep) { 207 PetscBool changed_y = PETSC_FALSE; 208 ierr = (*snes->ops->precheckstep)(snes,X,Y,snes->precheck,&changed_y);CHKERRQ(ierr); 209 } 210 211 if (PetscLogPrintInfo){ 212 ierr = SNESLSCheckResidual_Private(snes,snes->jacobian,F,Y,G,W);CHKERRQ(ierr); 213 } 214 215 /* Compute a (scaled) negative update in the line search routine: 216 Y <- X - lambda*Y 217 and evaluate G = function(Y) (depends on the line search). 218 */ 219 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 220 gnorm = fnorm; 221 ierr = PetscLineSearchApply(ls->linesearch, X, F, &fnorm, Y);CHKERRQ(ierr); 222 ierr = PetscLineSearchGetSuccess(ls->linesearch, &lssucceed);CHKERRQ(ierr); 223 ierr = PetscLineSearchGetNorms(ls->linesearch, &xnorm, &fnorm, &ynorm);CHKERRQ(ierr); 224 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); 225 if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break; 226 ierr = SNESGetFunctionDomainError(snes, &domainerror);CHKERRQ(ierr); 227 if (domainerror) { 228 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 229 PetscFunctionReturn(0); 230 } 231 if (!lssucceed) { 232 if (++snes->numFailures >= snes->maxFailures) { 233 PetscBool ismin; 234 snes->reason = SNES_DIVERGED_LINE_SEARCH; 235 ierr = SNESLSCheckLocalMin_Private(snes,snes->jacobian,G,W,gnorm,&ismin);CHKERRQ(ierr); 236 if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN; 237 break; 238 } 239 } 240 /* Monitor convergence */ 241 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 242 snes->iter = i+1; 243 snes->norm = fnorm; 244 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 245 SNESLogConvHistory(snes,snes->norm,lits); 246 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 247 /* Test for convergence, xnorm = || X || */ 248 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 249 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 250 if (snes->reason) break; 251 } 252 if (i == maxits) { 253 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 254 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 255 } 256 PetscFunctionReturn(0); 257 } 258 /* -------------------------------------------------------------------------- */ 259 /* 260 SNESSetUp_LS - Sets up the internal data structures for the later use 261 of the SNESLS nonlinear solver. 262 263 Input Parameter: 264 . snes - the SNES context 265 . x - the solution vector 266 267 Application Interface Routine: SNESSetUp() 268 269 Notes: 270 For basic use of the SNES solvers, the user need not explicitly call 271 SNESSetUp(), since these actions will automatically occur during 272 the call to SNESSolve(). 273 */ 274 #undef __FUNCT__ 275 #define __FUNCT__ "SNESSetUp_LS" 276 PetscErrorCode SNESSetUp_LS(SNES snes) 277 { 278 PetscErrorCode ierr; 279 SNES_LS *ls; 280 const char *optionsprefix; 281 282 PetscFunctionBegin; 283 ierr = SNESDefaultGetWork(snes,3);CHKERRQ(ierr); 284 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 285 ls = (SNES_LS *)snes->data; 286 287 /* set up the line search */ 288 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 289 ierr = PetscLineSearchCreate(((PetscObject)snes)->comm, &ls->linesearch);CHKERRQ(ierr); 290 ierr = PetscLineSearchSetSNES(ls->linesearch, snes);CHKERRQ(ierr); 291 ierr = PetscLineSearchSetType(ls->linesearch, PETSCLINESEARCHBT);CHKERRQ(ierr); 292 293 ierr = PetscLineSearchAppendOptionsPrefix(ls->linesearch, optionsprefix);CHKERRQ(ierr); 294 ierr = PetscLineSearchSetFromOptions(ls->linesearch);CHKERRQ(ierr); 295 296 PetscFunctionReturn(0); 297 } 298 /* -------------------------------------------------------------------------- */ 299 300 #undef __FUNCT__ 301 #define __FUNCT__ "SNESReset_LS" 302 PetscErrorCode SNESReset_LS(SNES snes) 303 { 304 SNES_LS *ls; 305 PetscErrorCode ierr; 306 PetscFunctionBegin; 307 ls = (SNES_LS *)snes->data; 308 ierr = PetscLineSearchDestroy(&ls->linesearch);CHKERRQ(ierr); 309 PetscFunctionReturn(0); 310 } 311 312 /* 313 SNESDestroy_LS - Destroys the private SNES_LS context that was created 314 with SNESCreate_LS(). 315 316 Input Parameter: 317 . snes - the SNES context 318 319 Application Interface Routine: SNESDestroy() 320 */ 321 #undef __FUNCT__ 322 #define __FUNCT__ "SNESDestroy_LS" 323 PetscErrorCode SNESDestroy_LS(SNES snes) 324 { 325 PetscErrorCode ierr; 326 327 PetscFunctionBegin; 328 ierr = SNESReset_LS(snes);CHKERRQ(ierr); 329 ierr = PetscFree(snes->data);CHKERRQ(ierr); 330 PetscFunctionReturn(0); 331 } 332 /* -------------------------------------------------------------------------- */ 333 334 /* 335 SNESView_LS - Prints info from the SNESLS data structure. 336 337 Input Parameters: 338 . SNES - the SNES context 339 . viewer - visualization context 340 341 Application Interface Routine: SNESView() 342 */ 343 #undef __FUNCT__ 344 #define __FUNCT__ "SNESView_LS" 345 static PetscErrorCode SNESView_LS(SNES snes,PetscViewer viewer) 346 { 347 const char *cstr; 348 PetscErrorCode ierr; 349 PetscBool iascii; 350 351 PetscFunctionBegin; 352 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 353 if (iascii) { 354 cstr = SNESLineSearchTypeName(snes->ls_type); 355 ierr = PetscViewerASCIIPrintf(viewer," line search variant: %s\n",cstr);CHKERRQ(ierr); 356 ierr = PetscViewerASCIIPrintf(viewer," alpha=%14.12e, maxstep=%14.12e, minlambda=%14.12e\n",(double)snes->ls_alpha,(double)snes->maxstep,(double)snes->steptol);CHKERRQ(ierr); 357 ierr = PetscViewerASCIIPrintf(viewer," damping factor=%14.12e\n",(double)snes->damping);CHKERRQ(ierr); 358 } 359 PetscFunctionReturn(0); 360 } 361 362 /* -------------------------------------------------------------------------- */ 363 /* 364 SNESSetFromOptions_LS - Sets various parameters for the SNESLS method. 365 366 Input Parameter: 367 . snes - the SNES context 368 369 Application Interface Routine: SNESSetFromOptions() 370 */ 371 #undef __FUNCT__ 372 #define __FUNCT__ "SNESSetFromOptions_LS" 373 static PetscErrorCode SNESSetFromOptions_LS(SNES snes) 374 { 375 PetscErrorCode ierr; 376 377 PetscFunctionBegin; 378 ierr = PetscOptionsHead("SNES Line search options");CHKERRQ(ierr); 379 ierr = PetscOptionsTail();CHKERRQ(ierr); 380 PetscFunctionReturn(0); 381 } 382 383 /* -------------------------------------------------------------------------- */ 384 /*MC 385 SNESLS - Newton based nonlinear solver that uses a line search 386 387 Options Database: 388 + -snes_ls [cubic,quadratic,basic,basicnonorms] - Selects line search 389 . -snes_ls_alpha <alpha> - Sets alpha used in determining if reduction in function norm is sufficient 390 . -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) 391 . -snes_ls_minlambda <minlambda> - Sets the minimum lambda the line search will use minlambda / max_i ( y[i]/x[i] ) 392 . -snes_ls_monitor - print information about progress of line searches 393 - -snes_ls_damping - damping factor used if -snes_ls is basic or basicnonorms 394 395 396 Notes: This is the default nonlinear solver in SNES 397 398 Level: beginner 399 400 .seealso: SNESCreate(), SNES, SNESSetType(), SNESTR, SNESLineSearchSet(), 401 SNESLineSearchSetPostCheck(), SNESLineSearchNo(), SNESLineSearchCubic(), SNESLineSearchQuadratic(), 402 SNESLineSearchSet(), SNESLineSearchNoNorms(), SNESLineSearchSetPreCheck(), SNESLineSearchSetParams(), SNESLineSearchGetParams() 403 404 M*/ 405 EXTERN_C_BEGIN 406 #undef __FUNCT__ 407 #define __FUNCT__ "SNESCreate_LS" 408 PetscErrorCode SNESCreate_LS(SNES snes) 409 { 410 PetscErrorCode ierr; 411 SNES_LS *neP; 412 413 PetscFunctionBegin; 414 snes->ops->setup = SNESSetUp_LS; 415 snes->ops->solve = SNESSolve_LS; 416 snes->ops->destroy = SNESDestroy_LS; 417 snes->ops->setfromoptions = SNESSetFromOptions_LS; 418 snes->ops->view = SNESView_LS; 419 snes->ops->reset = SNESReset_LS; 420 421 snes->usesksp = PETSC_TRUE; 422 snes->usespc = PETSC_FALSE; 423 ierr = PetscNewLog(snes,SNES_LS,&neP);CHKERRQ(ierr); 424 snes->data = (void*)neP; 425 PetscFunctionReturn(0); 426 } 427 EXTERN_C_END 428