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