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