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