1 #ifndef lint 2 static char vcid[] = "$Id: ls.c,v 1.11 1995/04/25 16:24:36 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,flg; 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)(snes,X,&snes->jacobian,&snes->jacobian_pre, 58 &flg,snes->jacP); 59 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg); 60 ierr = SLESSolve(snes->sles,F,Y,&lits); CHKERR(ierr); 61 ierr = (*neP->LineSearch)(snes, X, F, G, Y, W, fnorm, &ynorm, &gnorm ); 62 CHKERR(ierr); 63 64 TMP = F; F = G; G = TMP; 65 TMP = X; X = Y; Y = TMP; 66 fnorm = gnorm; 67 68 snes->norm = fnorm; 69 if (history && history_len > i+1) history[i+1] = fnorm; 70 VecNorm( X, &xnorm ); /* xnorm = || X || */ 71 if (snes->Monitor) (*snes->Monitor)(snes,i+1,fnorm,snes->monP); 72 73 /* Test for convergence */ 74 if ((*snes->Converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 75 if (X != snes->vec_sol) VecCopy( X, snes->vec_sol ); 76 break; 77 } 78 } 79 if (i == maxits) i--; 80 *outits = i+1; 81 return 0; 82 } 83 /* ------------------------------------------------------------ */ 84 /*ARGSUSED*/ 85 int SNESSetUp_LS(SNES snes ) 86 { 87 int ierr; 88 snes->nwork = 3; 89 ierr = VecGetVecs( snes->vec_sol, snes->nwork,&snes->work ); CHKERR(ierr); 90 PLogObjectParents(snes,snes->nwork,snes->work ); 91 return 0; 92 } 93 /* ------------------------------------------------------------ */ 94 /*ARGSUSED*/ 95 int SNESDestroy_LS(PetscObject obj) 96 { 97 SNES snes = (SNES) obj; 98 SLESDestroy(snes->sles); 99 VecFreeVecs(snes->work, snes->nwork ); 100 PLogObjectDestroy(obj); 101 PETSCHEADERDESTROY(obj); 102 return 0; 103 } 104 /*@ 105 SNESDefaultMonitor - Default SNES monitoring routine. Prints the 106 residual norm at each iteration. 107 108 Input Parameters: 109 . snes - the SNES context 110 . its - iteration number 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 SNESSetFunction(). 117 118 .keywords: SNES, nonlinear, default, monitor, residual, norm 119 120 .seealso: SNESSetMonitor() 121 @*/ 122 int SNESDefaultMonitor(SNES snes,int its, 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 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 216 VecNorm(y, ynorm ); /* ynorm = || y || */ 217 VecAXPY(&one, x, y ); /* y <- x + y */ 218 ierr = SNESComputeFunction(snes,y,g); CHKERR(ierr); 219 VecNorm( g, gnorm ); /* gnorm = || g || */ 220 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 221 return 1; 222 } 223 /* ------------------------------------------------------------------ */ 224 /*@ 225 SNESCubicLineSearch - This routine performs a cubic line search and 226 is the default line search method. 227 228 Input Parameters: 229 . snes - nonlinear context 230 . x - current iterate 231 . f - residual evaluated at x 232 . y - search direction (contains new iterate on output) 233 . w - work vector 234 . fnorm - 2-norm of f 235 236 Output parameters: 237 . g - residual evaluated at new iterate y 238 . y - new iterate (contains search direction on input) 239 . gnorm - 2-norm of g 240 . ynorm - 2-norm of search length 241 242 Returns: 243 1 if the line search succeeds; 0 if the line search fails. 244 245 Notes: 246 This line search is taken from "Numerical Methods for Unconstrained 247 Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 248 249 .keywords: SNES, nonlinear, line search, cubic 250 251 .seealso: SNESNoLineSearch(), SNESNoLineSearch(), SNESSetLineSearchRoutine() 252 @*/ 253 int SNESCubicLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 254 double fnorm, double *ynorm, double *gnorm ) 255 { 256 double steptol, initslope; 257 double lambdaprev, gnormprev; 258 double a, b, d, t1, t2; 259 #if defined(PETSC_COMPLEX) 260 Scalar cinitslope,clambda; 261 #endif 262 int ierr,count; 263 SNES_LS *neP = (SNES_LS *) snes->data; 264 Scalar one = 1.0,scale; 265 double maxstep,minlambda,alpha,lambda,lambdatemp; 266 267 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 268 alpha = neP->alpha; 269 maxstep = neP->maxstep; 270 steptol = neP->steptol; 271 272 VecNorm(y, ynorm ); 273 if (*ynorm > maxstep) { /* Step too big, so scale back */ 274 scale = maxstep/(*ynorm); 275 #if defined(PETSC_COMPLEX) 276 PLogInfo((PetscObject)snes,"Scaling step by %g\n",real(scale)); 277 #else 278 PLogInfo((PetscObject)snes,"Scaling step by %g\n",scale); 279 #endif 280 VecScale(&scale, y ); 281 *ynorm = maxstep; 282 } 283 minlambda = steptol/(*ynorm); 284 #if defined(PETSC_COMPLEX) 285 VecDot(f, y, &cinitslope ); 286 initslope = real(cinitslope); 287 #else 288 VecDot(f, y, &initslope ); 289 #endif 290 if (initslope > 0.0) initslope = -initslope; 291 if (initslope == 0.0) initslope = -1.0; 292 293 VecCopy(y, w ); 294 VecAXPY(&one, x, w ); 295 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 296 VecNorm(g, gnorm ); 297 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 298 VecCopy(w, y ); 299 PLogInfo((PetscObject)snes,"Using full step\n"); 300 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 301 return 0; 302 } 303 304 /* Fit points with quadratic */ 305 lambda = 1.0; count = 0; 306 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 307 lambdaprev = lambda; 308 gnormprev = *gnorm; 309 if (lambdatemp <= .1*lambda) { 310 lambda = .1*lambda; 311 } else lambda = lambdatemp; 312 VecCopy(x, w ); 313 #if defined(PETSC_COMPLEX) 314 clambda = lambda; VecAXPY(&clambda, y, w ); 315 #else 316 VecAXPY(&lambda, y, w ); 317 #endif 318 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 319 VecNorm(g, gnorm ); 320 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 321 VecCopy(w, y ); 322 PLogInfo((PetscObject)snes,"Quadratically determined step, lambda %g\n",lambda); 323 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 324 return 0; 325 } 326 327 /* Fit points with cubic */ 328 count = 1; 329 while (1) { 330 if (lambda <= minlambda) { /* bad luck; use full step */ 331 PLogInfo((PetscObject)snes,"Unable to find good step length! %d \n",count); 332 PLogInfo((PetscObject)snes, "f %g fnew %g ynorm %g lambda %g \n", 333 fnorm,*gnorm, *ynorm,lambda); 334 VecCopy(w, y ); 335 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 336 return -1; 337 } 338 t1 = *gnorm - fnorm - lambda*initslope; 339 t2 = gnormprev - fnorm - lambdaprev*initslope; 340 a = (t1/(lambda*lambda) - 341 t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 342 b = (-lambdaprev*t1/(lambda*lambda) + 343 lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 344 d = b*b - 3*a*initslope; 345 if (d < 0.0) d = 0.0; 346 if (a == 0.0) { 347 lambdatemp = -initslope/(2.0*b); 348 } else { 349 lambdatemp = (-b + sqrt(d))/(3.0*a); 350 } 351 if (lambdatemp > .5*lambda) { 352 lambdatemp = .5*lambda; 353 } 354 lambdaprev = lambda; 355 gnormprev = *gnorm; 356 if (lambdatemp <= .1*lambda) { 357 lambda = .1*lambda; 358 } 359 else lambda = lambdatemp; 360 VecCopy( x, w ); 361 #if defined(PETSC_COMPLEX) 362 VecAXPY(&clambda, y, w ); 363 #else 364 VecAXPY(&lambda, y, w ); 365 #endif 366 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 367 VecNorm(g, gnorm ); 368 if (*gnorm <= fnorm + alpha*initslope) { /* is reduction enough */ 369 VecCopy(w, y ); 370 PLogInfo((PetscObject)snes,"Cubically determined step, lambda %g\n",lambda); 371 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 372 return 0; 373 } 374 count++; 375 } 376 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 377 return 0; 378 } 379 /*@ 380 SNESQuadraticLineSearch - This routine performs a cubic line search. 381 382 Input Parameters: 383 . snes - the SNES context 384 . x - current iterate 385 . f - residual evaluated at x 386 . y - search direction (contains new iterate on output) 387 . w - work vector 388 . fnorm - 2-norm of f 389 390 Output parameters: 391 . g - residual evaluated at new iterate y 392 . y - new iterate (contains search direction on input) 393 . gnorm - 2-norm of g 394 . ynorm - 2-norm of search length 395 396 Returns: 397 1 if the line search succeeds; 0 if the line search fails. 398 399 Notes: 400 Use SNESSetLineSearchRoutine() 401 to set this routine within the SNES_NLS method. 402 403 .keywords: SNES, nonlinear, quadratic, line search 404 405 .seealso: SNESCubicLineSearch(), SNESNoLineSearch(), SNESSetLineSearchRoutine() 406 @*/ 407 int SNESQuadraticLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 408 double fnorm, double *ynorm, double *gnorm ) 409 { 410 double steptol, initslope; 411 double lambdaprev, gnormprev; 412 #if defined(PETSC_COMPLEX) 413 Scalar cinitslope,clambda; 414 #endif 415 int ierr,count; 416 SNES_LS *neP = (SNES_LS *) snes->data; 417 Scalar one = 1.0,scale; 418 double maxstep,minlambda,alpha,lambda,lambdatemp; 419 420 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 421 alpha = neP->alpha; 422 maxstep = neP->maxstep; 423 steptol = neP->steptol; 424 425 VecNorm(y, ynorm ); 426 if (*ynorm > maxstep) { /* Step too big, so scale back */ 427 scale = maxstep/(*ynorm); 428 VecScale(&scale, y ); 429 *ynorm = maxstep; 430 } 431 minlambda = steptol/(*ynorm); 432 #if defined(PETSC_COMPLEX) 433 VecDot(f, y, &cinitslope ); 434 initslope = real(cinitslope); 435 #else 436 VecDot(f, y, &initslope ); 437 #endif 438 if (initslope > 0.0) initslope = -initslope; 439 if (initslope == 0.0) initslope = -1.0; 440 441 VecCopy(y, w ); 442 VecAXPY(&one, x, w ); 443 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 444 VecNorm(g, gnorm ); 445 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 446 VecCopy(w, y ); 447 PLogInfo((PetscObject)snes,"Using full step\n"); 448 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 449 return 0; 450 } 451 452 /* Fit points with quadratic */ 453 lambda = 1.0; count = 0; 454 count = 1; 455 while (1) { 456 if (lambda <= minlambda) { /* bad luck; use full step */ 457 PLogInfo((PetscObject)snes,"Unable to find good step length! %d \n",count); 458 PLogInfo((PetscObject)snes, "f %g fnew %g ynorm %g lambda %g \n", 459 fnorm,*gnorm, *ynorm,lambda); 460 VecCopy(w, y ); 461 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 462 return 0; 463 } 464 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 465 lambdaprev = lambda; 466 gnormprev = *gnorm; 467 if (lambdatemp <= .1*lambda) { 468 lambda = .1*lambda; 469 } else lambda = lambdatemp; 470 VecCopy(x, w ); 471 #if defined(PETSC_COMPLEX) 472 clambda = lambda; VecAXPY(&clambda, y, w ); 473 #else 474 VecAXPY(&lambda, y, w ); 475 #endif 476 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 477 VecNorm(g, gnorm ); 478 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 479 VecCopy(w, y ); 480 PLogInfo((PetscObject)snes,"Quadratically determined step, lambda %g\n",lambda); 481 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 482 return 0; 483 } 484 count++; 485 } 486 487 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 488 return 0; 489 } 490 /* ------------------------------------------------------------ */ 491 /*@C 492 SNESSetLineSearchRoutine - Sets the line search routine to be used 493 by the method SNES_LS. 494 495 Input Parameters: 496 . snes - nonlinear context obtained from SNESCreate() 497 . func - pointer to int function 498 499 Possible routines: 500 . SNESCubicLineSearch() - default line search 501 . SNESQuadraticLineSearch() - quadratic line search 502 . SNESNoLineSearch() - the full Newton step (actually not a line search) 503 504 Calling sequence of func: 505 func (SNES snes, Vec x, Vec f, Vec g, Vec y, 506 Vec w, double fnorm, double *ynorm, double *gnorm) 507 508 Input parameters for func: 509 . snes - nonlinear context 510 . x - current iterate 511 . f - residual evaluated at x 512 . y - search direction (contains new iterate on output) 513 . w - work vector 514 . fnorm - 2-norm of f 515 516 Output parameters for func: 517 . g - residual evaluated at new iterate y 518 . y - new iterate (contains search direction on input) 519 . gnorm - 2-norm of g 520 . ynorm - 2-norm of search length 521 522 Returned by func: 523 1 if the line search succeeds; 0 if the line search fails. 524 525 .keywords: SNES, nonlinear, set, line search, routine 526 527 .seealso: SNESNoLineSearch(), SNESQuadraticLineSearch(), SNESCubicLineSearch() 528 @*/ 529 int SNESSetLineSearchRoutine(SNES snes,int (*func)(SNES,Vec,Vec,Vec,Vec,Vec, 530 double,double *,double*) ) 531 { 532 if ((snes)->type == SNES_NLS) 533 ((SNES_LS *)(snes->data))->LineSearch = func; 534 return 0; 535 } 536 537 static int SNESPrintHelp_LS(SNES snes) 538 { 539 SNES_LS *ls = (SNES_LS *)snes->data; 540 fprintf(stderr,"-snes_line_search [basic,quadratic,cubic]\n"); 541 fprintf(stderr,"-snes_line_search_alpha alpha (default %g)\n",ls->alpha); 542 fprintf(stderr,"-snes_line_search_maxstep max (default %g)\n",ls->maxstep); 543 fprintf(stderr,"-snes_line_search_steptol tol (default %g)\n",ls->steptol); 544 return 0; 545 } 546 547 #include "options.h" 548 static int SNESSetFromOptions_LS(SNES snes) 549 { 550 SNES_LS *ls = (SNES_LS *)snes->data; 551 char ver[16]; 552 double tmp; 553 554 if (OptionsGetDouble(0,snes->prefix,"-snes_line_search_alpa",&tmp)) { 555 ls->alpha = tmp; 556 } 557 if (OptionsGetDouble(0,snes->prefix,"-snes_line_search_maxstep",&tmp)) { 558 ls->maxstep = tmp; 559 } 560 if (OptionsGetDouble(0,snes->prefix,"-snes_line_search_steptol",&tmp)) { 561 ls->steptol = tmp; 562 } 563 if (OptionsGetString(0,snes->prefix,"-snes_line_search",ver,16)) { 564 if (!strcmp(ver,"basic")) { 565 SNESSetLineSearchRoutine(snes,SNESNoLineSearch); 566 } 567 else if (!strcmp(ver,"quadratic")) { 568 SNESSetLineSearchRoutine(snes,SNESQuadraticLineSearch); 569 } 570 else if (!strcmp(ver,"cubic")) { 571 SNESSetLineSearchRoutine(snes,SNESCubicLineSearch); 572 } 573 else {SETERR(1,"Unknown line search?");} 574 } 575 return 0; 576 } 577 578 /* ------------------------------------------------------------ */ 579 int SNESCreate_LS(SNES snes ) 580 { 581 SNES_LS *neP; 582 583 snes->type = SNES_NLS; 584 snes->Setup = SNESSetUp_LS; 585 snes->Solver = SNESSolve_LS; 586 snes->destroy = SNESDestroy_LS; 587 snes->Converged = SNESDefaultConverged; 588 snes->PrintHelp = SNESPrintHelp_LS; 589 snes->SetFromOptions = SNESSetFromOptions_LS; 590 591 neP = NEW(SNES_LS); CHKPTR(neP); 592 snes->data = (void *) neP; 593 neP->alpha = 1.e-4; 594 neP->maxstep = 1.e8; 595 neP->steptol = 1.e-12; 596 neP->LineSearch = SNESCubicLineSearch; 597 return 0; 598 } 599 600 601 602 603