1 #ifndef lint 2 static char vcid[] = "$Id: ls.c,v 1.7 1995/04/17 03:35:40 curfman Exp curfman $"; 3 #endif 4 5 #include <math.h> 6 #include "ls.h" 7 8 /* 9 Implements a line search variant of Newton's Method 10 for solving systems of nonlinear equations. 11 12 Input parameters: 13 . snes - nonlinear context obtained from SNESCreate() 14 15 Output Parameters: 16 . its - Number of global iterations until termination. 17 18 Notes: 19 See SNESCreate() and SNESSetUp() for information on the definition and 20 initialization of the nonlinear solver context. 21 22 This implements essentially a truncated Newton method with a 23 line search. By default a cubic backtracking line search 24 is employed, as described in the text "Numerical Methods for 25 Unconstrained Optimization and Nonlinear Equations" by Dennis 26 and Schnabel. See the examples in src/snes/examples. 27 */ 28 29 int SNESSolve_LS( SNES snes, int *outits ) 30 { 31 SNES_LS *neP = (SNES_LS *) snes->data; 32 int maxits, i, history_len,ierr,lits; 33 double fnorm, gnorm, xnorm, ynorm, *history; 34 Vec Y, X, F, G, W, TMP; 35 36 history = snes->conv_hist; /* convergence history */ 37 history_len = snes->conv_hist_len; /* convergence history length */ 38 maxits = snes->max_its; /* maximum number of iterations */ 39 X = snes->vec_sol; /* solution vector */ 40 F = snes->vec_res; /* residual vector */ 41 Y = snes->work[0]; /* work vectors */ 42 G = snes->work[1]; 43 W = snes->work[2]; 44 45 ierr = SNESComputeInitialGuess(snes,X); CHKERR(ierr); /* X <- X_0 */ 46 VecNorm( X, &xnorm ); /* xnorm = || X || */ 47 ierr = SNESComputeResidual(snes,X,F); CHKERR(ierr); /* (+/-) F(X) */ 48 VecNorm(F, &fnorm ); /* fnorm <- || F || */ 49 snes->norm = fnorm; 50 if (history && history_len > 0) history[0] = fnorm; 51 if (snes->Monitor) (*snes->Monitor)(snes,0,X,F,fnorm,snes->monP); 52 53 for ( i=0; i<maxits; i++ ) { 54 snes->iter = i+1; 55 56 /* Solve J Y = -F, where J is Jacobian matrix */ 57 (*snes->ComputeJacobian)(X,&snes->jacobian,snes->jacP); 58 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian,0); 59 ierr = SLESSolve(snes->sles,F,Y,&lits); CHKERR(ierr); 60 ierr = (*neP->LineSearch)(snes, X, F, G, Y, W, fnorm, &ynorm, &gnorm ); 61 62 TMP = F; F = G; G = TMP; 63 TMP = X; X = Y; Y = TMP; 64 fnorm = gnorm; 65 66 snes->norm = fnorm; 67 if (history && history_len > i+1) history[i+1] = fnorm; 68 VecNorm( X, &xnorm ); /* xnorm = || X || */ 69 if (snes->Monitor) (*snes->Monitor)(snes,i+1,X,F,fnorm,snes->monP); 70 71 /* Test for convergence */ 72 if ((*snes->Converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 73 if (X != snes->vec_sol) VecCopy( X, snes->vec_sol ); 74 break; 75 } 76 } 77 if (i == maxits) i--; 78 *outits = i+1; 79 return 0; 80 } 81 /* ------------------------------------------------------------ */ 82 /*ARGSUSED*/ 83 int SNESSetUp_LS(SNES snes ) 84 { 85 int ierr; 86 snes->nwork = 3; 87 ierr = VecGetVecs( snes->vec_sol, snes->nwork,&snes->work ); CHKERR(ierr); 88 PLogObjectParents(snes,snes->nwork,snes->work ); 89 return 0; 90 } 91 /* ------------------------------------------------------------ */ 92 /*ARGSUSED*/ 93 int SNESDestroy_LS(PetscObject obj) 94 { 95 SNES snes = (SNES) obj; 96 SLESDestroy(snes->sles); 97 VecFreeVecs(snes->work, snes->nwork ); 98 PLogObjectDestroy(obj); 99 PETSCHEADERDESTROY(obj); 100 return 0; 101 } 102 /*@ 103 SNESDefaultMonitor - Default SNES monitoring routine. Prints the 104 residual norm at each iteration. 105 106 Input Parameters: 107 . snes - the SNES context 108 . its - iteration number 109 . x - current iterate 110 . f - current residual (+/-) 111 . fnorm - 2-norm residual value (may be estimated) 112 . dummy - unused context 113 114 Notes: 115 f is either the residual or its negative, depending on the user's 116 preference, as set with SNESSetResidual(). 117 118 .keywords: SNES, nonlinear, default, monitor, residual, norm 119 120 .seealso: SNESSetMonitor() 121 @*/ 122 int SNESDefaultMonitor(SNES snes,int its, Vec x,Vec f,double fnorm,void *dummy) 123 { 124 fprintf( stdout, "iter = %d, residual norm %g \n",its,fnorm); 125 return 0; 126 } 127 128 /*@ 129 SNESDefaultConverged - Default test for monitoring the convergence 130 of the SNES solvers. 131 132 Input Parameters: 133 . snes - the SNES context 134 . xnorm - 2-norm of current iterate 135 . pnorm - 2-norm of current step 136 . fnorm - 2-norm of residual 137 . dummy - unused context 138 139 Returns: 140 $ 2 if ( fnorm < atol ), 141 $ 3 if ( pnorm < xtol*xnorm ), 142 $ -2 if ( nres > max_res ), 143 $ 0 otherwise, 144 145 where 146 $ atol - absolute residual norm tolerance, 147 $ set with SNESSetAbsoluteTolerance() 148 $ max_res - maximum number of residual evaluations, 149 $ set with SNESSetMaxResidualEvaluations() 150 $ nres - number of residual evaluations 151 $ xtol - relative residual norm tolerance, 152 $ set with SNESSetRelativeTolerance() 153 154 .keywords: SNES, nonlinear, default, converged, convergence 155 156 .seealso: SNESSetConvergenceTest() 157 @*/ 158 int SNESDefaultConverged(SNES snes,double xnorm,double pnorm,double fnorm, 159 void *dummy) 160 { 161 if (fnorm < snes->atol) { 162 PLogInfo((PetscObject)snes, 163 "Converged due to absolute residual norm %g < %g\n",fnorm,snes->atol); 164 return 2; 165 } 166 if (pnorm < snes->xtol*(xnorm)) { 167 PLogInfo((PetscObject)snes, 168 "Converged due to small update length %g < %g*%g\n", 169 pnorm,snes->xtol,xnorm); 170 return 3; 171 } 172 if (snes->nresids > snes->max_resids) { 173 PLogInfo((PetscObject)snes, 174 "Exceeded maximum number of residual evaluations: %d > %d\n", 175 snes->nresids, snes->max_resids ); 176 return -2; 177 } 178 return 0; 179 } 180 181 /* ------------------------------------------------------------ */ 182 /*ARGSUSED*/ 183 /*@ 184 SNESNoLineSearch - This routine is not a line search at all; 185 it simply uses the full Newton step. Thus, this routine is intended 186 to serve as a template and is not recommended for general use. 187 188 Input Parameters: 189 . snes - nonlinear context 190 . x - current iterate 191 . f - residual evaluated at x 192 . y - search direction (contains new iterate on output) 193 . w - work vector 194 . fnorm - 2-norm of f 195 196 Output parameters: 197 . g - residual evaluated at new iterate y 198 . y - new iterate (contains search direction on input) 199 . gnorm - 2-norm of g 200 . ynorm - 2-norm of search length 201 202 Returns: 203 1, indicating success of the step. 204 205 .keywords: SNES, nonlinear, line search, cubic 206 207 .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), 208 .seealso: SNESSetLineSearchRoutine() 209 @*/ 210 int SNESNoLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 211 double fnorm, double *ynorm, double *gnorm ) 212 { 213 int ierr; 214 Scalar one = 1.0; 215 VecNorm(y, ynorm ); /* ynorm = || y || */ 216 VecAXPY(&one, x, y ); /* y <- x + y */ 217 ierr = SNESComputeResidual(snes,y,g); CHKERR(ierr); 218 VecNorm( g, gnorm ); /* gnorm = || g || */ 219 return 1; 220 } 221 /* ------------------------------------------------------------------ */ 222 /*@ 223 SNESCubicLineSearch - This routine performs a cubic line search and 224 is the default line search method. 225 226 Input Parameters: 227 . snes - nonlinear context 228 . x - current iterate 229 . f - residual evaluated at x 230 . y - search direction (contains new iterate on output) 231 . w - work vector 232 . fnorm - 2-norm of f 233 234 Output parameters: 235 . g - residual evaluated at new iterate y 236 . y - new iterate (contains search direction on input) 237 . gnorm - 2-norm of g 238 . ynorm - 2-norm of search length 239 240 Returns: 241 1 if the line search succeeds; 0 if the line search fails. 242 243 Notes: 244 This line search is taken from "Numerical Methods for Unconstrained 245 Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 246 247 .keywords: SNES, nonlinear, line search, cubic 248 249 .seealso: SNESNoLineSearch(), SNESNoLineSearch(), SNESSetLineSearchRoutine() 250 @*/ 251 int SNESCubicLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 252 double fnorm, double *ynorm, double *gnorm ) 253 { 254 double steptol, initslope; 255 double lambdaprev, gnormprev; 256 double a, b, d, t1, t2; 257 #if defined(PETSC_COMPLEX) 258 Scalar cinitslope,clambda; 259 #endif 260 int ierr,count; 261 SNES_LS *neP = (SNES_LS *) snes->data; 262 Scalar one = 1.0,scale; 263 double maxstep,minlambda,alpha,lambda,lambdatemp; 264 265 alpha = neP->alpha; 266 maxstep = neP->maxstep; 267 steptol = neP->steptol; 268 269 VecNorm(y, ynorm ); 270 if (*ynorm > maxstep) { /* Step too big, so scale back */ 271 scale = maxstep/(*ynorm); 272 #if defined(PETSC_COMPLEX) 273 PLogInfo((PetscObject)snes,"Scaling step by %g\n",real(scale)); 274 #else 275 PLogInfo((PetscObject)snes,"Scaling step by %g\n",scale); 276 #endif 277 VecScale(&scale, y ); 278 *ynorm = maxstep; 279 } 280 minlambda = steptol/(*ynorm); 281 #if defined(PETSC_COMPLEX) 282 VecDot(f, y, &cinitslope ); 283 initslope = real(cinitslope); 284 #else 285 VecDot(f, y, &initslope ); 286 #endif 287 if (initslope > 0.0) initslope = -initslope; 288 if (initslope == 0.0) initslope = -1.0; 289 290 VecCopy(y, w ); 291 VecAXPY(&one, x, w ); 292 ierr = SNESComputeResidual(snes,w,g); CHKERR(ierr); 293 VecNorm(g, gnorm ); 294 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 295 VecCopy(w, y ); 296 PLogInfo((PetscObject)snes,"Using full step\n"); 297 return 0; 298 } 299 300 /* Fit points with quadratic */ 301 lambda = 1.0; count = 0; 302 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 303 lambdaprev = lambda; 304 gnormprev = *gnorm; 305 if (lambdatemp <= .1*lambda) { 306 lambda = .1*lambda; 307 } else lambda = lambdatemp; 308 VecCopy(x, w ); 309 #if defined(PETSC_COMPLEX) 310 clambda = lambda; VecAXPY(&clambda, y, w ); 311 #else 312 VecAXPY(&lambda, y, w ); 313 #endif 314 ierr = SNESComputeResidual(snes,w,g); CHKERR(ierr); 315 VecNorm(g, gnorm ); 316 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 317 VecCopy(w, y ); 318 PLogInfo((PetscObject)snes,"Quadratically determined step, lambda %g\n",lambda); 319 return 0; 320 } 321 322 /* Fit points with cubic */ 323 count = 1; 324 while (1) { 325 if (lambda <= minlambda) { /* bad luck; use full step */ 326 PLogInfo((PetscObject)snes,"Unable to find good step length! %d \n",count); 327 PLogInfo((PetscObject)snes, "f %g fnew %g ynorm %g lambda %g \n", 328 fnorm,*gnorm, *ynorm,lambda); 329 VecCopy(w, y ); 330 return 0; 331 } 332 t1 = *gnorm - fnorm - lambda*initslope; 333 t2 = gnormprev - fnorm - lambdaprev*initslope; 334 a = (t1/(lambda*lambda) - 335 t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 336 b = (-lambdaprev*t1/(lambda*lambda) + 337 lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 338 d = b*b - 3*a*initslope; 339 if (d < 0.0) d = 0.0; 340 if (a == 0.0) { 341 lambdatemp = -initslope/(2.0*b); 342 } else { 343 lambdatemp = (-b + sqrt(d))/(3.0*a); 344 } 345 if (lambdatemp > .5*lambda) { 346 lambdatemp = .5*lambda; 347 } 348 lambdaprev = lambda; 349 gnormprev = *gnorm; 350 if (lambdatemp <= .1*lambda) { 351 lambda = .1*lambda; 352 } 353 else lambda = lambdatemp; 354 VecCopy( x, w ); 355 #if defined(PETSC_COMPLEX) 356 VecAXPY(&clambda, y, w ); 357 #else 358 VecAXPY(&lambda, y, w ); 359 #endif 360 ierr = SNESComputeResidual(snes,w,g); CHKERR(ierr); 361 VecNorm(g, gnorm ); 362 if (*gnorm <= fnorm + alpha*initslope) { /* is reduction enough */ 363 VecCopy(w, y ); 364 PLogInfo((PetscObject)snes,"Cubically determined step, lambda %g\n",lambda); 365 return 0; 366 } 367 count++; 368 } 369 return 0; 370 } 371 /*@ 372 SNESQuadraticLineSearch - This routine performs a cubic line search. 373 374 Input Parameters: 375 . snes - the SNES context 376 . x - current iterate 377 . f - residual evaluated at x 378 . y - search direction (contains new iterate on output) 379 . w - work vector 380 . fnorm - 2-norm of f 381 382 Output parameters: 383 . g - residual evaluated at new iterate y 384 . y - new iterate (contains search direction on input) 385 . gnorm - 2-norm of g 386 . ynorm - 2-norm of search length 387 388 Returns: 389 1 if the line search succeeds; 0 if the line search fails. 390 391 Notes: 392 Use SNESSetLineSearchRoutine() 393 to set this routine within the SNES_NLS method. 394 395 .keywords: SNES, nonlinear, quadratic, line search 396 397 .seealso: SNESCubicLineSearch(), SNESNoLineSearch(), SNESSetLineSearchRoutine() 398 @*/ 399 int SNESQuadraticLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 400 double fnorm, double *ynorm, double *gnorm ) 401 { 402 double steptol, initslope; 403 double lambdaprev, gnormprev; 404 #if defined(PETSC_COMPLEX) 405 Scalar cinitslope,clambda; 406 #endif 407 int ierr,count; 408 SNES_LS *neP = (SNES_LS *) snes->data; 409 Scalar one = 1.0,scale; 410 double maxstep,minlambda,alpha,lambda,lambdatemp; 411 412 alpha = neP->alpha; 413 maxstep = neP->maxstep; 414 steptol = neP->steptol; 415 416 VecNorm(y, ynorm ); 417 if (*ynorm > maxstep) { /* Step too big, so scale back */ 418 scale = maxstep/(*ynorm); 419 VecScale(&scale, y ); 420 *ynorm = maxstep; 421 } 422 minlambda = steptol/(*ynorm); 423 #if defined(PETSC_COMPLEX) 424 VecDot(f, y, &cinitslope ); 425 initslope = real(cinitslope); 426 #else 427 VecDot(f, y, &initslope ); 428 #endif 429 if (initslope > 0.0) initslope = -initslope; 430 if (initslope == 0.0) initslope = -1.0; 431 432 VecCopy(y, w ); 433 VecAXPY(&one, x, w ); 434 ierr = SNESComputeResidual(snes,w,g); CHKERR(ierr); 435 VecNorm(g, gnorm ); 436 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 437 VecCopy(w, y ); 438 PLogInfo((PetscObject)snes,"Using full step\n"); 439 return 0; 440 } 441 442 /* Fit points with quadratic */ 443 lambda = 1.0; count = 0; 444 count = 1; 445 while (1) { 446 if (lambda <= minlambda) { /* bad luck; use full step */ 447 PLogInfo((PetscObject)snes,"Unable to find good step length! %d \n",count); 448 PLogInfo((PetscObject)snes, "f %g fnew %g ynorm %g lambda %g \n", 449 fnorm,*gnorm, *ynorm,lambda); 450 VecCopy(w, y ); 451 return 0; 452 } 453 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 454 lambdaprev = lambda; 455 gnormprev = *gnorm; 456 if (lambdatemp <= .1*lambda) { 457 lambda = .1*lambda; 458 } else lambda = lambdatemp; 459 VecCopy(x, w ); 460 #if defined(PETSC_COMPLEX) 461 clambda = lambda; VecAXPY(&clambda, y, w ); 462 #else 463 VecAXPY(&lambda, y, w ); 464 #endif 465 ierr = SNESComputeResidual(snes,w,g); CHKERR(ierr); 466 VecNorm(g, gnorm ); 467 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 468 VecCopy(w, y ); 469 PLogInfo((PetscObject)snes,"Quadratically determined step, lambda %g\n",lambda); 470 return 0; 471 } 472 count++; 473 } 474 475 return 0; 476 } 477 /* ------------------------------------------------------------ */ 478 /*@C 479 SNESSetLineSearchRoutine - Sets the line search routine to be used 480 by the method SNES_LS. 481 482 Input Parameters: 483 . snes - nonlinear context obtained from NLCreate() 484 . func - pointer to int function 485 486 Possible routines: 487 . SNESCubicLineSearch() - default line search 488 . SNESQuadraticLineSearch() - quadratic line search 489 . SNESNoLineSearch() - the full Newton step (actually not a line search) 490 491 Calling sequence of func: 492 func (SNES snes, Vec x, Vec f, Vec g, Vec y, 493 Vec w, double fnorm, double *ynorm, double *gnorm) 494 495 Input parameters for func: 496 . snes - nonlinear context 497 . x - current iterate 498 . f - residual evaluated at x 499 . y - search direction (contains new iterate on output) 500 . w - work vector 501 . fnorm - 2-norm of f 502 503 Output parameters for func: 504 . g - residual evaluated at new iterate y 505 . y - new iterate (contains search direction on input) 506 . gnorm - 2-norm of g 507 . ynorm - 2-norm of search length 508 509 Returned by func: 510 1 if the line search succeeds; 0 if the line search fails. 511 512 .keywords: SNES, nonlinear, set, line search, routine 513 514 .seealso: SNESNoLineSearch(), SNESQuadraticLineSearch(), SNESCubicLineSearch() 515 @*/ 516 int SNESSetLineSearchRoutine(SNES snes,int (*func)(SNES,Vec,Vec,Vec,Vec,Vec, 517 double,double *,double*) ) 518 { 519 if ((snes)->type == SNES_NLS) 520 ((SNES_LS *)(snes->data))->LineSearch = func; 521 return 0; 522 } 523 524 static int SNESPrintHelp_LS(SNES snes) 525 { 526 SNES_LS *ls = (SNES_LS *)snes->data; 527 fprintf(stderr,"-snes_line_search [basic,quadratic,cubic]\n"); 528 fprintf(stderr,"-snes_line_search_alpha alpha (default %g)\n",ls->alpha); 529 fprintf(stderr,"-snes_line_search_maxstep max (default %g)\n",ls->maxstep); 530 fprintf(stderr,"-snes_line_search_steptol tol (default %g)\n",ls->steptol); 531 return 0; 532 } 533 534 #include "options.h" 535 static int SNESSetFromOptions_LS(SNES snes) 536 { 537 SNES_LS *ls = (SNES_LS *)snes->data; 538 char ver[16]; 539 double tmp; 540 541 if (OptionsGetDouble(0,snes->prefix,"-snes_line_search_alpa",&tmp)) { 542 ls->alpha = tmp; 543 } 544 if (OptionsGetDouble(0,snes->prefix,"-snes_line_search_maxstep",&tmp)) { 545 ls->maxstep = tmp; 546 } 547 if (OptionsGetDouble(0,snes->prefix,"-snes_line_search_steptol",&tmp)) { 548 ls->steptol = tmp; 549 } 550 if (OptionsGetString(0,snes->prefix,"-snes_line_search",ver,16)) { 551 if (!strcmp(ver,"basic")) { 552 SNESSetLineSearchRoutine(snes,SNESNoLineSearch); 553 } 554 else if (!strcmp(ver,"quadratic")) { 555 SNESSetLineSearchRoutine(snes,SNESQuadraticLineSearch); 556 } 557 else if (!strcmp(ver,"cubic")) { 558 SNESSetLineSearchRoutine(snes,SNESCubicLineSearch); 559 } 560 else {SETERR(1,"Unknown line search?");} 561 } 562 return 0; 563 } 564 565 /* ------------------------------------------------------------ */ 566 int SNESCreate_LS(SNES snes ) 567 { 568 SNES_LS *neP; 569 570 snes->type = SNES_NLS; 571 snes->Setup = SNESSetUp_LS; 572 snes->Solver = SNESSolve_LS; 573 snes->destroy = SNESDestroy_LS; 574 snes->Monitor = 0; 575 snes->Converged = SNESDefaultConverged; 576 snes->PrintHelp = SNESPrintHelp_LS; 577 snes->SetFromOptions = SNESSetFromOptions_LS; 578 579 neP = NEW(SNES_LS); CHKPTR(neP); 580 snes->data = (void *) neP; 581 neP->alpha = 1.e-4; 582 neP->maxstep = 1.e8; 583 neP->steptol = 1.e-12; 584 neP->LineSearch = SNESCubicLineSearch; 585 return 0; 586 } 587 588 589 590 591