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