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