1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: ls.c,v 1.94 1997/08/03 16:36:19 curfman Exp bsmith $"; 3 #endif 4 5 #include <math.h> 6 #include "src/snes/impls/ls/ls.h" 7 #include "pinclude/pviewer.h" 8 9 /* 10 Implements a line search variant of Newton's Method 11 for solving systems of nonlinear equations. 12 13 Input parameters: 14 . snes - nonlinear context obtained from SNESCreate() 15 16 Output Parameters: 17 . outits - Number of global iterations until termination. 18 19 Notes: 20 This implements essentially a truncated Newton method with a 21 line search. By default a cubic backtracking line search 22 is employed, as described in the text "Numerical Methods for 23 Unconstrained Optimization and Nonlinear Equations" by Dennis 24 and Schnabel. 25 */ 26 27 #undef __FUNC__ 28 #define __FUNC__ "SNESSolve_EQ_LS" 29 int SNESSolve_EQ_LS(SNES snes,int *outits) 30 { 31 SNES_LS *neP = (SNES_LS *) snes->data; 32 int maxits, i, history_len, ierr, lits, lsfail; 33 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 34 double fnorm, gnorm, xnorm, ynorm, *history; 35 Vec Y, X, F, G, W, TMP; 36 37 history = snes->conv_hist; /* convergence history */ 38 history_len = snes->conv_hist_size; /* convergence history length */ 39 maxits = snes->max_its; /* maximum number of iterations */ 40 X = snes->vec_sol; /* solution vector */ 41 F = snes->vec_func; /* residual vector */ 42 Y = snes->work[0]; /* work vectors */ 43 G = snes->work[1]; 44 W = snes->work[2]; 45 46 snes->iter = 0; 47 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* F(X) */ 48 ierr = VecNorm(F,NORM_2,&fnorm); CHKERRQ(ierr); /* fnorm <- ||F|| */ 49 snes->norm = fnorm; 50 if (history) history[0] = fnorm; 51 SNESMonitor(snes,0,fnorm); 52 53 if (fnorm < snes->atol) {*outits = 0; return 0;} 54 55 /* set parameter for default relative tolerance convergence test */ 56 snes->ttol = fnorm*snes->rtol; 57 58 for ( i=0; i<maxits; i++ ) { 59 snes->iter = i+1; 60 61 /* Solve J Y = F, where J is Jacobian matrix */ 62 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg); CHKERRQ(ierr); 63 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg); CHKERRQ(ierr); 64 ierr = SLESSolve(snes->sles,F,Y,&lits); CHKERRQ(ierr); 65 snes->linear_its += PetscAbsInt(lits); 66 PLogInfo(snes,"SNESSolve_EQ_LS: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 67 68 /* Compute a (scaled) negative update in the line search routine: 69 Y <- X - lambda*Y 70 and evaluate G(Y) = function(Y)) 71 */ 72 ierr = VecCopy(Y,snes->vec_sol_update_always); CHKERRQ(ierr); 73 ierr = (*neP->LineSearch)(snes,X,F,G,Y,W,fnorm,&ynorm,&gnorm,&lsfail); CHKERRQ(ierr); 74 PLogInfo(snes,"SNESSolve_EQ_LS: fnorm=%g, gnorm=%g, ynorm=%g, lsfail=%d\n",fnorm,gnorm,ynorm,lsfail); 75 if (lsfail) snes->nfailures++; 76 77 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 78 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 79 fnorm = gnorm; 80 81 snes->norm = fnorm; 82 if (history && history_len > i+1) history[i+1] = fnorm; 83 SNESMonitor(snes,i+1,fnorm); 84 85 /* Test for convergence */ 86 if (snes->converged) { 87 ierr = VecNorm(X,NORM_2,&xnorm); CHKERRQ(ierr); /* xnorm = || X || */ 88 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 89 break; 90 } 91 } 92 } 93 if (X != snes->vec_sol) { 94 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 95 snes->vec_sol_always = snes->vec_sol; 96 snes->vec_func_always = snes->vec_func; 97 } 98 if (i == maxits) { 99 PLogInfo(snes, 100 "SNESSolve_EQ_LS: Maximum number of iterations has been reached: %d\n",maxits); 101 i--; 102 } 103 if (history) snes->conv_act_size = (history_len < i+1) ? history_len : i+1; 104 *outits = i+1; 105 return 0; 106 } 107 /* ------------------------------------------------------------ */ 108 #undef __FUNC__ 109 #define __FUNC__ "SNESSetUp_EQ_LS" 110 int SNESSetUp_EQ_LS(SNES snes ) 111 { 112 int ierr; 113 snes->nwork = 4; 114 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 115 PLogObjectParents(snes,snes->nwork,snes->work); 116 snes->vec_sol_update_always = snes->work[3]; 117 return 0; 118 } 119 /* ------------------------------------------------------------ */ 120 #undef __FUNC__ 121 #define __FUNC__ "SNESDestroy_EQ_LS" 122 int SNESDestroy_EQ_LS(PetscObject obj) 123 { 124 SNES snes = (SNES) obj; 125 int ierr; 126 if (snes->nwork) { 127 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 128 } 129 PetscFree(snes->data); 130 return 0; 131 } 132 /* ------------------------------------------------------------ */ 133 #undef __FUNC__ 134 #define __FUNC__ "SNESNoLineSearch" 135 /*ARGSUSED*/ 136 /*@C 137 SNESNoLineSearch - This routine is not a line search at all; 138 it simply uses the full Newton step. Thus, this routine is intended 139 to serve as a template and is not recommended for general use. 140 141 Input Parameters: 142 . snes - nonlinear context 143 . x - current iterate 144 . f - residual evaluated at x 145 . y - search direction (contains new iterate on output) 146 . w - work vector 147 . fnorm - 2-norm of f 148 149 Output Parameters: 150 . g - residual evaluated at new iterate y 151 . y - new iterate (contains search direction on input) 152 . gnorm - 2-norm of g 153 . ynorm - 2-norm of search length 154 . flag - set to 0, indicating a successful line search 155 156 Options Database Key: 157 $ -snes_eq_ls basic 158 159 .keywords: SNES, nonlinear, line search, cubic 160 161 .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), 162 SNESSetLineSearch() 163 @*/ 164 int SNESNoLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 165 double fnorm, double *ynorm, double *gnorm,int *flag ) 166 { 167 int ierr; 168 Scalar mone = -1.0; 169 170 *flag = 0; 171 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 172 ierr = VecNorm(y,NORM_2,ynorm); CHKERRQ(ierr); /* ynorm = || y || */ 173 ierr = VecAYPX(&mone,x,y); CHKERRQ(ierr); /* y <- y - x */ 174 ierr = SNESComputeFunction(snes,y,g); CHKERRQ(ierr); /* Compute F(y) */ 175 ierr = VecNorm(g,NORM_2,gnorm); CHKERRQ(ierr); /* gnorm = || g || */ 176 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 177 return 0; 178 } 179 /* ------------------------------------------------------------------ */ 180 #undef __FUNC__ 181 #define __FUNC__ "SNESNoLineSearchNoNorms" 182 /*ARGSUSED*/ 183 /*@C 184 SNESNoLineSearchNoNorms - This routine is not a line search at 185 all; it simply uses the full Newton step. This version does not 186 even compute the norm of the function or search direction; this 187 is intended only when you know the full step is fine and are 188 not checking for convergence of the nonlinear iteration (for 189 example, you are running always for a fixed number of Newton 190 steps). 191 192 Input Parameters: 193 . snes - nonlinear context 194 . x - current iterate 195 . f - residual evaluated at x 196 . y - search direction (contains new iterate on output) 197 . w - work vector 198 . fnorm - 2-norm of f 199 200 Output Parameters: 201 . g - residual evaluated at new iterate y 202 . gnorm - not changed 203 . ynorm - not changed 204 . flag - set to 0, indicating a successful line search 205 206 Options Database Key: 207 $ -snes_eq_ls basicnonorms 208 209 .keywords: SNES, nonlinear, line search, cubic 210 211 .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), 212 SNESSetLineSearch(), SNESNoLineSearch() 213 @*/ 214 int SNESNoLineSearchNoNorms(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 215 double fnorm, double *ynorm, double *gnorm,int *flag ) 216 { 217 int ierr; 218 Scalar mone = -1.0; 219 220 *flag = 0; 221 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 222 ierr = VecAYPX(&mone,x,y); CHKERRQ(ierr); /* y <- y - x */ 223 ierr = SNESComputeFunction(snes,y,g); CHKERRQ(ierr); /* Compute F(y) */ 224 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 225 return 0; 226 } 227 /* ------------------------------------------------------------------ */ 228 #undef __FUNC__ 229 #define __FUNC__ "SNESCubicLineSearch" 230 /*@C 231 SNESCubicLineSearch - Performs a cubic line search (default line search method). 232 233 Input Parameters: 234 . snes - nonlinear context 235 . x - current iterate 236 . f - residual evaluated at x 237 . y - search direction (contains new iterate on output) 238 . w - work vector 239 . fnorm - 2-norm of f 240 241 Output Parameters: 242 . g - residual evaluated at new iterate y 243 . y - new iterate (contains search direction on input) 244 . gnorm - 2-norm of g 245 . ynorm - 2-norm of search length 246 . flag - 0 if line search succeeds; -1 on failure. 247 248 Options Database Key: 249 $ -snes_eq_ls cubic 250 251 Notes: 252 This line search is taken from "Numerical Methods for Unconstrained 253 Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 254 255 .keywords: SNES, nonlinear, line search, cubic 256 257 .seealso: SNESNoLineSearch(), SNESNoLineSearch(), SNESSetLineSearch() 258 @*/ 259 int SNESCubicLineSearch(SNES snes,Vec x,Vec f,Vec g,Vec y,Vec w, 260 double fnorm,double *ynorm,double *gnorm,int *flag) 261 { 262 /* 263 Note that for line search purposes we work with with the related 264 minimization problem: 265 min z(x): R^n -> R, 266 where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2. 267 */ 268 269 double steptol, initslope, lambdaprev, gnormprev, a, b, d, t1, t2; 270 double maxstep, minlambda, alpha, lambda, lambdatemp, lambdaneg; 271 #if defined(PETSC_COMPLEX) 272 Scalar cinitslope, clambda; 273 #endif 274 int ierr, count; 275 SNES_LS *neP = (SNES_LS *) snes->data; 276 Scalar mone = -1.0,scale; 277 278 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 279 *flag = 0; 280 alpha = neP->alpha; 281 maxstep = neP->maxstep; 282 steptol = neP->steptol; 283 284 ierr = VecNorm(y,NORM_2,ynorm); CHKERRQ(ierr); 285 if (*ynorm < snes->atol) { 286 PLogInfo(snes,"SNESCubicLineSearch: Search direction and size are nearly 0\n"); 287 *gnorm = fnorm; 288 ierr = VecCopy(x,y); CHKERRQ(ierr); 289 ierr = VecCopy(f,g); CHKERRQ(ierr); 290 goto theend1; 291 } 292 if (*ynorm > maxstep) { /* Step too big, so scale back */ 293 scale = maxstep/(*ynorm); 294 #if defined(PETSC_COMPLEX) 295 PLogInfo(snes,"SNESCubicLineSearch: Scaling step by %g\n",real(scale)); 296 #else 297 PLogInfo(snes,"SNESCubicLineSearch: Scaling step by %g\n",scale); 298 #endif 299 ierr = VecScale(&scale,y); CHKERRQ(ierr); 300 *ynorm = maxstep; 301 } 302 minlambda = steptol/(*ynorm); 303 ierr = MatMult(snes->jacobian,y,w); CHKERRQ(ierr); 304 #if defined(PETSC_COMPLEX) 305 ierr = VecDot(f,w,&cinitslope); CHKERRQ(ierr); 306 initslope = real(cinitslope); 307 #else 308 ierr = VecDot(f,w,&initslope); CHKERRQ(ierr); 309 #endif 310 if (initslope > 0.0) initslope = -initslope; 311 if (initslope == 0.0) initslope = -1.0; 312 313 ierr = VecCopy(y,w); CHKERRQ(ierr); 314 ierr = VecAYPX(&mone,x,w); CHKERRQ(ierr); 315 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 316 ierr = VecNorm(g,NORM_2,gnorm); 317 if ((*gnorm)*(*gnorm)*0.5 <= fnorm*fnorm*0.5 + alpha*initslope) { /* Sufficient reduction */ 318 ierr = VecCopy(w,y); CHKERRQ(ierr); 319 PLogInfo(snes,"SNESCubicLineSearch: Using full step\n"); 320 goto theend1; 321 } 322 323 /* Fit points with quadratic */ 324 lambda = 1.0; count = 0; 325 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 326 lambdaprev = lambda; 327 gnormprev = *gnorm; 328 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 329 else lambda = lambdatemp; 330 ierr = VecCopy(x,w); CHKERRQ(ierr); 331 lambdaneg = -lambda; 332 #if defined(PETSC_COMPLEX) 333 clambda = lambdaneg; ierr = VecAXPY(&clambda,y,w); CHKERRQ(ierr); 334 #else 335 ierr = VecAXPY(&lambdaneg,y,w); CHKERRQ(ierr); 336 #endif 337 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 338 ierr = VecNorm(g,NORM_2,gnorm); CHKERRQ(ierr); 339 if ((*gnorm)*(*gnorm)*0.5 <= fnorm*fnorm*0.5 + alpha*initslope) { /* sufficient reduction */ 340 ierr = VecCopy(w,y); CHKERRQ(ierr); 341 PLogInfo(snes,"SNESCubicLineSearch: Quadratically determined step, lambda=%g\n",lambda); 342 goto theend1; 343 } 344 345 /* Fit points with cubic */ 346 count = 1; 347 while (1) { 348 if (lambda <= minlambda) { /* bad luck; use full step */ 349 PLogInfo(snes,"SNESCubicLineSearch:Unable to find good step length! %d \n",count); 350 PLogInfo(snes,"SNESCubicLineSearch:fnorm=%g, gnorm=%g, ynorm=%g, lambda=%g, initial slope=%g\n", 351 fnorm,*gnorm,*ynorm,lambda,initslope); 352 ierr = VecCopy(w,y); CHKERRQ(ierr); 353 *flag = -1; break; 354 } 355 t1 = ((*gnorm)*(*gnorm) - fnorm*fnorm)*0.5 - lambda*initslope; 356 t2 = (gnormprev*gnormprev - fnorm*fnorm)*0.5 - lambdaprev*initslope; 357 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 358 b = (-lambdaprev*t1/(lambda*lambda) + lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 359 d = b*b - 3*a*initslope; 360 if (d < 0.0) d = 0.0; 361 if (a == 0.0) { 362 lambdatemp = -initslope/(2.0*b); 363 } else { 364 lambdatemp = (-b + sqrt(d))/(3.0*a); 365 } 366 if (lambdatemp > .5*lambda) { 367 lambdatemp = .5*lambda; 368 } 369 lambdaprev = lambda; 370 gnormprev = *gnorm; 371 if (lambdatemp <= .1*lambda) { 372 lambda = .1*lambda; 373 } 374 else lambda = lambdatemp; 375 ierr = VecCopy( x, w ); CHKERRQ(ierr); 376 lambdaneg = -lambda; 377 #if defined(PETSC_COMPLEX) 378 clambda = lambdaneg; 379 ierr = VecAXPY(&clambda,y,w); CHKERRQ(ierr); 380 #else 381 ierr = VecAXPY(&lambdaneg,y,w); CHKERRQ(ierr); 382 #endif 383 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 384 ierr = VecNorm(g,NORM_2,gnorm); CHKERRQ(ierr); 385 if ((*gnorm)*(*gnorm)*0.5 <= fnorm*fnorm*0.5 + alpha*initslope) { /* is reduction enough? */ 386 ierr = VecCopy(w,y); CHKERRQ(ierr); 387 PLogInfo(snes,"SNESCubicLineSearch: Cubically determined step, lambda=%g\n",lambda); 388 goto theend1; 389 } else { 390 PLogInfo(snes,"SNESCubicLineSearch: Cubic step no good, shrinking lambda, lambda=%g\n",lambda); 391 } 392 count++; 393 } 394 theend1: 395 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 396 return 0; 397 } 398 /* ------------------------------------------------------------------ */ 399 #undef __FUNC__ 400 #define __FUNC__ "SNESQuadraticLineSearch" 401 /*@C 402 SNESQuadraticLineSearch - Performs a quadratic line search. 403 404 Input Parameters: 405 . snes - the SNES context 406 . x - current iterate 407 . f - residual evaluated at x 408 . y - search direction (contains new iterate on output) 409 . w - work vector 410 . fnorm - 2-norm of f 411 412 Output Parameters: 413 . g - residual evaluated at new iterate y 414 . y - new iterate (contains search direction on input) 415 . gnorm - 2-norm of g 416 . ynorm - 2-norm of search length 417 . flag - 0 if line search succeeds; -1 on failure. 418 419 Options Database Key: 420 $ -snes_eq_ls quadratic 421 422 Notes: 423 Use SNESSetLineSearch() 424 to set this routine within the SNES_EQ_LS method. 425 426 .keywords: SNES, nonlinear, quadratic, line search 427 428 .seealso: SNESCubicLineSearch(), SNESNoLineSearch(), SNESSetLineSearch() 429 @*/ 430 int SNESQuadraticLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 431 double fnorm, double *ynorm, double *gnorm,int *flag) 432 { 433 /* 434 Note that for line search purposes we work with with the related 435 minimization problem: 436 min z(x): R^n -> R, 437 where z(x) = .5 * fnorm*fnorm, and fnorm = || f ||_2. 438 */ 439 double steptol,initslope,lambdaprev,gnormprev,maxstep,minlambda,alpha,lambda,lambdatemp; 440 #if defined(PETSC_COMPLEX) 441 Scalar cinitslope,clambda; 442 #endif 443 int ierr,count; 444 SNES_LS *neP = (SNES_LS *) snes->data; 445 Scalar mone = -1.0,scale; 446 447 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 448 *flag = 0; 449 alpha = neP->alpha; 450 maxstep = neP->maxstep; 451 steptol = neP->steptol; 452 453 VecNorm(y, NORM_2,ynorm ); 454 if (*ynorm < snes->atol) { 455 PLogInfo(snes,"SNESQuadraticLineSearch: Search direction and size is 0\n"); 456 *gnorm = fnorm; 457 ierr = VecCopy(x,y); CHKERRQ(ierr); 458 ierr = VecCopy(f,g); CHKERRQ(ierr); 459 goto theend2; 460 } 461 if (*ynorm > maxstep) { /* Step too big, so scale back */ 462 scale = maxstep/(*ynorm); 463 ierr = VecScale(&scale,y); CHKERRQ(ierr); 464 *ynorm = maxstep; 465 } 466 minlambda = steptol/(*ynorm); 467 ierr = MatMult(snes->jacobian,y,w); CHKERRQ(ierr); 468 #if defined(PETSC_COMPLEX) 469 ierr = VecDot(f,w,&cinitslope); CHKERRQ(ierr); 470 initslope = real(cinitslope); 471 #else 472 ierr = VecDot(f,w,&initslope); CHKERRQ(ierr); 473 #endif 474 if (initslope > 0.0) initslope = -initslope; 475 if (initslope == 0.0) initslope = -1.0; 476 477 ierr = VecCopy(y,w); CHKERRQ(ierr); 478 ierr = VecAYPX(&mone,x,w); CHKERRQ(ierr); 479 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 480 ierr = VecNorm(g,NORM_2,gnorm); CHKERRQ(ierr); 481 if ((*gnorm)*(*gnorm)*0.5 <= fnorm*fnorm*0.5 + alpha*initslope) { /* Sufficient reduction */ 482 ierr = VecCopy(w,y); CHKERRQ(ierr); 483 PLogInfo(snes,"SNESQuadraticLineSearch: Using full step\n"); 484 goto theend2; 485 } 486 487 /* Fit points with quadratic */ 488 lambda = 1.0; count = 0; 489 count = 1; 490 while (1) { 491 if (lambda <= minlambda) { /* bad luck; use full step */ 492 PLogInfo(snes, 493 "SNESQuadraticLineSearch:Unable to find good step length! %d \n",count); 494 PLogInfo(snes, 495 "SNESQuadraticLineSearch:fnorm=%g, gnorm=%g, ynorm=%g, lambda=%g, initial slope=%g\n", 496 fnorm,*gnorm,*ynorm,lambda,initslope); 497 ierr = VecCopy(w,y); CHKERRQ(ierr); 498 *flag = -1; break; 499 } 500 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 501 lambdaprev = lambda; 502 gnormprev = *gnorm; 503 if (lambdatemp <= .1*lambda) { 504 lambda = .1*lambda; 505 } else lambda = lambdatemp; 506 ierr = VecCopy(x,w); CHKERRQ(ierr); 507 lambda = -lambda; 508 #if defined(PETSC_COMPLEX) 509 clambda = lambda; ierr = VecAXPY(&clambda,y,w); CHKERRQ(ierr); 510 #else 511 ierr = VecAXPY(&lambda,y,w); CHKERRQ(ierr); 512 #endif 513 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 514 ierr = VecNorm(g,NORM_2,gnorm); CHKERRQ(ierr); 515 if ((*gnorm)*(*gnorm)*0.5 <= fnorm*fnorm*0.5 + alpha*initslope) { /* sufficient reduction */ 516 ierr = VecCopy(w,y); CHKERRQ(ierr); 517 PLogInfo(snes, 518 "SNESQuadraticLineSearch:Quadratically determined step, lambda=%g\n",lambda); 519 break; 520 } 521 count++; 522 } 523 theend2: 524 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 525 return 0; 526 } 527 /* ------------------------------------------------------------ */ 528 #undef __FUNC__ 529 #define __FUNC__ "SNESSetLineSearch" 530 /*@C 531 SNESSetLineSearch - Sets the line search routine to be used 532 by the method SNES_EQ_LS. 533 534 Input Parameters: 535 . snes - nonlinear context obtained from SNESCreate() 536 . func - pointer to int function 537 538 Available Routines: 539 . SNESCubicLineSearch() - default line search 540 . SNESQuadraticLineSearch() - quadratic line search 541 . SNESNoLineSearch() - the full Newton step (actually not a line search) 542 543 Options Database Keys: 544 $ -snes_eq_ls [basic,quadratic,cubic] 545 $ -snes_eq_ls_alpha <alpha> 546 $ -snes_eq_ls_maxstep <max> 547 $ -snes_eq_ls_steptol <steptol> 548 549 Calling sequence of func: 550 func (SNES snes, Vec x, Vec f, Vec g, Vec y, 551 Vec w, double fnorm, double *ynorm, 552 double *gnorm, *flag) 553 554 Input parameters for func: 555 . snes - nonlinear context 556 . x - current iterate 557 . f - residual evaluated at x 558 . y - search direction (contains new iterate on output) 559 . w - work vector 560 . fnorm - 2-norm of f 561 562 Output parameters for func: 563 . g - residual evaluated at new iterate y 564 . y - new iterate (contains search direction on input) 565 . gnorm - 2-norm of g 566 . ynorm - 2-norm of search length 567 . flag - set to 0 if the line search succeeds; a nonzero integer 568 on failure. 569 570 .keywords: SNES, nonlinear, set, line search, routine 571 572 .seealso: SNESNoLineSearch(), SNESQuadraticLineSearch(), SNESCubicLineSearch() 573 @*/ 574 int SNESSetLineSearch(SNES snes,int (*func)(SNES,Vec,Vec,Vec,Vec,Vec, 575 double,double*,double*,int*)) 576 { 577 if ((snes)->type == SNES_EQ_LS) ((SNES_LS *)(snes->data))->LineSearch = func; 578 return 0; 579 } 580 /* ------------------------------------------------------------------ */ 581 #undef __FUNC__ 582 #define __FUNC__ "SNESPrintHelp_EQ_LS" 583 static int SNESPrintHelp_EQ_LS(SNES snes,char *p) 584 { 585 SNES_LS *ls = (SNES_LS *)snes->data; 586 587 PetscPrintf(snes->comm," method SNES_EQ_LS (ls) for systems of nonlinear equations:\n"); 588 PetscPrintf(snes->comm," %ssnes_eq_ls [basic,quadratic,cubic]\n",p); 589 PetscPrintf(snes->comm," %ssnes_eq_ls_alpha <alpha> (default %g)\n",p,ls->alpha); 590 PetscPrintf(snes->comm," %ssnes_eq_ls_maxstep <max> (default %g)\n",p,ls->maxstep); 591 PetscPrintf(snes->comm," %ssnes_eq_ls_steptol <tol> (default %g)\n",p,ls->steptol); 592 return 0; 593 } 594 /* ------------------------------------------------------------------ */ 595 #undef __FUNC__ 596 #define __FUNC__ "SNESView_EQ_LS" 597 static int SNESView_EQ_LS(PetscObject obj,Viewer viewer) 598 { 599 SNES snes = (SNES)obj; 600 SNES_LS *ls = (SNES_LS *)snes->data; 601 FILE *fd; 602 char *cstr; 603 int ierr; 604 ViewerType vtype; 605 606 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 607 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 608 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 609 if (ls->LineSearch == SNESNoLineSearch) cstr = "SNESNoLineSearch"; 610 else if (ls->LineSearch == SNESQuadraticLineSearch) cstr = "SNESQuadraticLineSearch"; 611 else if (ls->LineSearch == SNESCubicLineSearch) cstr = "SNESCubicLineSearch"; 612 else cstr = "unknown"; 613 PetscFPrintf(snes->comm,fd," line search variant: %s\n",cstr); 614 PetscFPrintf(snes->comm,fd," alpha=%g, maxstep=%g, steptol=%g\n", 615 ls->alpha,ls->maxstep,ls->steptol); 616 } 617 return 0; 618 } 619 /* ------------------------------------------------------------------ */ 620 #undef __FUNC__ 621 #define __FUNC__ "SNESSetFromOptions_EQ_LS" 622 static int SNESSetFromOptions_EQ_LS(SNES snes) 623 { 624 SNES_LS *ls = (SNES_LS *)snes->data; 625 char ver[16]; 626 double tmp; 627 int ierr,flg; 628 629 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_ls_alpha",&tmp, &flg);CHKERRQ(ierr); 630 if (flg) { 631 ls->alpha = tmp; 632 } 633 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_ls_maxstep",&tmp, &flg);CHKERRQ(ierr); 634 if (flg) { 635 ls->maxstep = tmp; 636 } 637 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_ls_steptol",&tmp, &flg);CHKERRQ(ierr); 638 if (flg) { 639 ls->steptol = tmp; 640 } 641 ierr = OptionsGetString(snes->prefix,"-snes_eq_ls",ver,16, &flg); CHKERRQ(ierr); 642 if (flg) { 643 if (!PetscStrcmp(ver,"basic")) { 644 SNESSetLineSearch(snes,SNESNoLineSearch); 645 } 646 else if (!PetscStrcmp(ver,"basicnonorms")) { 647 SNESSetLineSearch(snes,SNESNoLineSearchNoNorms); 648 } 649 else if (!PetscStrcmp(ver,"quadratic")) { 650 SNESSetLineSearch(snes,SNESQuadraticLineSearch); 651 } 652 else if (!PetscStrcmp(ver,"cubic")) { 653 SNESSetLineSearch(snes,SNESCubicLineSearch); 654 } 655 else {SETERRQ(1,0,"Unknown line search");} 656 } 657 return 0; 658 } 659 /* ------------------------------------------------------------ */ 660 #undef __FUNC__ 661 #define __FUNC__ "SNESCreate_EQ_LS" 662 int SNESCreate_EQ_LS(SNES snes ) 663 { 664 SNES_LS *neP; 665 666 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) SETERRQ(1,0,"For SNES_NONLINEAR_EQUATIONS only"); 667 snes->type = SNES_EQ_LS; 668 snes->setup = SNESSetUp_EQ_LS; 669 snes->solve = SNESSolve_EQ_LS; 670 snes->destroy = SNESDestroy_EQ_LS; 671 snes->converged = SNESConverged_EQ_LS; 672 snes->printhelp = SNESPrintHelp_EQ_LS; 673 snes->setfromoptions = SNESSetFromOptions_EQ_LS; 674 snes->view = SNESView_EQ_LS; 675 snes->nwork = 0; 676 677 neP = PetscNew(SNES_LS); CHKPTRQ(neP); 678 PLogObjectMemory(snes,sizeof(SNES_LS)); 679 snes->data = (void *) neP; 680 neP->alpha = 1.e-4; 681 neP->maxstep = 1.e8; 682 neP->steptol = 1.e-12; 683 neP->LineSearch = SNESCubicLineSearch; 684 return 0; 685 } 686 687