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