1 #ifndef lint 2 static char vcid[] = "$Id: ls.c,v 1.20 1995/05/30 22:52:06 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; 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); CHKERR(ierr); /* X <- X_0 */ 47 VecNorm( X, &xnorm ); /* xnorm = || X || */ 48 ierr = SNESComputeFunction(snes,X,F); CHKERR(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); CHKERR(ierr); 60 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg); 61 ierr = SLESSolve(snes->sles,F,Y,&lits); CHKERR(ierr); 62 ierr = (*neP->LineSearch)(snes, X, F, G, Y, W, fnorm, &ynorm, &gnorm ); 63 CHKERR(ierr); 64 65 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 66 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 67 fnorm = gnorm; 68 69 snes->norm = fnorm; 70 if (history && history_len > i+1) history[i+1] = fnorm; 71 VecNorm( X, &xnorm ); /* xnorm = || X || */ 72 if (snes->Monitor) (*snes->Monitor)(snes,i+1,fnorm,snes->monP); 73 74 /* Test for convergence */ 75 if ((*snes->Converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 76 if (X != snes->vec_sol) { 77 VecCopy( X, snes->vec_sol ); 78 snes->vec_sol_always = snes->vec_sol; 79 snes->vec_func_always = snes->vec_func; 80 } 81 break; 82 } 83 } 84 if (i == maxits) i--; 85 *outits = i+1; 86 return 0; 87 } 88 /* ------------------------------------------------------------ */ 89 int SNESSetUp_LS(SNES snes ) 90 { 91 int ierr; 92 snes->nwork = 3; 93 ierr = VecGetVecs( snes->vec_sol, snes->nwork,&snes->work ); CHKERR(ierr); 94 PLogObjectParents(snes,snes->nwork,snes->work ); 95 return 0; 96 } 97 /* ------------------------------------------------------------ */ 98 int SNESDestroy_LS(PetscObject obj) 99 { 100 SNES snes = (SNES) obj; 101 VecFreeVecs(snes->work, snes->nwork ); 102 FREE(snes->data); 103 return 0; 104 } 105 /*@ 106 SNESDefaultMonitor - Default SNES monitoring routine. Prints the 107 residual norm at each iteration. 108 109 Input Parameters: 110 . snes - the SNES context 111 . its - iteration number 112 . fnorm - 2-norm residual value (may be estimated) 113 . dummy - unused context 114 115 Notes: 116 f is either the residual or its negative, depending on the user's 117 preference, as set with SNESSetFunction(). 118 119 .keywords: SNES, nonlinear, default, monitor, residual, norm 120 121 .seealso: SNESSetMonitor() 122 @*/ 123 int SNESDefaultMonitor(SNES snes,int its, double fnorm,void *dummy) 124 { 125 MPIU_printf(snes->comm, "iter = %d, Function norm %g \n",its,fnorm); 126 return 0; 127 } 128 129 int SNESDefaultSMonitor(SNES snes,int its, double fnorm,void *dummy) 130 { 131 if (fnorm > 1.e-9 || fnorm == 0.0) { 132 MPIU_printf(snes->comm, "iter = %d, Function norm %g \n",its,fnorm); 133 } 134 else if (fnorm > 1.e-11){ 135 MPIU_printf(snes->comm, "iter = %d, Function norm %5.3e \n",its,fnorm); 136 } 137 else { 138 MPIU_printf(snes->comm, "iter = %d, Function norm < 1.e-11\n",its); 139 } 140 return 0; 141 } 142 143 /*@ 144 SNESDefaultConverged - Default test for monitoring the convergence 145 of the SNES solvers. 146 147 Input Parameters: 148 . snes - the SNES context 149 . xnorm - 2-norm of current iterate 150 . pnorm - 2-norm of current step 151 . fnorm - 2-norm of residual 152 . dummy - unused context 153 154 Returns: 155 $ 2 if ( fnorm < atol ), 156 $ 3 if ( pnorm < xtol*xnorm ), 157 $ -2 if ( nres > max_res ), 158 $ 0 otherwise, 159 160 where 161 $ atol - absolute residual norm tolerance, 162 $ set with SNESSetAbsoluteTolerance() 163 $ max_res - maximum number of residual evaluations, 164 $ set with SNESSetMaxResidualEvaluations() 165 $ nres - number of residual evaluations 166 $ xtol - relative residual norm tolerance, 167 $ set with SNESSetRelativeTolerance() 168 169 .keywords: SNES, nonlinear, default, converged, convergence 170 171 .seealso: SNESSetConvergenceTest() 172 @*/ 173 int SNESDefaultConverged(SNES snes,double xnorm,double pnorm,double fnorm, 174 void *dummy) 175 { 176 if (fnorm < snes->atol) { 177 PLogInfo((PetscObject)snes, 178 "Converged due to absolute residual norm %g < %g\n",fnorm,snes->atol); 179 return 2; 180 } 181 if (pnorm < snes->xtol*(xnorm)) { 182 PLogInfo((PetscObject)snes, 183 "Converged due to small update length %g < %g*%g\n", 184 pnorm,snes->xtol,xnorm); 185 return 3; 186 } 187 if (snes->nfuncs > snes->max_funcs) { 188 PLogInfo((PetscObject)snes, 189 "Exceeded maximum number of residual evaluations: %d > %d\n", 190 snes->nfuncs, snes->max_funcs ); 191 return -2; 192 } 193 return 0; 194 } 195 196 /* ------------------------------------------------------------ */ 197 /*ARGSUSED*/ 198 /*@ 199 SNESNoLineSearch - This routine is not a line search at all; 200 it simply uses the full Newton step. Thus, this routine is intended 201 to serve as a template and is not recommended for general use. 202 203 Input Parameters: 204 . snes - nonlinear context 205 . x - current iterate 206 . f - residual evaluated at x 207 . y - search direction (contains new iterate on output) 208 . w - work vector 209 . fnorm - 2-norm of f 210 211 Output Parameters: 212 . g - residual evaluated at new iterate y 213 . y - new iterate (contains search direction on input) 214 . gnorm - 2-norm of g 215 . ynorm - 2-norm of search length 216 217 Options Database Key: 218 $ -snes_line_search basic 219 220 Returns: 221 1, indicating success of the step. 222 223 .keywords: SNES, nonlinear, line search, cubic 224 225 .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), 226 .seealso: SNESSetLineSearchRoutine() 227 @*/ 228 int SNESNoLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 229 double fnorm, double *ynorm, double *gnorm ) 230 { 231 int ierr; 232 Scalar one = 1.0; 233 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 234 VecNorm(y, ynorm ); /* ynorm = || y || */ 235 VecAXPY(&one, x, y ); /* y <- x + y */ 236 ierr = SNESComputeFunction(snes,y,g); CHKERR(ierr); 237 VecNorm( g, gnorm ); /* gnorm = || g || */ 238 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 239 return 1; 240 } 241 /* ------------------------------------------------------------------ */ 242 /*@ 243 SNESCubicLineSearch - This routine performs a cubic line search and 244 is the default line search method. 245 246 Input Parameters: 247 . snes - nonlinear context 248 . x - current iterate 249 . f - residual evaluated at x 250 . y - search direction (contains new iterate on output) 251 . w - work vector 252 . fnorm - 2-norm of f 253 254 Output parameters: 255 . g - residual evaluated at new iterate y 256 . y - new iterate (contains search direction on input) 257 . gnorm - 2-norm of g 258 . ynorm - 2-norm of search length 259 260 Returns: 261 1 if the line search succeeds; 0 if the line search fails. 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 ) 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 alpha = neP->alpha; 290 maxstep = neP->maxstep; 291 steptol = neP->steptol; 292 293 VecNorm(y, ynorm ); 294 if (*ynorm > maxstep) { /* Step too big, so scale back */ 295 scale = maxstep/(*ynorm); 296 #if defined(PETSC_COMPLEX) 297 PLogInfo((PetscObject)snes,"Scaling step by %g\n",real(scale)); 298 #else 299 PLogInfo((PetscObject)snes,"Scaling step by %g\n",scale); 300 #endif 301 VecScale(&scale, y ); 302 *ynorm = maxstep; 303 } 304 minlambda = steptol/(*ynorm); 305 #if defined(PETSC_COMPLEX) 306 VecDot(f, y, &cinitslope ); 307 initslope = real(cinitslope); 308 #else 309 VecDot(f, y, &initslope ); 310 #endif 311 if (initslope > 0.0) initslope = -initslope; 312 if (initslope == 0.0) initslope = -1.0; 313 314 VecCopy(y, w ); 315 VecAXPY(&one, x, w ); 316 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 317 VecNorm(g, gnorm ); 318 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 319 VecCopy(w, y ); 320 PLogInfo((PetscObject)snes,"Using full step\n"); 321 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 322 return 0; 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); CHKERR(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 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 345 return 0; 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 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 357 return -1; 358 } 359 t1 = *gnorm - fnorm - lambda*initslope; 360 t2 = gnormprev - fnorm - lambdaprev*initslope; 361 a = (t1/(lambda*lambda) - 362 t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 363 b = (-lambdaprev*t1/(lambda*lambda) + 364 lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 365 d = b*b - 3*a*initslope; 366 if (d < 0.0) d = 0.0; 367 if (a == 0.0) { 368 lambdatemp = -initslope/(2.0*b); 369 } else { 370 lambdatemp = (-b + sqrt(d))/(3.0*a); 371 } 372 if (lambdatemp > .5*lambda) { 373 lambdatemp = .5*lambda; 374 } 375 lambdaprev = lambda; 376 gnormprev = *gnorm; 377 if (lambdatemp <= .1*lambda) { 378 lambda = .1*lambda; 379 } 380 else lambda = lambdatemp; 381 VecCopy( x, w ); 382 #if defined(PETSC_COMPLEX) 383 VecAXPY(&clambda, y, w ); 384 #else 385 VecAXPY(&lambda, y, w ); 386 #endif 387 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 388 VecNorm(g, gnorm ); 389 if (*gnorm <= fnorm + alpha*initslope) { /* is reduction enough */ 390 VecCopy(w, y ); 391 PLogInfo((PetscObject)snes,"Cubically determined step, lambda %g\n",lambda); 392 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 393 return 0; 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); CHKERR(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 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 486 return 0; 487 } 488 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 489 lambdaprev = lambda; 490 gnormprev = *gnorm; 491 if (lambdatemp <= .1*lambda) { 492 lambda = .1*lambda; 493 } else lambda = lambdatemp; 494 VecCopy(x, w ); 495 #if defined(PETSC_COMPLEX) 496 clambda = lambda; VecAXPY(&clambda, y, w ); 497 #else 498 VecAXPY(&lambda, y, w ); 499 #endif 500 ierr = SNESComputeFunction(snes,w,g); CHKERR(ierr); 501 VecNorm(g, gnorm ); 502 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 503 VecCopy(w, y ); 504 PLogInfo((PetscObject)snes,"Quadratically determined step, lambda %g\n",lambda); 505 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 506 return 0; 507 } 508 count++; 509 } 510 511 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 512 return 0; 513 } 514 /* ------------------------------------------------------------ */ 515 /*@C 516 SNESSetLineSearchRoutine - Sets the line search routine to be used 517 by the method SNES_LS. 518 519 Input Parameters: 520 . snes - nonlinear context obtained from SNESCreate() 521 . func - pointer to int function 522 523 Available Routines: 524 . SNESCubicLineSearch() - default line search 525 . SNESQuadraticLineSearch() - quadratic line search 526 . SNESNoLineSearch() - the full Newton step (actually not a line search) 527 528 Options Database Keys: 529 $ -snes_line_search [basic,quadratic,cubic] 530 531 Calling sequence of func: 532 func (SNES snes, Vec x, Vec f, Vec g, Vec y, 533 Vec w, double fnorm, double *ynorm, double *gnorm) 534 535 Input parameters for func: 536 . snes - nonlinear context 537 . x - current iterate 538 . f - residual evaluated at x 539 . y - search direction (contains new iterate on output) 540 . w - work vector 541 . fnorm - 2-norm of f 542 543 Output parameters for func: 544 . g - residual evaluated at new iterate y 545 . y - new iterate (contains search direction on input) 546 . gnorm - 2-norm of g 547 . ynorm - 2-norm of search length 548 549 Returned by func: 550 1 if the line search succeeds; 0 if the line search fails. 551 552 .keywords: SNES, nonlinear, set, line search, routine 553 554 .seealso: SNESNoLineSearch(), SNESQuadraticLineSearch(), SNESCubicLineSearch() 555 @*/ 556 int SNESSetLineSearchRoutine(SNES snes,int (*func)(SNES,Vec,Vec,Vec,Vec,Vec, 557 double,double *,double*) ) 558 { 559 if ((snes)->type == SNES_NLS) 560 ((SNES_LS *)(snes->data))->LineSearch = func; 561 return 0; 562 } 563 564 static int SNESPrintHelp_LS(SNES snes) 565 { 566 SNES_LS *ls = (SNES_LS *)snes->data; 567 fprintf(stderr,"-snes_line_search [basic,quadratic,cubic]\n"); 568 fprintf(stderr,"-snes_line_search_alpha alpha (default %g)\n",ls->alpha); 569 fprintf(stderr,"-snes_line_search_maxstep max (default %g)\n",ls->maxstep); 570 fprintf(stderr,"-snes_line_search_steptol tol (default %g)\n",ls->steptol); 571 return 0; 572 } 573 574 static int SNESSetFromOptions_LS(SNES snes) 575 { 576 SNES_LS *ls = (SNES_LS *)snes->data; 577 char ver[16]; 578 double tmp; 579 580 if (OptionsGetDouble(snes->prefix,"-snes_line_search_alpa",&tmp)) { 581 ls->alpha = tmp; 582 } 583 if (OptionsGetDouble(snes->prefix,"-snes_line_search_maxstep",&tmp)) { 584 ls->maxstep = tmp; 585 } 586 if (OptionsGetDouble(snes->prefix,"-snes_line_search_steptol",&tmp)) { 587 ls->steptol = tmp; 588 } 589 if (OptionsGetString(snes->prefix,"-snes_line_search",ver,16)) { 590 if (!strcmp(ver,"basic")) { 591 SNESSetLineSearchRoutine(snes,SNESNoLineSearch); 592 } 593 else if (!strcmp(ver,"quadratic")) { 594 SNESSetLineSearchRoutine(snes,SNESQuadraticLineSearch); 595 } 596 else if (!strcmp(ver,"cubic")) { 597 SNESSetLineSearchRoutine(snes,SNESCubicLineSearch); 598 } 599 else {SETERR(1,"Unknown line search?");} 600 } 601 return 0; 602 } 603 604 /* ------------------------------------------------------------ */ 605 int SNESCreate_LS(SNES snes ) 606 { 607 SNES_LS *neP; 608 609 snes->type = SNES_NLS; 610 snes->setup = SNESSetUp_LS; 611 snes->solve = SNESSolve_LS; 612 snes->destroy = SNESDestroy_LS; 613 snes->Converged = SNESDefaultConverged; 614 snes->printhelp = SNESPrintHelp_LS; 615 snes->setfromoptions = SNESSetFromOptions_LS; 616 617 neP = NEW(SNES_LS); CHKPTR(neP); 618 snes->data = (void *) neP; 619 neP->alpha = 1.e-4; 620 neP->maxstep = 1.e8; 621 neP->steptol = 1.e-12; 622 neP->LineSearch = SNESCubicLineSearch; 623 return 0; 624 } 625 626 627 628 629