1 #ifndef lint 2 static char vcid[] = "$Id: ls.c,v 1.10 1995/04/19 03:01:24 bsmith Exp bsmith $"; 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 = SNESComputeFunction(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,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,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 . fnorm - 2-norm residual value (may be estimated) 110 . dummy - unused context 111 112 Notes: 113 f is either the residual or its negative, depending on the user's 114 preference, as set with SNESSetFunction(). 115 116 .keywords: SNES, nonlinear, default, monitor, residual, norm 117 118 .seealso: SNESSetMonitor() 119 @*/ 120 int SNESDefaultMonitor(SNES snes,int its, double fnorm,void *dummy) 121 { 122 fprintf( stdout, "iter = %d, residual norm %g \n",its,fnorm); 123 return 0; 124 } 125 126 /*@ 127 SNESDefaultConverged - Default test for monitoring the convergence 128 of the SNES solvers. 129 130 Input Parameters: 131 . snes - the SNES context 132 . xnorm - 2-norm of current iterate 133 . pnorm - 2-norm of current step 134 . fnorm - 2-norm of residual 135 . dummy - unused context 136 137 Returns: 138 $ 2 if ( fnorm < atol ), 139 $ 3 if ( pnorm < xtol*xnorm ), 140 $ -2 if ( nres > max_res ), 141 $ 0 otherwise, 142 143 where 144 $ atol - absolute residual norm tolerance, 145 $ set with SNESSetAbsoluteTolerance() 146 $ max_res - maximum number of residual evaluations, 147 $ set with SNESSetMaxResidualEvaluations() 148 $ nres - number of residual evaluations 149 $ xtol - relative residual norm tolerance, 150 $ set with SNESSetRelativeTolerance() 151 152 .keywords: SNES, nonlinear, default, converged, convergence 153 154 .seealso: SNESSetConvergenceTest() 155 @*/ 156 int SNESDefaultConverged(SNES snes,double xnorm,double pnorm,double fnorm, 157 void *dummy) 158 { 159 if (fnorm < snes->atol) { 160 PLogInfo((PetscObject)snes, 161 "Converged due to absolute residual norm %g < %g\n",fnorm,snes->atol); 162 return 2; 163 } 164 if (pnorm < snes->xtol*(xnorm)) { 165 PLogInfo((PetscObject)snes, 166 "Converged due to small update length %g < %g*%g\n", 167 pnorm,snes->xtol,xnorm); 168 return 3; 169 } 170 if (snes->nresids > snes->max_resids) { 171 PLogInfo((PetscObject)snes, 172 "Exceeded maximum number of residual evaluations: %d > %d\n", 173 snes->nresids, snes->max_resids ); 174 return -2; 175 } 176 return 0; 177 } 178 179 /* ------------------------------------------------------------ */ 180 /*ARGSUSED*/ 181 /*@ 182 SNESNoLineSearch - This routine is not a line search at all; 183 it simply uses the full Newton step. Thus, this routine is intended 184 to serve as a template and is not recommended for general use. 185 186 Input Parameters: 187 . snes - nonlinear context 188 . x - current iterate 189 . f - residual evaluated at x 190 . y - search direction (contains new iterate on output) 191 . w - work vector 192 . fnorm - 2-norm of f 193 194 Output parameters: 195 . g - residual evaluated at new iterate y 196 . y - new iterate (contains search direction on input) 197 . gnorm - 2-norm of g 198 . ynorm - 2-norm of search length 199 200 Returns: 201 1, indicating success of the step. 202 203 .keywords: SNES, nonlinear, line search, cubic 204 205 .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), 206 .seealso: SNESSetLineSearchRoutine() 207 @*/ 208 int SNESNoLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 209 double fnorm, double *ynorm, double *gnorm ) 210 { 211 int ierr; 212 Scalar one = 1.0; 213 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 214 VecNorm(y, ynorm ); /* ynorm = || y || */ 215 VecAXPY(&one, x, y ); /* y <- x + y */ 216 ierr = SNESComputeFunction(snes,y,g); CHKERR(ierr); 217 VecNorm( g, gnorm ); /* gnorm = || g || */ 218 PLogEventEnd(SNES_LineSearch,snes,x,f,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 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 266 alpha = neP->alpha; 267 maxstep = neP->maxstep; 268 steptol = neP->steptol; 269 270 VecNorm(y, ynorm ); 271 if (*ynorm > maxstep) { /* Step too big, so scale back */ 272 scale = maxstep/(*ynorm); 273 #if defined(PETSC_COMPLEX) 274 PLogInfo((PetscObject)snes,"Scaling step by %g\n",real(scale)); 275 #else 276 PLogInfo((PetscObject)snes,"Scaling step by %g\n",scale); 277 #endif 278 VecScale(&scale, y ); 279 *ynorm = maxstep; 280 } 281 minlambda = steptol/(*ynorm); 282 #if defined(PETSC_COMPLEX) 283 VecDot(f, y, &cinitslope ); 284 initslope = real(cinitslope); 285 #else 286 VecDot(f, y, &initslope ); 287 #endif 288 if (initslope > 0.0) initslope = -initslope; 289 if (initslope == 0.0) initslope = -1.0; 290 291 VecCopy(y, w ); 292 VecAXPY(&one, x, w ); 293 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 294 VecNorm(g, gnorm ); 295 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 296 VecCopy(w, y ); 297 PLogInfo((PetscObject)snes,"Using full step\n"); 298 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 299 return 0; 300 } 301 302 /* Fit points with quadratic */ 303 lambda = 1.0; count = 0; 304 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 305 lambdaprev = lambda; 306 gnormprev = *gnorm; 307 if (lambdatemp <= .1*lambda) { 308 lambda = .1*lambda; 309 } else lambda = lambdatemp; 310 VecCopy(x, w ); 311 #if defined(PETSC_COMPLEX) 312 clambda = lambda; VecAXPY(&clambda, y, w ); 313 #else 314 VecAXPY(&lambda, y, w ); 315 #endif 316 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 317 VecNorm(g, gnorm ); 318 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 319 VecCopy(w, y ); 320 PLogInfo((PetscObject)snes,"Quadratically determined step, lambda %g\n",lambda); 321 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 322 return 0; 323 } 324 325 /* Fit points with cubic */ 326 count = 1; 327 while (1) { 328 if (lambda <= minlambda) { /* bad luck; use full step */ 329 PLogInfo((PetscObject)snes,"Unable to find good step length! %d \n",count); 330 PLogInfo((PetscObject)snes, "f %g fnew %g ynorm %g lambda %g \n", 331 fnorm,*gnorm, *ynorm,lambda); 332 VecCopy(w, y ); 333 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 334 return 0; 335 } 336 t1 = *gnorm - fnorm - lambda*initslope; 337 t2 = gnormprev - fnorm - lambdaprev*initslope; 338 a = (t1/(lambda*lambda) - 339 t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 340 b = (-lambdaprev*t1/(lambda*lambda) + 341 lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 342 d = b*b - 3*a*initslope; 343 if (d < 0.0) d = 0.0; 344 if (a == 0.0) { 345 lambdatemp = -initslope/(2.0*b); 346 } else { 347 lambdatemp = (-b + sqrt(d))/(3.0*a); 348 } 349 if (lambdatemp > .5*lambda) { 350 lambdatemp = .5*lambda; 351 } 352 lambdaprev = lambda; 353 gnormprev = *gnorm; 354 if (lambdatemp <= .1*lambda) { 355 lambda = .1*lambda; 356 } 357 else lambda = lambdatemp; 358 VecCopy( x, w ); 359 #if defined(PETSC_COMPLEX) 360 VecAXPY(&clambda, y, w ); 361 #else 362 VecAXPY(&lambda, y, w ); 363 #endif 364 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 365 VecNorm(g, gnorm ); 366 if (*gnorm <= fnorm + alpha*initslope) { /* is reduction enough */ 367 VecCopy(w, y ); 368 PLogInfo((PetscObject)snes,"Cubically determined step, lambda %g\n",lambda); 369 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 370 return 0; 371 } 372 count++; 373 } 374 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 375 return 0; 376 } 377 /*@ 378 SNESQuadraticLineSearch - This routine performs a cubic line search. 379 380 Input Parameters: 381 . snes - the SNES context 382 . x - current iterate 383 . f - residual evaluated at x 384 . y - search direction (contains new iterate on output) 385 . w - work vector 386 . fnorm - 2-norm of f 387 388 Output parameters: 389 . g - residual evaluated at new iterate y 390 . y - new iterate (contains search direction on input) 391 . gnorm - 2-norm of g 392 . ynorm - 2-norm of search length 393 394 Returns: 395 1 if the line search succeeds; 0 if the line search fails. 396 397 Notes: 398 Use SNESSetLineSearchRoutine() 399 to set this routine within the SNES_NLS method. 400 401 .keywords: SNES, nonlinear, quadratic, line search 402 403 .seealso: SNESCubicLineSearch(), SNESNoLineSearch(), SNESSetLineSearchRoutine() 404 @*/ 405 int SNESQuadraticLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 406 double fnorm, double *ynorm, double *gnorm ) 407 { 408 double steptol, initslope; 409 double lambdaprev, gnormprev; 410 #if defined(PETSC_COMPLEX) 411 Scalar cinitslope,clambda; 412 #endif 413 int ierr,count; 414 SNES_LS *neP = (SNES_LS *) snes->data; 415 Scalar one = 1.0,scale; 416 double maxstep,minlambda,alpha,lambda,lambdatemp; 417 418 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 419 alpha = neP->alpha; 420 maxstep = neP->maxstep; 421 steptol = neP->steptol; 422 423 VecNorm(y, ynorm ); 424 if (*ynorm > maxstep) { /* Step too big, so scale back */ 425 scale = maxstep/(*ynorm); 426 VecScale(&scale, y ); 427 *ynorm = maxstep; 428 } 429 minlambda = steptol/(*ynorm); 430 #if defined(PETSC_COMPLEX) 431 VecDot(f, y, &cinitslope ); 432 initslope = real(cinitslope); 433 #else 434 VecDot(f, y, &initslope ); 435 #endif 436 if (initslope > 0.0) initslope = -initslope; 437 if (initslope == 0.0) initslope = -1.0; 438 439 VecCopy(y, w ); 440 VecAXPY(&one, x, w ); 441 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 442 VecNorm(g, gnorm ); 443 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 444 VecCopy(w, y ); 445 PLogInfo((PetscObject)snes,"Using full step\n"); 446 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 447 return 0; 448 } 449 450 /* Fit points with quadratic */ 451 lambda = 1.0; count = 0; 452 count = 1; 453 while (1) { 454 if (lambda <= minlambda) { /* bad luck; use full step */ 455 PLogInfo((PetscObject)snes,"Unable to find good step length! %d \n",count); 456 PLogInfo((PetscObject)snes, "f %g fnew %g ynorm %g lambda %g \n", 457 fnorm,*gnorm, *ynorm,lambda); 458 VecCopy(w, y ); 459 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 460 return 0; 461 } 462 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 463 lambdaprev = lambda; 464 gnormprev = *gnorm; 465 if (lambdatemp <= .1*lambda) { 466 lambda = .1*lambda; 467 } else lambda = lambdatemp; 468 VecCopy(x, w ); 469 #if defined(PETSC_COMPLEX) 470 clambda = lambda; VecAXPY(&clambda, y, w ); 471 #else 472 VecAXPY(&lambda, y, w ); 473 #endif 474 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 475 VecNorm(g, gnorm ); 476 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 477 VecCopy(w, y ); 478 PLogInfo((PetscObject)snes,"Quadratically determined step, lambda %g\n",lambda); 479 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 480 return 0; 481 } 482 count++; 483 } 484 485 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 486 return 0; 487 } 488 /* ------------------------------------------------------------ */ 489 /*@C 490 SNESSetLineSearchRoutine - Sets the line search routine to be used 491 by the method SNES_LS. 492 493 Input Parameters: 494 . snes - nonlinear context obtained from SNESCreate() 495 . func - pointer to int function 496 497 Possible routines: 498 . SNESCubicLineSearch() - default line search 499 . SNESQuadraticLineSearch() - quadratic line search 500 . SNESNoLineSearch() - the full Newton step (actually not a line search) 501 502 Calling sequence of func: 503 func (SNES snes, Vec x, Vec f, Vec g, Vec y, 504 Vec w, double fnorm, double *ynorm, double *gnorm) 505 506 Input parameters for func: 507 . snes - nonlinear context 508 . x - current iterate 509 . f - residual evaluated at x 510 . y - search direction (contains new iterate on output) 511 . w - work vector 512 . fnorm - 2-norm of f 513 514 Output parameters for func: 515 . g - residual evaluated at new iterate y 516 . y - new iterate (contains search direction on input) 517 . gnorm - 2-norm of g 518 . ynorm - 2-norm of search length 519 520 Returned by func: 521 1 if the line search succeeds; 0 if the line search fails. 522 523 .keywords: SNES, nonlinear, set, line search, routine 524 525 .seealso: SNESNoLineSearch(), SNESQuadraticLineSearch(), SNESCubicLineSearch() 526 @*/ 527 int SNESSetLineSearchRoutine(SNES snes,int (*func)(SNES,Vec,Vec,Vec,Vec,Vec, 528 double,double *,double*) ) 529 { 530 if ((snes)->type == SNES_NLS) 531 ((SNES_LS *)(snes->data))->LineSearch = func; 532 return 0; 533 } 534 535 static int SNESPrintHelp_LS(SNES snes) 536 { 537 SNES_LS *ls = (SNES_LS *)snes->data; 538 fprintf(stderr,"-snes_line_search [basic,quadratic,cubic]\n"); 539 fprintf(stderr,"-snes_line_search_alpha alpha (default %g)\n",ls->alpha); 540 fprintf(stderr,"-snes_line_search_maxstep max (default %g)\n",ls->maxstep); 541 fprintf(stderr,"-snes_line_search_steptol tol (default %g)\n",ls->steptol); 542 return 0; 543 } 544 545 #include "options.h" 546 static int SNESSetFromOptions_LS(SNES snes) 547 { 548 SNES_LS *ls = (SNES_LS *)snes->data; 549 char ver[16]; 550 double tmp; 551 552 if (OptionsGetDouble(0,snes->prefix,"-snes_line_search_alpa",&tmp)) { 553 ls->alpha = tmp; 554 } 555 if (OptionsGetDouble(0,snes->prefix,"-snes_line_search_maxstep",&tmp)) { 556 ls->maxstep = tmp; 557 } 558 if (OptionsGetDouble(0,snes->prefix,"-snes_line_search_steptol",&tmp)) { 559 ls->steptol = tmp; 560 } 561 if (OptionsGetString(0,snes->prefix,"-snes_line_search",ver,16)) { 562 if (!strcmp(ver,"basic")) { 563 SNESSetLineSearchRoutine(snes,SNESNoLineSearch); 564 } 565 else if (!strcmp(ver,"quadratic")) { 566 SNESSetLineSearchRoutine(snes,SNESQuadraticLineSearch); 567 } 568 else if (!strcmp(ver,"cubic")) { 569 SNESSetLineSearchRoutine(snes,SNESCubicLineSearch); 570 } 571 else {SETERR(1,"Unknown line search?");} 572 } 573 return 0; 574 } 575 576 /* ------------------------------------------------------------ */ 577 int SNESCreate_LS(SNES snes ) 578 { 579 SNES_LS *neP; 580 581 snes->type = SNES_NLS; 582 snes->Setup = SNESSetUp_LS; 583 snes->Solver = SNESSolve_LS; 584 snes->destroy = SNESDestroy_LS; 585 snes->Converged = SNESDefaultConverged; 586 snes->PrintHelp = SNESPrintHelp_LS; 587 snes->SetFromOptions = SNESSetFromOptions_LS; 588 589 neP = NEW(SNES_LS); CHKPTR(neP); 590 snes->data = (void *) neP; 591 neP->alpha = 1.e-4; 592 neP->maxstep = 1.e8; 593 neP->steptol = 1.e-12; 594 neP->LineSearch = SNESCubicLineSearch; 595 return 0; 596 } 597 598 599 600 601