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