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