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