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