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