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