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