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