1 #ifndef lint 2 static char vcid[] = "$Id: ls.c,v 1.55 1996/01/12 22:09:42 bsmith Exp balay $"; 3 #endif 4 5 #include <math.h> 6 #include "ls.h" 7 #include "pinclude/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. 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 ierr = VecNorm(X,NORM_2,&xnorm); CHKERRQ(ierr); /* xnorm = || X || */ 46 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* F(X) */ 47 ierr = VecNorm(F,NORM_2,&fnorm); CHKERRQ(ierr); /* fnorm <- ||F|| */ 48 snes->norm = fnorm; 49 if (history && history_len > 0) history[0] = fnorm; 50 if (snes->monitor) {ierr = (*snes->monitor)(snes,0,fnorm,snes->monP); CHKERRQ(ierr);} 51 52 for ( i=0; i<maxits; i++ ) { 53 snes->iter = i+1; 54 55 /* Solve J Y = -F, where J is Jacobian matrix */ 56 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 57 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg); CHKERRQ(ierr); 58 ierr = SLESSolve(snes->sles,F,Y,&lits); CHKERRQ(ierr); 59 ierr = VecCopy(Y,snes->vec_sol_update_always); CHKERRQ(ierr); 60 ierr = (*neP->LineSearch)(snes,X,F,G,Y,W,fnorm,&ynorm,&gnorm,&lsfail);CHKERRQ(ierr); 61 if (lsfail) snes->nfailures++; 62 63 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 64 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 65 fnorm = gnorm; 66 67 snes->norm = fnorm; 68 if (history && history_len > i+1) history[i+1] = fnorm; 69 ierr = VecNorm(X,NORM_2,&xnorm); CHKERRQ(ierr); /* xnorm = || X || */ 70 if (snes->monitor) {ierr = (*snes->monitor)(snes,i+1,fnorm,snes->monP);CHKERRQ(ierr);} 71 72 /* Test for convergence */ 73 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 74 if (X != snes->vec_sol) { 75 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 76 snes->vec_sol_always = snes->vec_sol; 77 snes->vec_func_always = snes->vec_func; 78 } 79 break; 80 } 81 } 82 if (i == maxits) { 83 PLogInfo((PetscObject)snes, 84 "SNESSolve_LS: Maximum number of iterations has been reached: %d\n",maxits); 85 i--; 86 } 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 = VecDuplicateVecs(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 int ierr; 105 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 106 PetscFree(snes->data); 107 return 0; 108 } 109 /* ------------------------------------------------------------ */ 110 /*ARGSUSED*/ 111 /*@C 112 SNESNoLineSearch - This routine is not a line search at all; 113 it simply uses the full Newton step. Thus, this routine is intended 114 to serve as a template and is not recommended for general use. 115 116 Input Parameters: 117 . snes - nonlinear context 118 . x - current iterate 119 . f - residual evaluated at x 120 . y - search direction (contains new iterate on output) 121 . w - work vector 122 . fnorm - 2-norm of f 123 124 Output Parameters: 125 . g - residual evaluated at new iterate y 126 . y - new iterate (contains search direction on input) 127 . gnorm - 2-norm of g 128 . ynorm - 2-norm of search length 129 . flag - set to 0, indicating a successful line search 130 131 Options Database Key: 132 $ -snes_line_search basic 133 134 .keywords: SNES, nonlinear, line search, cubic 135 136 .seealso: SNESCubicLineSearch(), SNESQuadraticLineSearch(), 137 SNESSetLineSearchRoutine() 138 @*/ 139 int SNESNoLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 140 double fnorm, double *ynorm, double *gnorm,int *flag ) 141 { 142 int ierr; 143 Scalar mone = -1.0; 144 145 *flag = 0; 146 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 147 ierr = VecNorm(y,NORM_2,ynorm); CHKERRQ(ierr); /* ynorm = || y || */ 148 ierr = VecAYPX(&mone,x,y); CHKERRQ(ierr); /* y <- x + y */ 149 ierr = SNESComputeFunction(snes,y,g); CHKERRQ(ierr); /* F(y) */ 150 ierr = VecNorm(g,NORM_2,gnorm); CHKERRQ(ierr); /* gnorm = || g || */ 151 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 152 return 0; 153 } 154 /* ------------------------------------------------------------------ */ 155 /*@C 156 SNESCubicLineSearch - This routine performs a cubic line search and 157 is the default line search method. 158 159 Input Parameters: 160 . snes - nonlinear context 161 . x - current iterate 162 . f - residual evaluated at x 163 . y - search direction (contains new iterate on output) 164 . w - work vector 165 . fnorm - 2-norm of f 166 167 Output Parameters: 168 . g - residual evaluated at new iterate y 169 . y - new iterate (contains search direction on input) 170 . gnorm - 2-norm of g 171 . ynorm - 2-norm of search length 172 . flag - 0 if line search succeeds; -1 on failure. 173 174 Options Database Key: 175 $ -snes_line_search cubic 176 177 Notes: 178 This line search is taken from "Numerical Methods for Unconstrained 179 Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 180 181 .keywords: SNES, nonlinear, line search, cubic 182 183 .seealso: SNESNoLineSearch(), SNESNoLineSearch(), SNESSetLineSearchRoutine() 184 @*/ 185 int SNESCubicLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 186 double fnorm, double *ynorm, double *gnorm,int *flag) 187 { 188 double steptol, initslope,lambdaprev, gnormprev,a, b, d, t1, t2; 189 double maxstep,minlambda,alpha,lambda,lambdatemp; 190 #if defined(PETSC_COMPLEX) 191 Scalar cinitslope,clambda; 192 #endif 193 int ierr, count; 194 SNES_LS *neP = (SNES_LS *) snes->data; 195 Scalar mone = -1.0,scale; 196 197 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 198 *flag = 0; 199 alpha = neP->alpha; 200 maxstep = neP->maxstep; 201 steptol = neP->steptol; 202 203 ierr = VecNorm(y,NORM_2,ynorm); CHKERRQ(ierr); 204 if (*ynorm > maxstep) { /* Step too big, so scale back */ 205 scale = maxstep/(*ynorm); 206 #if defined(PETSC_COMPLEX) 207 PLogInfo((PetscObject)snes,"SNESCubicLineSearch: Scaling step by %g\n",real(scale)); 208 #else 209 PLogInfo((PetscObject)snes,"SNESCubicLineSearch: Scaling step by %g\n",scale); 210 #endif 211 ierr = VecScale(&scale,y); CHKERRQ(ierr); 212 *ynorm = maxstep; 213 } 214 minlambda = steptol/(*ynorm); 215 ierr = MatMult(snes->jacobian,y,w); CHKERRQ(ierr); 216 #if defined(PETSC_COMPLEX) 217 ierr = VecDot(f,w,&cinitslope); CHKERRQ(ierr); 218 initslope = real(cinitslope); 219 #else 220 ierr = VecDot(f,w,&initslope); CHKERRQ(ierr); 221 #endif 222 if (initslope > 0.0) initslope = -initslope; 223 if (initslope == 0.0) initslope = -1.0; 224 225 ierr = VecCopy(y,w); CHKERRQ(ierr); 226 ierr = VecAYPX(&mone,x,w); CHKERRQ(ierr); 227 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 228 ierr = VecNorm(g,NORM_2,gnorm); 229 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 230 ierr = VecCopy(w,y); CHKERRQ(ierr); 231 PLogInfo((PetscObject)snes,"SNESCubicLineSearch: Using full step\n"); 232 goto theend; 233 } 234 235 /* Fit points with quadratic */ 236 lambda = 1.0; count = 0; 237 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 238 lambdaprev = lambda; 239 gnormprev = *gnorm; 240 if (lambdatemp <= .1*lambda) lambda = .1*lambda; 241 else lambda = lambdatemp; 242 ierr = VecCopy(x,w); CHKERRQ(ierr); 243 lambda = -lambda; 244 #if defined(PETSC_COMPLEX) 245 clambda = lambda; ierr = VecAXPY(&clambda,y,w); CHKERRQ(ierr); 246 #else 247 ierr = VecAXPY(&lambda,y,w); CHKERRQ(ierr); 248 #endif 249 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 250 ierr = VecNorm(g,NORM_2,gnorm); CHKERRQ(ierr); 251 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 252 ierr = VecCopy(w,y); CHKERRQ(ierr); 253 PLogInfo((PetscObject)snes, 254 "SNESCubicLineSearch: Quadratically determined step, lambda %g\n",lambda); 255 goto theend; 256 } 257 258 /* Fit points with cubic */ 259 count = 1; 260 while (1) { 261 if (lambda <= minlambda) { /* bad luck; use full step */ 262 PLogInfo((PetscObject)snes, 263 "SNESCubicLineSearch:Unable to find good step length! %d \n",count); 264 PLogInfo((PetscObject)snes, 265 "SNESCubicLineSearch:f %g fnew %g ynorm %g lambda %g\n",fnorm,*gnorm,*ynorm,lambda); 266 ierr = VecCopy(w,y); CHKERRQ(ierr); 267 *flag = -1; break; 268 } 269 t1 = *gnorm - fnorm - lambda*initslope; 270 t2 = gnormprev - fnorm - lambdaprev*initslope; 271 a = (t1/(lambda*lambda) - t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 272 b = (-lambdaprev*t1/(lambda*lambda) + 273 lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 274 d = b*b - 3*a*initslope; 275 if (d < 0.0) d = 0.0; 276 if (a == 0.0) { 277 lambdatemp = -initslope/(2.0*b); 278 } else { 279 lambdatemp = (-b + sqrt(d))/(3.0*a); 280 } 281 if (lambdatemp > .5*lambda) { 282 lambdatemp = .5*lambda; 283 } 284 lambdaprev = lambda; 285 gnormprev = *gnorm; 286 if (lambdatemp <= .1*lambda) { 287 lambda = .1*lambda; 288 } 289 else lambda = lambdatemp; 290 ierr = VecCopy( x, w ); CHKERRQ(ierr); 291 lambda = -lambda; 292 #if defined(PETSC_COMPLEX) 293 clambda = lambda; 294 ierr = VecAXPY(&clambda,y,w); CHKERRQ(ierr); 295 #else 296 ierr = VecAXPY(&lambda,y,w); CHKERRQ(ierr); 297 #endif 298 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 299 ierr = VecNorm(g,NORM_2,gnorm); CHKERRQ(ierr); 300 if (*gnorm <= fnorm + alpha*initslope) { /* is reduction enough */ 301 ierr = VecCopy(w,y); CHKERRQ(ierr); 302 PLogInfo((PetscObject)snes, 303 "SNESCubicLineSearch: Cubically determined step, lambda %g\n",lambda); 304 *flag = -1; break; 305 } 306 count++; 307 } 308 theend: 309 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 310 return 0; 311 } 312 /* ------------------------------------------------------------------ */ 313 /*@C 314 SNESQuadraticLineSearch - This routine performs a cubic line search. 315 316 Input Parameters: 317 . snes - the SNES context 318 . x - current iterate 319 . f - residual evaluated at x 320 . y - search direction (contains new iterate on output) 321 . w - work vector 322 . fnorm - 2-norm of f 323 324 Output Parameters: 325 . g - residual evaluated at new iterate y 326 . y - new iterate (contains search direction on input) 327 . gnorm - 2-norm of g 328 . ynorm - 2-norm of search length 329 . flag - 0 if line search succeeds; -1 on failure. 330 331 Options Database Key: 332 $ -snes_line_search quadratic 333 334 Notes: 335 Use SNESSetLineSearchRoutine() 336 to set this routine within the SNES_EQ_NLS method. 337 338 .keywords: SNES, nonlinear, quadratic, line search 339 340 .seealso: SNESCubicLineSearch(), SNESNoLineSearch(), SNESSetLineSearchRoutine() 341 @*/ 342 int SNESQuadraticLineSearch(SNES snes, Vec x, Vec f, Vec g, Vec y, Vec w, 343 double fnorm, double *ynorm, double *gnorm,int *flag) 344 { 345 double steptol,initslope,lambdaprev,gnormprev,maxstep,minlambda,alpha,lambda,lambdatemp; 346 #if defined(PETSC_COMPLEX) 347 Scalar cinitslope,clambda; 348 #endif 349 int ierr,count; 350 SNES_LS *neP = (SNES_LS *) snes->data; 351 Scalar mone = -1.0,scale; 352 353 PLogEventBegin(SNES_LineSearch,snes,x,f,g); 354 *flag = 0; 355 alpha = neP->alpha; 356 maxstep = neP->maxstep; 357 steptol = neP->steptol; 358 359 VecNorm(y, NORM_2,ynorm ); 360 if (*ynorm > maxstep) { /* Step too big, so scale back */ 361 scale = maxstep/(*ynorm); 362 ierr = VecScale(&scale,y); CHKERRQ(ierr); 363 *ynorm = maxstep; 364 } 365 minlambda = steptol/(*ynorm); 366 ierr = MatMult(snes->jacobian,y,w); CHKERRQ(ierr); 367 #if defined(PETSC_COMPLEX) 368 ierr = VecDot(f,w,&cinitslope); CHKERRQ(ierr); 369 initslope = real(cinitslope); 370 #else 371 ierr = VecDot(f,w,&initslope); CHKERRQ(ierr); 372 #endif 373 if (initslope > 0.0) initslope = -initslope; 374 if (initslope == 0.0) initslope = -1.0; 375 376 ierr = VecCopy(y,w); CHKERRQ(ierr); 377 ierr = VecAYPX(&mone,x,w); CHKERRQ(ierr); 378 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 379 ierr = VecNorm(g,NORM_2,gnorm); CHKERRQ(ierr); 380 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 381 ierr = VecCopy(w,y); CHKERRQ(ierr); 382 PLogInfo((PetscObject)snes,"SNESQuadraticLineSearch: Using full step\n"); 383 goto theend; 384 } 385 386 /* Fit points with quadratic */ 387 lambda = 1.0; count = 0; 388 count = 1; 389 while (1) { 390 if (lambda <= minlambda) { /* bad luck; use full step */ 391 PLogInfo((PetscObject)snes, 392 "SNESQuadraticLineSearch:Unable to find good step length! %d \n",count); 393 PLogInfo((PetscObject)snes, 394 "SNESQuadraticLineSearch:f %g fnew %g ynorm %g lambda %g\n",fnorm,*gnorm,*ynorm,lambda); 395 ierr = VecCopy(w,y); CHKERRQ(ierr); 396 *flag = -1; break; 397 } 398 lambdatemp = -initslope/((*gnorm)*(*gnorm) - fnorm*fnorm - 2.0*initslope); 399 lambdaprev = lambda; 400 gnormprev = *gnorm; 401 if (lambdatemp <= .1*lambda) { 402 lambda = .1*lambda; 403 } else lambda = lambdatemp; 404 ierr = VecCopy(x,w); CHKERRQ(ierr); 405 lambda = -lambda; 406 #if defined(PETSC_COMPLEX) 407 clambda = lambda; ierr = VecAXPY(&clambda,y,w); CHKERRQ(ierr); 408 #else 409 ierr = VecAXPY(&lambda,y,w); CHKERRQ(ierr); 410 #endif 411 ierr = SNESComputeFunction(snes,w,g); CHKERRQ(ierr); 412 ierr = VecNorm(g,NORM_2,gnorm); CHKERRQ(ierr); 413 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 414 ierr = VecCopy(w,y); CHKERRQ(ierr); 415 PLogInfo((PetscObject)snes, 416 "SNESQuadraticLineSearch:Quadratically determined step, lambda %g\n",lambda); 417 break; 418 } 419 count++; 420 } 421 theend: 422 PLogEventEnd(SNES_LineSearch,snes,x,f,g); 423 return 0; 424 } 425 /* ------------------------------------------------------------ */ 426 /*@C 427 SNESSetLineSearchRoutine - Sets the line search routine to be used 428 by the method SNES_LS. 429 430 Input Parameters: 431 . snes - nonlinear context obtained from SNESCreate() 432 . func - pointer to int function 433 434 Available Routines: 435 . SNESCubicLineSearch() - default line search 436 . SNESQuadraticLineSearch() - quadratic line search 437 . SNESNoLineSearch() - the full Newton step (actually not a line search) 438 439 Options Database Keys: 440 $ -snes_line_search [basic,quadratic,cubic] 441 442 Calling sequence of func: 443 func (SNES snes, Vec x, Vec f, Vec g, Vec y, 444 Vec w, double fnorm, double *ynorm, 445 double *gnorm, *flag) 446 447 Input parameters for func: 448 . snes - nonlinear context 449 . x - current iterate 450 . f - residual evaluated at x 451 . y - search direction (contains new iterate on output) 452 . w - work vector 453 . fnorm - 2-norm of f 454 455 Output parameters for func: 456 . g - residual evaluated at new iterate y 457 . y - new iterate (contains search direction on input) 458 . gnorm - 2-norm of g 459 . ynorm - 2-norm of search length 460 . flag - set to 0 if the line search succeeds; a nonzero integer 461 on failure. 462 463 .keywords: SNES, nonlinear, set, line search, routine 464 465 .seealso: SNESNoLineSearch(), SNESQuadraticLineSearch(), SNESCubicLineSearch() 466 @*/ 467 int SNESSetLineSearchRoutine(SNES snes,int (*func)(SNES,Vec,Vec,Vec,Vec,Vec, 468 double,double*,double*,int*)) 469 { 470 if ((snes)->type == SNES_EQ_NLS) ((SNES_LS *)(snes->data))->LineSearch = func; 471 return 0; 472 } 473 /* ------------------------------------------------------------------ */ 474 static int SNESPrintHelp_LS(SNES snes,char *p) 475 { 476 SNES_LS *ls = (SNES_LS *)snes->data; 477 478 MPIU_printf(snes->comm," method ls (system of nonlinear equations):\n"); 479 MPIU_printf(snes->comm," %ssnes_line_search [basic,quadratic,cubic]\n",p); 480 MPIU_printf(snes->comm," %ssnes_line_search_alpha alpha (default %g)\n",p,ls->alpha); 481 MPIU_printf(snes->comm," %ssnes_line_search_maxstep max (default %g)\n",p,ls->maxstep); 482 MPIU_printf(snes->comm," %ssnes_line_search_steptol tol (default %g)\n",p,ls->steptol); 483 return 0; 484 } 485 /* ------------------------------------------------------------------ */ 486 static int SNESView_LS(PetscObject obj,Viewer viewer) 487 { 488 SNES snes = (SNES)obj; 489 SNES_LS *ls = (SNES_LS *)snes->data; 490 FILE *fd; 491 char *cstring; 492 int ierr; 493 494 ierr = ViewerFileGetPointer_Private(viewer,&fd); CHKERRQ(ierr); 495 if (ls->LineSearch == SNESNoLineSearch) cstring = "SNESNoLineSearch"; 496 else if (ls->LineSearch == SNESQuadraticLineSearch) cstring = "SNESQuadraticLineSearch"; 497 else if (ls->LineSearch == SNESCubicLineSearch) cstring = "SNESCubicLineSearch"; 498 else cstring = "unknown"; 499 MPIU_fprintf(snes->comm,fd," line search variant: %s\n",cstring); 500 MPIU_fprintf(snes->comm,fd," alpha=%g, maxstep=%g, steptol=%g\n", 501 ls->alpha,ls->maxstep,ls->steptol); 502 return 0; 503 } 504 /* ------------------------------------------------------------------ */ 505 static int SNESSetFromOptions_LS(SNES snes) 506 { 507 SNES_LS *ls = (SNES_LS *)snes->data; 508 char ver[16]; 509 double tmp; 510 int ierr,flg; 511 512 ierr = OptionsGetDouble(snes->prefix,"-snes_line_search_alpa",&tmp, &flg);CHKERRQ(ierr); 513 if (flg) { 514 ls->alpha = tmp; 515 } 516 ierr = OptionsGetDouble(snes->prefix,"-snes_line_search_maxstep",&tmp, &flg);CHKERRQ(ierr); 517 if (flg) { 518 ls->maxstep = tmp; 519 } 520 ierr = OptionsGetDouble(snes->prefix,"-snes_line_search_steptol",&tmp, &flg);CHKERRQ(ierr); 521 if (flg) { 522 ls->steptol = tmp; 523 } 524 ierr = OptionsGetString(snes->prefix,"-snes_line_search",ver,16, &flg); CHKERRQ(ierr); 525 if (flg) { 526 if (!PetscStrcmp(ver,"basic")) { 527 SNESSetLineSearchRoutine(snes,SNESNoLineSearch); 528 } 529 else if (!PetscStrcmp(ver,"quadratic")) { 530 SNESSetLineSearchRoutine(snes,SNESQuadraticLineSearch); 531 } 532 else if (!PetscStrcmp(ver,"cubic")) { 533 SNESSetLineSearchRoutine(snes,SNESCubicLineSearch); 534 } 535 else {SETERRQ(1,"SNESSetFromOptions_LS:Unknown line search");} 536 } 537 return 0; 538 } 539 /* ------------------------------------------------------------ */ 540 int SNESCreate_LS(SNES snes ) 541 { 542 SNES_LS *neP; 543 544 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 545 SETERRQ(1,"SNESCreate_LS:For SNES_NONLINEAR_EQUATIONS only"); 546 snes->type = SNES_EQ_NLS; 547 snes->setup = SNESSetUp_LS; 548 snes->solve = SNESSolve_LS; 549 snes->destroy = SNESDestroy_LS; 550 snes->converged = SNESDefaultConverged; 551 snes->printhelp = SNESPrintHelp_LS; 552 snes->setfromoptions = SNESSetFromOptions_LS; 553 snes->view = SNESView_LS; 554 555 neP = PetscNew(SNES_LS); CHKPTRQ(neP); 556 PLogObjectMemory(snes,sizeof(SNES_LS)); 557 snes->data = (void *) neP; 558 neP->alpha = 1.e-4; 559 neP->maxstep = 1.e8; 560 neP->steptol = 1.e-12; 561 neP->LineSearch = SNESCubicLineSearch; 562 return 0; 563 } 564 565 566 567 568