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