00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "NOX_LineSearch_MoreThuente.H"
00043
00044 #include "NOX_Abstract_Vector.H"
00045 #include "NOX_Abstract_Group.H"
00046 #include "NOX_Solver_Generic.H"
00047 #include "Teuchos_ParameterList.hpp"
00048 #include "NOX_Utils.H"
00049 #include "NOX_MeritFunction_Generic.H"
00050 #include "NOX_GlobalData.H"
00051
00052 NOX::LineSearch::MoreThuente::
00053 MoreThuente(const Teuchos::RCP<NOX::GlobalData>& gd,
00054 Teuchos::ParameterList& params) :
00055 globalDataPtr(gd),
00056 print(gd->getUtils()),
00057 slope(gd),
00058 paramsPtr(0)
00059 {
00060 reset(gd, params);
00061 }
00062
00063 NOX::LineSearch::MoreThuente::~MoreThuente()
00064 {
00065 }
00066
00067 bool NOX::LineSearch::MoreThuente::
00068 reset(const Teuchos::RCP<NOX::GlobalData>& gd,
00069 Teuchos::ParameterList& params)
00070 {
00071 globalDataPtr = gd;
00072 meritFuncPtr = gd->getMeritFunction();
00073 print.reset(gd->getUtils());
00074 slope.reset(gd);
00075
00076 paramsPtr = ¶ms;
00077 Teuchos::ParameterList& p = params.sublist("More'-Thuente");
00078 ftol = p.get("Sufficient Decrease", 1.0e-4);
00079 gtol = p.get("Curvature Condition", 0.9999);
00080 xtol = p.get("Interval Width", 1.0e-15);
00081 stpmin = p.get("Minimum Step", 1.0e-12);
00082 stpmax = p.get("Maximum Step", 1.0e+6);
00083 maxfev = p.get("Max Iters", 20);
00084 defaultstep = p.get("Default Step", 1.0);
00085 recoverystep = p.get("Recovery Step", defaultstep);
00086
00087
00088 if ((ftol < 0.0) ||
00089 (gtol < 0.0) ||
00090 (xtol < 0.0) ||
00091 (stpmin < 0.0) ||
00092 (stpmax < stpmin) ||
00093 (maxfev <= 0) ||
00094 (defaultstep <= 0))
00095 {
00096 print.out() << "NOX::LineSearch::MoreThuente::reset - Error in Input Parameter!" << endl;
00097 throw "NOX Error";
00098 }
00099
00100 counter.reset();
00101
00102 string choice = p.get("Sufficient Decrease Condition", "Armijo-Goldstein");
00103 if (choice == "Ared/Pred")
00104 suffDecrCond = AredPred;
00105 else if (choice == "Armijo-Goldstein")
00106 suffDecrCond = ArmijoGoldstein;
00107 else {
00108 print.out() << "ERROR: NOX::LineSearch::MoreThuente::reset() - the choice of "
00109 << "\"Sufficient Decrease Condition\" is invalid." << endl;
00110 throw "NOX Error";
00111 }
00112
00113 choice = p.get("Recovery Step Type", "Constant");
00114 if (choice == "Constant")
00115 recoveryStepType = Constant;
00116 else if (choice == "Last Computed Step") {
00117 recoveryStepType = LastComputedStep;
00118 }
00119 else {
00120 print.out() << "NOX::LineSearch::MoreThuente::reset - Invalid "
00121 << "\"Recovery Step Type\"" << endl;
00122 throw "NOX Error";
00123 }
00124
00125 useOptimizedSlopeCalc = p.get("Optimize Slope Calculation", false);
00126
00127 return true;
00128 }
00129
00130
00131 bool NOX::LineSearch::MoreThuente::compute(Abstract::Group& grp, double& step,
00132 const Abstract::Vector& dir,
00133 const Solver::Generic& s)
00134 {
00135 counter.incrementNumLineSearches();
00136 const Abstract::Group& oldGrp = s.getPreviousSolutionGroup();
00137 int info = cvsrch(grp, step, oldGrp, dir, s);
00138
00139 if (step != 1.0)
00140 counter.incrementNumNonTrivialLineSearches();
00141
00142 counter.setValues(*paramsPtr);
00143
00144 return (info == 1);
00145 }
00146
00147 int NOX::LineSearch::MoreThuente::
00148 cvsrch(Abstract::Group& newgrp, double& stp, const Abstract::Group& oldgrp,
00149 const Abstract::Vector& dir, const Solver::Generic& s)
00150 {
00151 if (print.isPrintType(NOX::Utils::InnerIteration))
00152 {
00153 print.out() << "\n" << NOX::Utils::fill(72) << "\n" << "-- More'-Thuente Line Search -- \n";
00154 }
00155
00156
00157 stp = defaultstep;
00158
00159 int info = 0;
00160 int infoc = 1;
00161
00162
00163
00164 double dginit = 0.0;
00165 if (useOptimizedSlopeCalc)
00166 dginit = slope.computeSlopeWithOutJac(dir, oldgrp);
00167 else
00168 dginit = slope.computeSlope(dir, oldgrp);
00169
00170 if (dginit >= 0.0)
00171 {
00172 if (print.isPrintType(NOX::Utils::Warning))
00173 {
00174 print.out() << "NOX::LineSearch::MoreThuente::cvsrch - Non-descent direction (dginit = " << dginit << ")" << endl;
00175 }
00176 stp = recoverystep;
00177 newgrp.computeX(oldgrp, dir, stp);
00178 return 7;
00179 }
00180
00181
00182
00183 bool brackt = false;
00184 bool stage1 = true;
00185 int nfev = 0;
00186 double dgtest = ftol * dginit;
00187 double width = stpmax - stpmin;
00188 double width1 = 2 * width;
00189
00190
00191 double finit = meritFuncPtr->computef(oldgrp);
00192
00193
00194
00195
00196
00197
00198
00199
00200 double stx = 0.0;
00201 double fx = finit;
00202 double dgx = dginit;
00203 double sty = 0.0;
00204 double fy = finit;
00205 double dgy = dginit;
00206
00207
00208 const Teuchos::ParameterList& p = s.getList();
00209 double eta_original = 0.0;
00210 double eta = 0.0;
00211 eta_original = const_cast<Teuchos::ParameterList&>(p).sublist("Direction").
00212 sublist("Newton").sublist("Linear Solver").
00213 get(std::string("Tolerance"), -1.0);
00214 eta = eta_original;
00215
00216
00217
00218
00219 double stmin, stmax;
00220 double fm, fxm, fym, dgm, dgxm, dgym;
00221
00222 while (1)
00223 {
00224
00225
00226
00227
00228 if (brackt)
00229 {
00230 stmin = min(stx, sty);
00231 stmax = max(stx, sty);
00232 }
00233 else
00234 {
00235 stmin = stx;
00236 stmax = stp + 4 * (stp - stx);
00237 }
00238
00239
00240 stp = max(stp, stpmin);
00241 stp = min(stp, stpmax);
00242
00243
00244
00245
00246 if ((brackt && ((stp <= stmin) || (stp >= stmax))) ||
00247 (nfev >= maxfev - 1) || (infoc == 0) ||
00248 (brackt && (stmax - stmin <= xtol * stmax)))
00249 {
00250 stp = stx;
00251 }
00252
00253
00254
00255
00256 newgrp.computeX(oldgrp, dir, stp);
00257
00258 NOX::Abstract::Group::ReturnType rtype;
00259 rtype = newgrp.computeF();
00260 if (rtype != NOX::Abstract::Group::Ok)
00261 {
00262 print.err() << "NOX::LineSearch::MoreThuente::cvrch - Unable to compute F" << endl;
00263 throw "NOX Error";
00264 }
00265
00266 double f = meritFuncPtr->computef(newgrp);
00267
00268 if (!useOptimizedSlopeCalc) {
00269
00270 rtype = newgrp.computeJacobian();
00271 if (rtype != NOX::Abstract::Group::Ok)
00272 {
00273 print.err() << "NOX::LineSearch::MoreThuente::cvrch - Unable to compute Jacobian" << endl;
00274 throw "NOX Error";
00275 }
00276
00277 rtype = newgrp.computeGradient();
00278 if (rtype != NOX::Abstract::Group::Ok)
00279 {
00280 print.err() << "NOX::LineSearch::MoreThuente::cvrch - Unable to compute Gradient" << endl;
00281 throw "NOX Error";
00282 }
00283 }
00284
00285 nfev ++;
00286 string message = "";
00287
00288 double dg = 0.0;
00289 if (useOptimizedSlopeCalc)
00290 dg = slope.computeSlopeWithOutJac(dir, newgrp);
00291 else
00292 dg = slope.computeSlope(dir, newgrp);
00293
00294
00295 double ftest1 = finit + stp * dgtest;
00296
00297
00298 double ftest2 = 0.0;
00299 ftest2 = oldgrp.getNormF() * (1.0 - ftol * (1.0 - eta));
00300
00301
00302
00303 if ((brackt && ((stp <= stmin) || (stp >= stmax))) || (infoc == 0))
00304 info = 6;
00305
00306 if ((stp == stpmax) && (f <= ftest1) && (dg <= dgtest))
00307 info = 5;
00308
00309 if ((stp == stpmin) && ((f > ftest1) || (dg >= dgtest)))
00310 info = 4;
00311
00312 if (nfev >= maxfev)
00313 info = 3;
00314
00315 if (brackt && (stmax-stmin <= xtol*stmax))
00316 info = 2;
00317
00318
00319
00320
00321
00322
00323 bool sufficientDecreaseTest = false;
00324 if (suffDecrCond == ArmijoGoldstein)
00325 sufficientDecreaseTest = (f <= ftest1);
00326 else {
00327 double ap_normF = 0.0;
00328 ap_normF = newgrp.getNormF();
00329
00330 sufficientDecreaseTest = (ap_normF <= ftest2);
00331 }
00332
00333 if ((sufficientDecreaseTest) && (fabs(dg) <= gtol*(-dginit)))
00334 info = 1;
00335
00336 if (info != 0)
00337 {
00338 if (info != 1)
00339 {
00340
00341 counter.incrementNumFailedLineSearches();
00342
00343 if (recoveryStepType == Constant)
00344 stp = recoverystep;
00345
00346 newgrp.computeX(oldgrp, dir, stp);
00347
00348 message = "(USING RECOVERY STEP!)";
00349
00350
00351
00352
00353
00354
00355 }
00356 else
00357 {
00358 message = "(STEP ACCEPTED!)";
00359 }
00360
00361 print.printStep(nfev, stp, finit, f, message);
00362
00363
00364 return info;
00365
00366 }
00367
00368 print.printStep(nfev, stp, finit, f, message);
00369
00370
00371 counter.incrementNumIterations();
00372
00373
00374
00375
00376 if (stage1 && (f <= ftest1) && (dg >= min(ftol, gtol) * dginit))
00377 stage1 = false;
00378
00379
00380
00381
00382
00383
00384
00385 if (stage1 && (f <= fx) && (f > ftest1))
00386 {
00387
00388
00389
00390 fm = f - stp * dgtest;
00391 fxm = fx - stx * dgtest;
00392 fym = fy - sty * dgtest;
00393 dgm = dg - dgtest;
00394 dgxm = dgx - dgtest;
00395 dgym = dgy - dgtest;
00396
00397
00398
00399
00400 infoc = cstep(stx,fxm,dgxm,sty,fym,dgym,stp,fm,dgm,
00401 brackt,stmin,stmax);
00402
00403
00404
00405 fx = fxm + stx*dgtest;
00406 fy = fym + sty*dgtest;
00407 dgx = dgxm + dgtest;
00408 dgy = dgym + dgtest;
00409
00410 }
00411
00412 else
00413 {
00414
00415
00416
00417
00418 infoc = cstep(stx,fx,dgx,sty,fy,dgy,stp,f,dg,
00419 brackt,stmin,stmax);
00420
00421 }
00422
00423
00424
00425
00426 if (brackt)
00427 {
00428 if (fabs(sty - stx) >= 0.66 * width1)
00429 stp = stx + 0.5 * (sty - stx);
00430 width1 = width;
00431 width = fabs(sty-stx);
00432 }
00433
00434 }
00435
00436 }
00437
00438
00439 int NOX::LineSearch::MoreThuente::cstep(double& stx, double& fx, double& dx,
00440 double& sty, double& fy, double& dy,
00441 double& stp, double& fp, double& dp,
00442 bool& brackt, double stmin, double stmax)
00443 {
00444 int info = 0;
00445
00446
00447
00448 if ((brackt && ((stp <= min(stx, sty)) || (stp >= max(stx, sty)))) ||
00449 (dx * (stp - stx) >= 0.0) || (stmax < stmin))
00450 return info;
00451
00452
00453
00454 double sgnd = dp * (dx / fabs(dx));
00455
00456
00457
00458
00459
00460
00461 bool bound;
00462 double theta;
00463 double s;
00464 double gamma;
00465 double p,q,r;
00466 double stpc, stpq, stpf;
00467
00468 if (fp > fx)
00469 {
00470 info = 1;
00471 bound = 1;
00472 theta = 3 * (fx - fp) / (stp - stx) + dx + dp;
00473 s = absmax(theta, dx, dp);
00474 gamma = s * sqrt(((theta / s) * (theta / s)) - (dx / s) * (dp / s));
00475 if (stp < stx)
00476 gamma = -gamma;
00477
00478 p = (gamma - dx) + theta;
00479 q = ((gamma - dx) + gamma) + dp;
00480 r = p / q;
00481 stpc = stx + r * (stp - stx);
00482 stpq = stx + ((dx / ((fx - fp) / (stp - stx) + dx)) / 2) * (stp - stx);
00483 if (fabs(stpc - stx) < fabs(stpq - stx))
00484 stpf = stpc;
00485 else
00486 stpf = stpc + (stpq - stpc) / 2;
00487
00488 brackt = true;
00489 }
00490
00491
00492
00493
00494
00495
00496 else if (sgnd < 0.0)
00497 {
00498 info = 2;
00499 bound = false;
00500 theta = 3 * (fx - fp) / (stp - stx) + dx + dp;
00501 s = absmax(theta,dx,dp);
00502 gamma = s * sqrt(((theta/s) * (theta/s)) - (dx / s) * (dp / s));
00503 if (stp > stx)
00504 gamma = -gamma;
00505 p = (gamma - dp) + theta;
00506 q = ((gamma - dp) + gamma) + dx;
00507 r = p / q;
00508 stpc = stp + r * (stx - stp);
00509 stpq = stp + (dp / (dp - dx)) * (stx - stp);
00510 if (fabs(stpc - stp) > fabs(stpq - stp))
00511 stpf = stpc;
00512 else
00513 stpf = stpq;
00514 brackt = true;
00515 }
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526 else if (fabs(dp) < fabs(dx))
00527 {
00528 info = 3;
00529 bound = true;
00530 theta = 3 * (fx - fp) / (stp - stx) + dx + dp;
00531 s = absmax(theta, dx, dp);
00532
00533
00534
00535
00536 gamma = s * sqrt(max(0,(theta / s) * (theta / s) - (dx / s) * (dp / s)));
00537 if (stp > stx)
00538 gamma = -gamma;
00539
00540 p = (gamma - dp) + theta;
00541 q = (gamma + (dx - dp)) + gamma;
00542 r = p / q;
00543 if ((r < 0.0) && (gamma != 0.0))
00544 stpc = stp + r * (stx - stp);
00545 else if (stp > stx)
00546 stpc = stmax;
00547 else
00548 stpc = stmin;
00549
00550 stpq = stp + (dp/ (dp - dx)) * (stx - stp);
00551 if (brackt)
00552 {
00553 if (fabs(stp - stpc) < fabs(stp - stpq))
00554 stpf = stpc;
00555 else
00556 stpf = stpq;
00557 }
00558 else
00559 {
00560 if (fabs(stp - stpc) > fabs(stp - stpq))
00561 stpf = stpc;
00562 else
00563 stpf = stpq;
00564 }
00565 }
00566
00567
00568
00569
00570
00571
00572 else {
00573 info = 4;
00574 bound = false;
00575 if (brackt)
00576 {
00577 theta = 3 * (fp - fy) / (sty - stp) + dy + dp;
00578 s = absmax(theta, dy, dp);
00579 gamma = s * sqrt(((theta/s)*(theta/s)) - (dy / s) * (dp / s));
00580 if (stp > sty)
00581 gamma = -gamma;
00582 p = (gamma - dp) + theta;
00583 q = ((gamma - dp) + gamma) + dy;
00584 r = p / q;
00585 stpc = stp + r * (sty - stp);
00586 stpf = stpc;
00587 }
00588 else if (stp > stx)
00589 stpf = stmax;
00590 else
00591 stpf = stmin;
00592 }
00593
00594
00595
00596
00597 if (fp > fx)
00598 {
00599 sty = stp;
00600 fy = fp;
00601 dy = dp;
00602 }
00603 else
00604 {
00605 if (sgnd < 0.0)
00606 {
00607 sty = stx;
00608 fy = fx;
00609 dy = dx;
00610 }
00611 stx = stp;
00612 fx = fp;
00613 dx = dp;
00614 }
00615
00616
00617
00618 stpf = min(stmax, stpf);
00619 stpf = max(stmin, stpf);
00620 stp = stpf;
00621 if (brackt && bound)
00622 {
00623 if (sty > stx)
00624 stp = min(stx + 0.66 * (sty - stx), stp);
00625 else
00626 stp = max(stx + 0.66 * (sty - stx), stp);
00627 }
00628
00629 return info;
00630
00631 }
00632
00633 double NOX::LineSearch::MoreThuente::min(double a, double b)
00634 {
00635 return (a < b ? a : b);
00636 }
00637
00638 double NOX::LineSearch::MoreThuente::max(double a, double b)
00639 {
00640 return (a > b ? a : b);
00641 }
00642
00643
00644 double NOX::LineSearch::MoreThuente::absmax(double a, double b, double c)
00645 {
00646 a = fabs(a);
00647 b = fabs(b);
00648 c = fabs(c);
00649
00650 if (a > b)
00651 return (a > c) ? a : c;
00652 else
00653 return (b > c) ? b : c;
00654 }
00655
00656