1 #ifndef lint 2 static char vcid[] = "$Id: ls.c,v 1.35 1995/07/27 03:01:17 curfman Exp curfman $"; 3 #endif 4 5 #include <math.h> 6 #include "ls.h" 7 #include "pviewer.h" 8 #if defined(HAVE_STRING_H) 9 #include <string.h> 10 #endif 11 12 /* 13 Implements a line search variant of Newton's Method 14 for solving systems of nonlinear equations. 15 16 Input parameters: 17 . snes - nonlinear context obtained from SNESCreate() 18 19 Output Parameters: 20 . outits - Number of global iterations until termination. 21 22 Notes: 23 This implements essentially a truncated Newton method with a 24 line search. By default a cubic backtracking line search 25 is employed, as described in the text "Numerical Methods for 26 Unconstrained Optimization and Nonlinear Equations" by Dennis 27 and Schnabel. See the examples in src/snes/examples. 28 */ 29 30 int SNESSolve_LS(SNES snes,int *outits) 31 { 32 SNES_LS *neP = (SNES_LS *) snes->data; 33 int maxits, i, history_len, ierr, lits, lsfail; 34 MatStructure flg = ALLMAT_DIFFERENT_NONZERO_PATTERN; 35 double fnorm, gnorm, xnorm, ynorm, *history; 36 Vec Y, X, F, G, W, TMP; 37 SLES sles; 38 KSP ksp; 39 40 history = snes->conv_hist; /* convergence history */ 41 history_len = snes->conv_hist_len; /* convergence history length */ 42 maxits = snes->max_its; /* maximum number of iterations */ 43 X = snes->vec_sol; /* solution vector */ 44 F = snes->vec_func; /* residual vector */ 45 Y = snes->work[0]; /* work vectors */ 46 G = snes->work[1]; 47 W = snes->work[2]; 48 49 ierr = SNESComputeInitialGuess(snes,X); CHKERRQ(ierr); /* X <- X_0 */ 50 VecNorm(X,&xnorm); /* xnorm = || X || */ 51 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* (+/-) F(X) */ 52 VecNorm(F,&fnorm); /* fnorm <- ||F|| */ 53 snes->norm = fnorm; 54 if (history && history_len > 0) history[0] = fnorm; 55 if (snes->monitor) 56 {ierr = (*snes->monitor)(snes,0,fnorm,snes->monP); CHKERRQ(ierr);} 57 58 /* Set the KSP stopping criteria to use the Eisenstat-Walker method */ 59 if (snes->ksp_ewconv) { 60 ierr = SNESGetSLES(snes,&sles); CHKERRQ(ierr); 61 ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr); 62 ierr = KSPSetConvergenceTest(ksp,SNES_KSP_EW_Converged_Private, 63 (void *)snes); CHKERRQ(ierr); 64 } 65 66 for ( i=0; i<maxits; i++ ) { 67 snes->iter = i+1; 68 69 /* Solve J Y = -F, where J is Jacobian matrix */ 70 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre, 71 &flg); CHKERRQ(ierr); 72 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre, 73 flg); CHKERRQ(ierr); 74 ierr = SLESSolve(snes->sles,F,Y,&lits); CHKERRQ(ierr); 75 ierr = VecCopy(Y,snes->vec_sol_update_always); CHKERRQ(ierr); 76 ierr = (*neP->LineSearch)(snes,X,F,G,Y,W,fnorm,&ynorm,&gnorm,&lsfail); 77 if (lsfail) snes->nfailures++; 78 CHKERRQ(ierr); 79 80 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 81 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 82 fnorm = gnorm; 83 84 snes->norm = fnorm; 85 if (history && history_len > i+1) history[i+1] = fnorm; 86 VecNorm(X,&xnorm); /* xnorm = || X || */ 87 if (snes->monitor) 88 {(*snes->monitor)(snes,i+1,fnorm,snes->monP); CHKERRQ(ierr);} 89 90 /* Test for convergence */ 91 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 92 if (X != snes->vec_sol) { 93 VecCopy(X,snes->vec_sol); 94 snes->vec_sol_always = snes->vec_sol; 95 snes->vec_func_always = snes->vec_func; 96 } 97 break; 98 } 99 } 100 if (i == maxits) { 101 PLogInfo((PetscObject)snes, 102 "Maximum number of iterations has been reached: %d\n",maxits); 103 i--; 104 } 105 *outits = i+1; 106 return 0; 107 } 108 /* ------------------------------------------------------------ */ 109 int SNESSetUp_LS(SNES snes ) 110 { 111 int ierr; 112 snes->nwork = 4; 113 ierr = VecGetVecs(snes->vec_sol,snes->nwork,&snes->work); CHKERRQ(ierr); 114 PLogObjectParents(snes,snes->nwork,snes->work ); 115 snes->vec_sol_update_always = snes->work[3]; 116 return 0; 117 } 118 /* ------------------------------------------------------------ */ 119 int SNESDestroy_LS(PetscObject obj) 120 { 121 SNES snes = (SNES) obj; 122 VecFreeVecs(snes->work,snes->nwork); 123 PETSCFREE(snes->data); 124 return 0; 125 } 126 /* ------------------------------------------------------------ */ 127 /*ARGSUSED*/ 128 /*@ 129 SNESNoLineSearch - This routine is not a line search at all; 130 it simply uses the full Newton step. Thus, this routine is intended 131 to serve as a template and is not recommended for general use. 132 133 Input Parameters: 134 . snes - nonlinear context 135 . x - current iterate 136 . f - residual evaluated at x 137 . y - search direction (contains new iterate on output) 138 . w - work vector 139 . fnorm - 2-norm of f 140 141 Output Parameters: 142 . g - residual evaluated at new iterate y 143 . y - new iterate (contains search direction on input) 144 . gnorm - 2-norm of g 145 . ynorm - 2-norm of search length 146 . flag - set to 0, indicating a successful line search 147 148 Options Database Key: 149 $ -snes_line_search basic 150 151 .keywords: SNES, nonlinear, line search, cubic 152 153 .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), 154 SNESSetLineSearchRoutine() 155 @*/ 156 int SNESNoLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 157 double fnorm, double *ynorm, double *gnorm,int *flag ) 158 { 159 int ierr; 160 Scalar one = 1.0; 161 *flag = 0; 162 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 163 VecNorm(y,ynorm); /* ynorm = || y || */ 164 VecAXPY(&one,x,y); /* y <- x + y */ 165 ierr = SNESComputeFunction(snes,y,g); CHKERRQ(ierr); 166 VecNorm(g,gnorm); /* gnorm = || g || */ 167 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 168 return 0; 169 } 170 /* ------------------------------------------------------------------ */ 171 /*@ 172 SNESCubicLineSearch - This routine performs a cubic line search and 173 is the default line search method. 174 175 Input Parameters: 176 . snes - nonlinear context 177 . x - current iterate 178 . f - residual evaluated at x 179 . y - search direction (contains new iterate on output) 180 . w - work vector 181 . fnorm - 2-norm of f 182 183 Output parameters: 184 . g - residual evaluated at new iterate y 185 . y - new iterate (contains search direction on input) 186 . gnorm - 2-norm of g 187 . ynorm - 2-norm of search length 188 . flag - 0 if line search succeeds; -1 on failure. 189 190 Options Database Key: 191 $ -snes_line_search cubic 192 193 Notes: 194 This line search is taken from "Numerical Methods for Unconstrained 195 Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 196 197 .keywords: SNES, nonlinear, line search, cubic 198 199 .seealso: SNESNoLineSearch(), SNESNoLineSearch(), SNESSetLineSearchRoutine() 200 @*/ 201 int SNESCubicLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 202 double fnorm, double *ynorm, double *gnorm,int *flag) 203 { 204 double steptol, initslope; 205 double lambdaprev, gnormprev; 206 double a, b, d, t1, t2; 207 #if defined(PETSC_COMPLEX) 208 Scalar cinitslope,clambda; 209 #endif 210 int ierr,count; 211 SNES_LS *neP = (SNES_LS *) snes->data; 212 Scalar one = 1.0,scale; 213 double maxstep,minlambda,alpha,lambda,lambdatemp; 214 215 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 216 *flag = 0; 217 alpha = neP->alpha; 218 maxstep = neP->maxstep; 219 steptol = neP->steptol; 220 221 VecNorm(y, ynorm ); 222 if (*ynorm > maxstep) { /* Step too big, so scale back */ 223 scale = maxstep/(*ynorm); 224 #if defined(PETSC_COMPLEX) 225 PLogInfo((PetscObject)snes,"Scaling step by %g\n",real(scale)); 226 #else 227 PLogInfo((PetscObject)snes,"Scaling step by %g\n",scale); 228 #endif 229 VecScale(&scale, y ); 230 *ynorm = maxstep; 231 } 232 minlambda = steptol/(*ynorm); 233 #if defined(PETSC_COMPLEX) 234 VecDot(f, y, &cinitslope ); 235 initslope = real(cinitslope); 236 #else 237 VecDot(f, y, &initslope ); 238 #endif 239 if (initslope > 0.0) initslope = -initslope; 240 if (initslope == 0.0) initslope = -1.0; 241 242 VecCopy(y, w ); 243 VecAXPY(&one, x, w ); 244 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 245 VecNorm(g, gnorm ); 246 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 247 VecCopy(w, y ); 248 PLogInfo((PetscObject)snes,"Using full step\n"); 249 goto theend; 250 } 251 252 /* Fit points with quadratic */ 253 lambda = 1.0; count = 0; 254 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 255 lambdaprev = lambda; 256 gnormprev = *gnorm; 257 if (lambdatemp <= .1*lambda) { 258 lambda = .1*lambda; 259 } else lambda = lambdatemp; 260 VecCopy(x, w ); 261 #if defined(PETSC_COMPLEX) 262 clambda = lambda; VecAXPY(&clambda, y, w ); 263 #else 264 VecAXPY(&lambda, y, w ); 265 #endif 266 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 267 VecNorm(g, gnorm ); 268 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 269 VecCopy(w, y ); 270 PLogInfo((PetscObject)snes,"Quadratically determined step, lambda %g\n",lambda); 271 goto theend; 272 } 273 274 /* Fit points with cubic */ 275 count = 1; 276 while (1) { 277 if (lambda <= minlambda) { /* bad luck; use full step */ 278 PLogInfo((PetscObject)snes,"Unable to find good step length! %d \n",count); 279 PLogInfo((PetscObject)snes, "f %g fnew %g ynorm %g lambda %g \n", 280 fnorm,*gnorm, *ynorm,lambda); 281 VecCopy(w, y ); 282 *flag = -1; break; 283 } 284 t1 = *gnorm - fnorm - lambda*initslope; 285 t2 = gnormprev - fnorm - lambdaprev*initslope; 286 a = (t1/(lambda*lambda) - 287 t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 288 b = (-lambdaprev*t1/(lambda*lambda) + 289 lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 290 d = b*b - 3*a*initslope; 291 if (d < 0.0) d = 0.0; 292 if (a == 0.0) { 293 lambdatemp = -initslope/(2.0*b); 294 } else { 295 lambdatemp = (-b + sqrt(d))/(3.0*a); 296 } 297 if (lambdatemp > .5*lambda) { 298 lambdatemp = .5*lambda; 299 } 300 lambdaprev = lambda; 301 gnormprev = *gnorm; 302 if (lambdatemp <= .1*lambda) { 303 lambda = .1*lambda; 304 } 305 else lambda = lambdatemp; 306 VecCopy( x, w ); 307 #if defined(PETSC_COMPLEX) 308 VecAXPY(&clambda, y, w ); 309 #else 310 VecAXPY(&lambda, y, w ); 311 #endif 312 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 313 VecNorm(g, gnorm ); 314 if (*gnorm <= fnorm + alpha*initslope) { /* is reduction enough */ 315 VecCopy(w, y ); 316 PLogInfo((PetscObject)snes,"Cubically determined step, lambda %g\n",lambda); 317 *flag = -1; break; 318 } 319 count++; 320 } 321 theend: 322 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 323 return 0; 324 } 325 /* ------------------------------------------------------------------ */ 326 /*@ 327 SNESQuadraticLineSearch - This routine performs a cubic line search. 328 329 Input Parameters: 330 . snes - the SNES context 331 . x - current iterate 332 . f - residual evaluated at x 333 . y - search direction (contains new iterate on output) 334 . w - work vector 335 . fnorm - 2-norm of f 336 337 Output Parameters: 338 . g - residual evaluated at new iterate y 339 . y - new iterate (contains search direction on input) 340 . gnorm - 2-norm of g 341 . ynorm - 2-norm of search length 342 . flag - 0 if line search succeeds; -1 on failure. 343 344 Options Database Key: 345 $ -snes_line_search quadratic 346 347 Notes: 348 Use SNESSetLineSearchRoutine() 349 to set this routine within the SNES_NLS method. 350 351 .keywords: SNES, nonlinear, quadratic, line search 352 353 .seealso: SNESCubicLineSearch(), SNESNoLineSearch(), SNESSetLineSearchRoutine() 354 @*/ 355 int SNESQuadraticLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 356 double fnorm, double *ynorm, double *gnorm,int *flag) 357 { 358 double steptol, initslope; 359 double lambdaprev, gnormprev; 360 #if defined(PETSC_COMPLEX) 361 Scalar cinitslope,clambda; 362 #endif 363 int ierr,count; 364 SNES_LS *neP = (SNES_LS *) snes->data; 365 Scalar one = 1.0,scale; 366 double maxstep,minlambda,alpha,lambda,lambdatemp; 367 368 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 369 *flag = 0; 370 alpha = neP->alpha; 371 maxstep = neP->maxstep; 372 steptol = neP->steptol; 373 374 VecNorm(y, ynorm ); 375 if (*ynorm > maxstep) { /* Step too big, so scale back */ 376 scale = maxstep/(*ynorm); 377 VecScale(&scale, y ); 378 *ynorm = maxstep; 379 } 380 minlambda = steptol/(*ynorm); 381 #if defined(PETSC_COMPLEX) 382 VecDot(f, y, &cinitslope ); 383 initslope = real(cinitslope); 384 #else 385 VecDot(f, y, &initslope ); 386 #endif 387 if (initslope > 0.0) initslope = -initslope; 388 if (initslope == 0.0) initslope = -1.0; 389 390 VecCopy(y, w ); 391 VecAXPY(&one, x, w ); 392 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 393 VecNorm(g, gnorm ); 394 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 395 VecCopy(w, y ); 396 PLogInfo((PetscObject)snes,"Using full step\n"); 397 goto theend; 398 } 399 400 /* Fit points with quadratic */ 401 lambda = 1.0; count = 0; 402 count = 1; 403 while (1) { 404 if (lambda <= minlambda) { /* bad luck; use full step */ 405 PLogInfo((PetscObject)snes,"Unable to find good step length! %d \n",count); 406 PLogInfo((PetscObject)snes, "f %g fnew %g ynorm %g lambda %g \n", 407 fnorm,*gnorm, *ynorm,lambda); 408 VecCopy(w, y ); 409 *flag = -1; break; 410 } 411 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 412 lambdaprev = lambda; 413 gnormprev = *gnorm; 414 if (lambdatemp <= .1*lambda) { 415 lambda = .1*lambda; 416 } else lambda = lambdatemp; 417 VecCopy(x, w ); 418 #if defined(PETSC_COMPLEX) 419 clambda = lambda; VecAXPY(&clambda, y, w ); 420 #else 421 VecAXPY(&lambda, y, w ); 422 #endif 423 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 424 VecNorm(g, gnorm ); 425 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 426 VecCopy(w, y ); 427 PLogInfo((PetscObject)snes,"Quadratically determined step, lambda %g\n",lambda); 428 break; 429 } 430 count++; 431 } 432 theend: 433 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 434 return 0; 435 } 436 /* ------------------------------------------------------------ */ 437 /*@ 438 SNESSetLineSearchRoutine - Sets the line search routine to be used 439 by the method SNES_LS. 440 441 Input Parameters: 442 . snes - nonlinear context obtained from SNESCreate() 443 . func - pointer to int function 444 445 Available Routines: 446 . SNESCubicLineSearch() - default line search 447 . SNESQuadraticLineSearch() - quadratic line search 448 . SNESNoLineSearch() - the full Newton step (actually not a line search) 449 450 Options Database Keys: 451 $ -snes_line_search [basic,quadratic,cubic] 452 453 Calling sequence of func: 454 func (SNES snes, Vec x, Vec f, Vec g, Vec y, 455 Vec w, double fnorm, double *ynorm, 456 double *gnorm, *flag) 457 458 Input parameters for func: 459 . snes - nonlinear context 460 . x - current iterate 461 . f - residual evaluated at x 462 . y - search direction (contains new iterate on output) 463 . w - work vector 464 . fnorm - 2-norm of f 465 466 Output parameters for func: 467 . g - residual evaluated at new iterate y 468 . y - new iterate (contains search direction on input) 469 . gnorm - 2-norm of g 470 . ynorm - 2-norm of search length 471 . flag - set to 0 if the line search succeeds; a nonzero integer 472 on failure. 473 474 .keywords: SNES, nonlinear, set, line search, routine 475 476 .seealso: SNESNoLineSearch(), SNESQuadraticLineSearch(), SNESCubicLineSearch() 477 @*/ 478 int SNESSetLineSearchRoutine(SNES snes,int (*func)(SNES,Vec,Vec,Vec,Vec,Vec, 479 double,double *,double*,int*) ) 480 { 481 if ((snes)->type == SNES_NLS) 482 ((SNES_LS *)(snes->data))->LineSearch = func; 483 return 0; 484 } 485 /* ------------------------------------------------------------------ */ 486 static int SNESPrintHelp_LS(SNES snes) 487 { 488 SNES_LS *ls = (SNES_LS *)snes->data; 489 char *p; 490 if (snes->prefix) p = snes->prefix; else p = "-"; 491 MPIU_printf(snes->comm," method ls (system of nonlinear equations):\n"); 492 MPIU_printf(snes->comm," %ssnes_line_search [basic,quadratic,cubic]\n",p); 493 MPIU_printf(snes->comm," %ssnes_line_search_alpha alpha (default %g)\n",p,ls->alpha); 494 MPIU_printf(snes->comm," %ssnes_line_search_maxstep max (default %g)\n",p,ls->maxstep); 495 MPIU_printf(snes->comm," %ssnes_line_search_steptol tol (default %g)\n",p,ls->steptol); 496 return 0; 497 } 498 /* ------------------------------------------------------------------ */ 499 static int SNESView_LS(PetscObject obj,Viewer viewer) 500 { 501 SNES snes = (SNES)obj; 502 SNES_LS *ls = (SNES_LS *)snes->data; 503 FILE *fd = ViewerFileGetPointer_Private(viewer); 504 char *cstring; 505 506 if (ls->LineSearch == SNESNoLineSearch) 507 cstring = "SNESNoLineSearch"; 508 else if (ls->LineSearch == SNESQuadraticLineSearch) 509 cstring = "SNESQuadraticLineSearch"; 510 else if (ls->LineSearch == SNESCubicLineSearch) 511 cstring = "SNESCubicLineSearch"; 512 else 513 cstring = "unknown"; 514 MPIU_fprintf(snes->comm,fd," line search variant: %s\n",cstring); 515 MPIU_fprintf(snes->comm,fd," alpha=%g, maxstep=%g, steptol=%g\n", 516 ls->alpha,ls->maxstep,ls->steptol); 517 return 0; 518 } 519 /* ------------------------------------------------------------------ */ 520 static int SNESSetFromOptions_LS(SNES snes) 521 { 522 SNES_LS *ls = (SNES_LS *)snes->data; 523 char ver[16]; 524 double tmp; 525 526 if (OptionsGetDouble(snes->prefix,"-snes_line_search_alpa",&tmp)) { 527 ls->alpha = tmp; 528 } 529 if (OptionsGetDouble(snes->prefix,"-snes_line_search_maxstep",&tmp)) { 530 ls->maxstep = tmp; 531 } 532 if (OptionsGetDouble(snes->prefix,"-snes_line_search_steptol",&tmp)) { 533 ls->steptol = tmp; 534 } 535 if (OptionsGetString(snes->prefix,"-snes_line_search",ver,16)) { 536 if (!strcmp(ver,"basic")) { 537 SNESSetLineSearchRoutine(snes,SNESNoLineSearch); 538 } 539 else if (!strcmp(ver,"quadratic")) { 540 SNESSetLineSearchRoutine(snes,SNESQuadraticLineSearch); 541 } 542 else if (!strcmp(ver,"cubic")) { 543 SNESSetLineSearchRoutine(snes,SNESCubicLineSearch); 544 } 545 else {SETERRQ(1,"SNESSetFromOptions_LS: Unknown line search");} 546 } 547 return 0; 548 } 549 /* ------------------------------------------------------------ */ 550 int SNESCreate_LS(SNES snes ) 551 { 552 SNES_LS *neP; 553 554 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) SETERRQ(1, 555 "SNESCreate_LS: Valid for SNES_NONLINEAR_EQUATIONS problems only"); 556 snes->type = SNES_EQ_NLS; 557 snes->setup = SNESSetUp_LS; 558 snes->solve = SNESSolve_LS; 559 snes->destroy = SNESDestroy_LS; 560 snes->converged = SNESDefaultConverged; 561 snes->printhelp = SNESPrintHelp_LS; 562 snes->setfromoptions = SNESSetFromOptions_LS; 563 snes->view = SNESView_LS; 564 565 neP = PETSCNEW(SNES_LS); CHKPTRQ(neP); 566 snes->data = (void *) neP; 567 neP->alpha = 1.e-4; 568 neP->maxstep = 1.e8; 569 neP->steptol = 1.e-12; 570 neP->LineSearch = SNESCubicLineSearch; 571 return 0; 572 } 573 574 575 576 577