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