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