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