1 #ifndef lint 2 static char vcid[] = "$Id: ls.c,v 1.23 1995/06/13 01:18:08 bsmith Exp curfman $"; 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 MatStructure flg = ALLMAT_DIFFERENT_NONZERO_PATTERN; 34 double fnorm, gnorm, xnorm, ynorm, *history; 35 Vec Y, X, F, G, W, TMP; 36 37 history = snes->conv_hist; /* convergence history */ 38 history_len = snes->conv_hist_len; /* convergence history length */ 39 maxits = snes->max_its; /* maximum number of iterations */ 40 X = snes->vec_sol; /* solution vector */ 41 F = snes->vec_func; /* residual vector */ 42 Y = snes->work[0]; /* work vectors */ 43 G = snes->work[1]; 44 W = snes->work[2]; 45 46 ierr = SNESComputeInitialGuess(snes,X); CHKERRQ(ierr); /* X <- X_0 */ 47 VecNorm( X, &xnorm ); /* xnorm = || X || */ 48 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* (+/-) F(X) */ 49 VecNorm(F, &fnorm ); /* fnorm <- || F || */ 50 snes->norm = fnorm; 51 if (history && history_len > 0) history[0] = fnorm; 52 if (snes->Monitor) (*snes->Monitor)(snes,0,fnorm,snes->monP); 53 54 for ( i=0; i<maxits; i++ ) { 55 snes->iter = i+1; 56 57 /* Solve J Y = -F, where J is Jacobian matrix */ 58 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre, 59 &flg,snes->jacP); CHKERRQ(ierr); 60 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg); 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 ); 64 CHKERRQ(ierr); 65 66 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 67 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 68 fnorm = gnorm; 69 70 snes->norm = fnorm; 71 if (history && history_len > i+1) history[i+1] = fnorm; 72 VecNorm( X, &xnorm ); /* xnorm = || X || */ 73 if (snes->Monitor) (*snes->Monitor)(snes,i+1,fnorm,snes->monP); 74 75 /* Test for convergence */ 76 if ((*snes->Converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 77 if (X != snes->vec_sol) { 78 VecCopy( X, snes->vec_sol ); 79 snes->vec_sol_always = snes->vec_sol; 80 snes->vec_func_always = snes->vec_func; 81 } 82 break; 83 } 84 } 85 if (i == maxits) i--; 86 *outits = i+1; 87 return 0; 88 } 89 /* ------------------------------------------------------------ */ 90 int SNESSetUp_LS(SNES snes ) 91 { 92 int ierr; 93 snes->nwork = 4; 94 ierr = VecGetVecs( snes->vec_sol, snes->nwork,&snes->work ); CHKERRQ(ierr); 95 PLogObjectParents(snes,snes->nwork,snes->work ); 96 snes->vec_sol_update_always = snes->work[3]; 97 return 0; 98 } 99 /* ------------------------------------------------------------ */ 100 int SNESDestroy_LS(PetscObject obj) 101 { 102 SNES snes = (SNES) obj; 103 VecFreeVecs(snes->work, snes->nwork ); 104 PETSCFREE(snes->data); 105 return 0; 106 } 107 /*@ 108 SNESDefaultMonitor - Default SNES monitoring routine. Prints the 109 residual norm at each iteration. 110 111 Input Parameters: 112 . snes - the SNES context 113 . its - iteration number 114 . fnorm - 2-norm residual value (may be estimated) 115 . dummy - unused context 116 117 Notes: 118 f is either the residual or its negative, depending on the user's 119 preference, as set with SNESSetFunction(). 120 121 .keywords: SNES, nonlinear, default, monitor, residual, norm 122 123 .seealso: SNESSetMonitor() 124 @*/ 125 int SNESDefaultMonitor(SNES snes,int its, double fnorm,void *dummy) 126 { 127 MPIU_printf(snes->comm, "iter = %d, Function norm %g \n",its,fnorm); 128 return 0; 129 } 130 131 int SNESDefaultSMonitor(SNES snes,int its, double fnorm,void *dummy) 132 { 133 if (fnorm > 1.e-9 || fnorm == 0.0) { 134 MPIU_printf(snes->comm, "iter = %d, Function norm %g \n",its,fnorm); 135 } 136 else if (fnorm > 1.e-11){ 137 MPIU_printf(snes->comm, "iter = %d, Function norm %5.3e \n",its,fnorm); 138 } 139 else { 140 MPIU_printf(snes->comm, "iter = %d, Function norm < 1.e-11\n",its); 141 } 142 return 0; 143 } 144 145 /*@ 146 SNESDefaultConverged - Default test for monitoring the convergence 147 of the SNES solvers. 148 149 Input Parameters: 150 . snes - the SNES context 151 . xnorm - 2-norm of current iterate 152 . pnorm - 2-norm of current step 153 . fnorm - 2-norm of residual 154 . dummy - unused context 155 156 Returns: 157 $ 2 if ( fnorm < atol ), 158 $ 3 if ( pnorm < xtol*xnorm ), 159 $ -2 if ( nres > max_res ), 160 $ 0 otherwise, 161 162 where 163 $ atol - absolute residual norm tolerance, 164 $ set with SNESSetAbsoluteTolerance() 165 $ max_res - maximum number of residual evaluations, 166 $ set with SNESSetMaxResidualEvaluations() 167 $ nres - number of residual evaluations 168 $ xtol - relative residual norm tolerance, 169 $ set with SNESSetRelativeTolerance() 170 171 .keywords: SNES, nonlinear, default, converged, convergence 172 173 .seealso: SNESSetConvergenceTest() 174 @*/ 175 int SNESDefaultConverged(SNES snes,double xnorm,double pnorm,double fnorm, 176 void *dummy) 177 { 178 if (fnorm < snes->atol) { 179 PLogInfo((PetscObject)snes, 180 "Converged due to absolute residual norm %g < %g\n",fnorm,snes->atol); 181 return 2; 182 } 183 if (pnorm < snes->xtol*(xnorm)) { 184 PLogInfo((PetscObject)snes, 185 "Converged due to small update length %g < %g*%g\n", 186 pnorm,snes->xtol,xnorm); 187 return 3; 188 } 189 if (snes->nfuncs > snes->max_funcs) { 190 PLogInfo((PetscObject)snes, 191 "Exceeded maximum number of residual evaluations: %d > %d\n", 192 snes->nfuncs, snes->max_funcs ); 193 return -2; 194 } 195 return 0; 196 } 197 198 /* ------------------------------------------------------------ */ 199 /*ARGSUSED*/ 200 /*@ 201 SNESNoLineSearch - This routine is not a line search at all; 202 it simply uses the full Newton step. Thus, this routine is intended 203 to serve as a template and is not recommended for general use. 204 205 Input Parameters: 206 . snes - nonlinear context 207 . x - current iterate 208 . f - residual evaluated at x 209 . y - search direction (contains new iterate on output) 210 . w - work vector 211 . fnorm - 2-norm of f 212 213 Output Parameters: 214 . g - residual evaluated at new iterate y 215 . y - new iterate (contains search direction on input) 216 . gnorm - 2-norm of g 217 . ynorm - 2-norm of search length 218 219 Options Database Key: 220 $ -snes_line_search basic 221 222 Returns: 223 1, indicating success of the step. 224 225 .keywords: SNES, nonlinear, line search, cubic 226 227 .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), 228 .seealso: SNESSetLineSearchRoutine() 229 @*/ 230 int SNESNoLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 231 double fnorm, double *ynorm, double *gnorm ) 232 { 233 int ierr; 234 Scalar one = 1.0; 235 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 236 VecNorm(y, ynorm ); /* ynorm = || y || */ 237 VecAXPY(&one, x, y ); /* y <- x + y */ 238 ierr = SNESComputeFunction(snes,y,g); CHKERRQ(ierr); 239 VecNorm( g, gnorm ); /* gnorm = || g || */ 240 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 241 return 1; 242 } 243 /* ------------------------------------------------------------------ */ 244 /*@ 245 SNESCubicLineSearch - This routine performs a cubic line search and 246 is the default line search method. 247 248 Input Parameters: 249 . snes - nonlinear context 250 . x - current iterate 251 . f - residual evaluated at x 252 . y - search direction (contains new iterate on output) 253 . w - work vector 254 . fnorm - 2-norm of f 255 256 Output parameters: 257 . g - residual evaluated at new iterate y 258 . y - new iterate (contains search direction on input) 259 . gnorm - 2-norm of g 260 . ynorm - 2-norm of search length 261 262 Returns: 263 1 if the line search succeeds; 0 if the line search fails. 264 265 Options Database Key: 266 $ -snes_line_search cubic 267 268 Notes: 269 This line search is taken from "Numerical Methods for Unconstrained 270 Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 271 272 .keywords: SNES, nonlinear, line search, cubic 273 274 .seealso: SNESNoLineSearch(), SNESNoLineSearch(), SNESSetLineSearchRoutine() 275 @*/ 276 int SNESCubicLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 277 double fnorm, double *ynorm, double *gnorm ) 278 { 279 double steptol, initslope; 280 double lambdaprev, gnormprev; 281 double a, b, d, t1, t2; 282 #if defined(PETSC_COMPLEX) 283 Scalar cinitslope,clambda; 284 #endif 285 int ierr,count; 286 SNES_LS *neP = (SNES_LS *) snes->data; 287 Scalar one = 1.0,scale; 288 double maxstep,minlambda,alpha,lambda,lambdatemp; 289 290 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 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 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 324 return 0; 325 } 326 327 /* Fit points with quadratic */ 328 lambda = 1.0; count = 0; 329 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 330 lambdaprev = lambda; 331 gnormprev = *gnorm; 332 if (lambdatemp <= .1*lambda) { 333 lambda = .1*lambda; 334 } else lambda = lambdatemp; 335 VecCopy(x, w ); 336 #if defined(PETSC_COMPLEX) 337 clambda = lambda; VecAXPY(&clambda, y, w ); 338 #else 339 VecAXPY(&lambda, y, w ); 340 #endif 341 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 342 VecNorm(g, gnorm ); 343 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 344 VecCopy(w, y ); 345 PLogInfo((PetscObject)snes,"Quadratically determined step, lambda %g\n",lambda); 346 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 347 return 0; 348 } 349 350 /* Fit points with cubic */ 351 count = 1; 352 while (1) { 353 if (lambda <= minlambda) { /* bad luck; use full step */ 354 PLogInfo((PetscObject)snes,"Unable to find good step length! %d \n",count); 355 PLogInfo((PetscObject)snes, "f %g fnew %g ynorm %g lambda %g \n", 356 fnorm,*gnorm, *ynorm,lambda); 357 VecCopy(w, y ); 358 break; 359 } 360 t1 = *gnorm - fnorm - lambda*initslope; 361 t2 = gnormprev - fnorm - lambdaprev*initslope; 362 a = (t1/(lambda*lambda) - 363 t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 364 b = (-lambdaprev*t1/(lambda*lambda) + 365 lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 366 d = b*b - 3*a*initslope; 367 if (d < 0.0) d = 0.0; 368 if (a == 0.0) { 369 lambdatemp = -initslope/(2.0*b); 370 } else { 371 lambdatemp = (-b + sqrt(d))/(3.0*a); 372 } 373 if (lambdatemp > .5*lambda) { 374 lambdatemp = .5*lambda; 375 } 376 lambdaprev = lambda; 377 gnormprev = *gnorm; 378 if (lambdatemp <= .1*lambda) { 379 lambda = .1*lambda; 380 } 381 else lambda = lambdatemp; 382 VecCopy( x, w ); 383 #if defined(PETSC_COMPLEX) 384 VecAXPY(&clambda, y, w ); 385 #else 386 VecAXPY(&lambda, y, w ); 387 #endif 388 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 389 VecNorm(g, gnorm ); 390 if (*gnorm <= fnorm + alpha*initslope) { /* is reduction enough */ 391 VecCopy(w, y ); 392 PLogInfo((PetscObject)snes,"Cubically determined step, lambda %g\n",lambda); 393 break; 394 } 395 count++; 396 } 397 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 398 return 0; 399 } 400 /*@ 401 SNESQuadraticLineSearch - This routine performs a cubic line search. 402 403 Input Parameters: 404 . snes - the SNES context 405 . x - current iterate 406 . f - residual evaluated at x 407 . y - search direction (contains new iterate on output) 408 . w - work vector 409 . fnorm - 2-norm of f 410 411 Output Parameters: 412 . g - residual evaluated at new iterate y 413 . y - new iterate (contains search direction on input) 414 . gnorm - 2-norm of g 415 . ynorm - 2-norm of search length 416 417 Returns: 418 1 if the line search succeeds; 0 if the line search fails. 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 ) 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 alpha = neP->alpha; 446 maxstep = neP->maxstep; 447 steptol = neP->steptol; 448 449 VecNorm(y, ynorm ); 450 if (*ynorm > maxstep) { /* Step too big, so scale back */ 451 scale = maxstep/(*ynorm); 452 VecScale(&scale, y ); 453 *ynorm = maxstep; 454 } 455 minlambda = steptol/(*ynorm); 456 #if defined(PETSC_COMPLEX) 457 VecDot(f, y, &cinitslope ); 458 initslope = real(cinitslope); 459 #else 460 VecDot(f, y, &initslope ); 461 #endif 462 if (initslope > 0.0) initslope = -initslope; 463 if (initslope == 0.0) initslope = -1.0; 464 465 VecCopy(y, w ); 466 VecAXPY(&one, x, w ); 467 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 468 VecNorm(g, gnorm ); 469 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 470 VecCopy(w, y ); 471 PLogInfo((PetscObject)snes,"Using full step\n"); 472 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 473 return 0; 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 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 509 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 510 return 0; 511 } 512 /* ------------------------------------------------------------ */ 513 /*@C 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, double *gnorm) 532 533 Input parameters for func: 534 . snes - nonlinear context 535 . x - current iterate 536 . f - residual evaluated at x 537 . y - search direction (contains new iterate on output) 538 . w - work vector 539 . fnorm - 2-norm of f 540 541 Output parameters for func: 542 . g - residual evaluated at new iterate y 543 . y - new iterate (contains search direction on input) 544 . gnorm - 2-norm of g 545 . ynorm - 2-norm of search length 546 547 Returned by func: 548 1 if the line search succeeds; 0 if the line search fails. 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*) ) 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 SNESSetFromOptions_LS(SNES snes) 573 { 574 SNES_LS *ls = (SNES_LS *)snes->data; 575 char ver[16]; 576 double tmp; 577 578 if (OptionsGetDouble(snes->prefix,"-snes_line_search_alpa",&tmp)) { 579 ls->alpha = tmp; 580 } 581 if (OptionsGetDouble(snes->prefix,"-snes_line_search_maxstep",&tmp)) { 582 ls->maxstep = tmp; 583 } 584 if (OptionsGetDouble(snes->prefix,"-snes_line_search_steptol",&tmp)) { 585 ls->steptol = tmp; 586 } 587 if (OptionsGetString(snes->prefix,"-snes_line_search",ver,16)) { 588 if (!strcmp(ver,"basic")) { 589 SNESSetLineSearchRoutine(snes,SNESNoLineSearch); 590 } 591 else if (!strcmp(ver,"quadratic")) { 592 SNESSetLineSearchRoutine(snes,SNESQuadraticLineSearch); 593 } 594 else if (!strcmp(ver,"cubic")) { 595 SNESSetLineSearchRoutine(snes,SNESCubicLineSearch); 596 } 597 else {SETERRQ(1,"Unknown line search?");} 598 } 599 return 0; 600 } 601 602 /* ------------------------------------------------------------ */ 603 int SNESCreate_LS(SNES snes ) 604 { 605 SNES_LS *neP; 606 607 snes->type = SNES_NLS; 608 snes->setup = SNESSetUp_LS; 609 snes->solve = SNESSolve_LS; 610 snes->destroy = SNESDestroy_LS; 611 snes->Converged = SNESDefaultConverged; 612 snes->printhelp = SNESPrintHelp_LS; 613 snes->setfromoptions = SNESSetFromOptions_LS; 614 615 neP = PETSCNEW(SNES_LS); CHKPTRQ(neP); 616 snes->data = (void *) neP; 617 neP->alpha = 1.e-4; 618 neP->maxstep = 1.e8; 619 neP->steptol = 1.e-12; 620 neP->LineSearch = SNESCubicLineSearch; 621 return 0; 622 } 623 624 625 626 627