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