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_Polynomial.H"
00043
00044 #include "NOX_LineSearch_Utils_Printing.H"
00045 #include "NOX_LineSearch_Utils_Slope.H"
00046 #include "NOX_Abstract_Vector.H"
00047 #include "NOX_Abstract_Group.H"
00048 #include "NOX_Solver_Generic.H"
00049 #include "Teuchos_ParameterList.hpp"
00050 #include "NOX_MeritFunction_Generic.H"
00051 #include "NOX_StatusTest_FiniteValue.H"
00052 #include "NOX_GlobalData.H"
00053
00054 NOX::LineSearch::Polynomial::
00055 Polynomial(const Teuchos::RCP<NOX::GlobalData>& gd,
00056 Teuchos::ParameterList& params) :
00057 globalDataPtr(gd),
00058 paramsPtr(NULL),
00059 print(gd->getUtils()),
00060 slopeUtil(gd)
00061 {
00062 reset(gd, params);
00063 }
00064
00065 NOX::LineSearch::Polynomial::~Polynomial()
00066 {
00067
00068 }
00069
00070 bool NOX::LineSearch::Polynomial::
00071 reset(const Teuchos::RCP<NOX::GlobalData>& gd,
00072 Teuchos::ParameterList& params)
00073 {
00074 globalDataPtr = gd;
00075 meritFuncPtr = gd->getMeritFunction();
00076 print.reset(gd->getUtils());
00077 paramsPtr = ¶ms;
00078 slopeUtil.reset(gd);
00079
00080 Teuchos::ParameterList& p = params.sublist("Polynomial");
00081
00082 string choice = p.get("Sufficient Decrease Condition", "Armijo-Goldstein");
00083
00084 if (choice == "Armijo-Goldstein")
00085 suffDecrCond = ArmijoGoldstein;
00086 else if (choice == "Ared/Pred")
00087 suffDecrCond = AredPred;
00088 else if (choice == "None")
00089 suffDecrCond = None;
00090 else
00091 {
00092 print.err() << "NOX::LineSearch::Polynomial::reset - Invalid \"Sufficient Decrease Condition\"" << endl;
00093 throw "NOX Error";
00094 }
00095
00096 choice = p.get("Interpolation Type", "Cubic");
00097
00098 if (choice == "Cubic")
00099 interpolationType = Cubic;
00100 else if (choice == "Quadratic")
00101 interpolationType = Quadratic;
00102 else if (choice == "Quadratic3")
00103 interpolationType = Quadratic3;
00104 else
00105 {
00106 print.err() << "NOX::LineSearch::Polynomial::reset - Invalid \"Interpolation Type\"" << endl;
00107 throw "NOX Error";
00108 }
00109
00110 choice = p.get("Recovery Step Type", "Constant");
00111
00112 if (choice == "Constant")
00113 recoveryStepType = Constant;
00114 else if (choice == "Last Computed Step") {
00115 recoveryStepType = LastComputedStep;
00116 }
00117 else {
00118 print.err() << "NOX::LineSearch::Polynomial::reset - Invalid \"Recovery Step Type\"" << endl;
00119 throw "NOX Error";
00120 }
00121
00122 minStep = p.get("Minimum Step", 1.0e-12);
00123 defaultStep = p.get("Default Step", 1.0);
00124 recoveryStep = p.get("Recovery Step", defaultStep);
00125 maxIters = p.get("Max Iters", 100);
00126 alpha = p.get("Alpha Factor", 1.0e-4);
00127 minBoundFactor = p.get("Min Bounds Factor", 0.1);
00128 maxBoundFactor = p.get("Max Bounds Factor", 0.5);
00129 doForceInterpolation = p.get("Force Interpolation", false);
00130 useCounter = p.get("Use Counters", true);
00131 maxIncreaseIter = p.get("Maximum Iteration for Increase", 0);
00132 maxRelativeIncrease = p.get("Allowed Relative Increase", 1.e2);
00133
00134
00135 doAllowIncrease = (maxIncreaseIter > 0);
00136
00137
00138 if (useCounter)
00139 counter.reset();
00140
00141 return true;
00142 }
00143
00144 bool NOX::LineSearch::Polynomial::compute(Abstract::Group& newGrp,
00145 double& step,
00146 const Abstract::Vector& dir,
00147 const Solver::Generic& s)
00148 {
00149 printOpeningRemarks();
00150
00151 int nNonlinearIters = s.getNumIterations();
00152
00153 if (useCounter)
00154 counter.incrementNumLineSearches();
00155
00156
00157 string direction = const_cast<Teuchos::ParameterList&>(s.getList()).
00158 sublist("Direction").get("Method", "Newton");
00159 double eta = (suffDecrCond == AredPred) ?
00160 const_cast<Teuchos::ParameterList&>(s.getList()).
00161 sublist("Direction").sublist(direction).sublist("Linear Solver").
00162 get("Tolerance", -1.0) : 0.0;
00163
00164
00165 const Abstract::Group& oldGrp = s.getPreviousSolutionGroup();
00166 double oldPhi = meritFuncPtr->computef(oldGrp);
00167 double oldValue = computeValue(oldGrp, oldPhi);
00168 double oldSlope = meritFuncPtr->computeSlope(dir, oldGrp);
00169
00170
00171 step = defaultStep;
00172 updateGrp(newGrp, oldGrp, dir, step);
00173 double newPhi = meritFuncPtr->computef(newGrp);
00174 double newValue = computeValue(newGrp, newPhi);
00175
00176 bool isConverged = false;
00177 bool isFailed = false;
00178 int nIters = 1;
00179
00180 if (oldSlope >= 0.0)
00181 {
00182 printBadSlopeWarning(oldSlope);
00183 isFailed = true;
00184 }
00185 else
00186 isConverged = checkConvergence(newValue, oldValue, oldSlope, step,
00187 eta, nIters, nNonlinearIters);
00188
00189
00190 if ((useCounter) && (!isConverged))
00191 counter.incrementNumNonTrivialLineSearches();
00192
00193 double prevPhi = 0.0;
00194 double prevPrevPhi = 0.0;
00195 double prevStep = 0.0;
00196 double prevPrevStep = 0.0;
00197
00198 while ((!isConverged) && (!isFailed))
00199 {
00200 print.printStep(nIters, step, oldValue, newValue,
00201 "", (suffDecrCond != AredPred));
00202
00203 if (nIters > maxIters)
00204 {
00205 isFailed = true;
00206 break;
00207 }
00208
00209 if (interpolationType == Quadratic3)
00210 {
00211
00212
00213 prevPrevPhi = prevPhi;
00214 prevPhi = newPhi;
00215 prevPrevStep = prevStep;
00216 prevStep = step;
00217
00218 if (nIters == 1)
00219 {
00220 step = 0.5 * step;
00221 }
00222 else
00223 {
00224 double c1 = prevStep * prevStep * (prevPrevPhi - oldPhi) -
00225 prevPrevStep * prevPrevStep * (prevPhi - oldPhi);
00226 double c2 = prevPrevStep * (prevPhi - oldPhi) -
00227 prevStep * (prevPrevPhi - oldPhi);
00228
00229 if (c1 < 0)
00230 step = -0.5 * c1 / c2;
00231 }
00232 }
00233
00234 else if ((nIters == 1) || (interpolationType == Quadratic))
00235 {
00236
00237
00238 prevPhi = newPhi;
00239 prevStep = step;
00240
00241 step = - (oldSlope * prevStep * prevStep) /
00242 (2.0 * (prevPhi - oldPhi - prevStep * oldSlope)) ;
00243
00244 }
00245
00246 else
00247 {
00248
00249
00250 prevPrevPhi = prevPhi;
00251 prevPhi = newPhi;
00252 prevPrevStep = prevStep;
00253 prevStep = step;
00254
00255 double term1 = prevPhi - oldPhi - prevStep * oldSlope ;
00256 double term2 = prevPrevPhi - oldPhi - prevPrevStep * oldSlope ;
00257
00258 double a = 1.0 / (prevStep - prevPrevStep) *
00259 (term1 / (prevStep * prevStep) - term2 /
00260 (prevPrevStep * prevPrevStep)) ;
00261
00262 double b = 1.0 / (prevStep - prevPrevStep) *
00263 (-1.0 * term1 * prevPrevStep / (prevStep * prevStep) +
00264 term2 * prevStep / (prevPrevStep * prevPrevStep)) ;
00265
00266 double disc = b * b - 3.0 * a * oldSlope;
00267
00268 if (disc < 0)
00269 {
00270 isFailed = true;
00271 break;
00272 }
00273
00274 if (b > 0.0)
00275 {
00276 step = -oldSlope / (b + sqrt(disc));
00277 }
00278 else
00279 {
00280 if (fabs(a) < 1.e-12)
00281 {
00282 step = -oldSlope / (2.0 * b);
00283 }
00284 else
00285 {
00286 step = (-b + sqrt(disc))/ (3.0 * a);
00287 }
00288 }
00289 }
00290
00291
00292 if (step < minBoundFactor * prevStep)
00293 step = minBoundFactor * prevStep;
00294 else if (step > maxBoundFactor * prevStep)
00295 step = maxBoundFactor * prevStep;
00296
00297
00298 if (step < minStep)
00299 {
00300 isFailed = true;
00301 break;
00302 }
00303
00304
00305 updateGrp(newGrp, oldGrp, dir, step);
00306 newPhi = meritFuncPtr->computef(newGrp);
00307 newValue = computeValue(newGrp, newPhi);
00308
00309 nIters ++;
00310
00311 if (useCounter)
00312 counter.incrementNumIterations();
00313
00314 isConverged = checkConvergence(newValue, oldValue, oldSlope, step,
00315 eta, nIters, nNonlinearIters);
00316
00317 }
00318
00319
00320 if (isFailed)
00321 {
00322 if (useCounter)
00323 counter.incrementNumFailedLineSearches();
00324
00325 if (recoveryStepType == Constant)
00326 step = recoveryStep;
00327
00328 if (step == 0.0)
00329 {
00330 newGrp = oldGrp;
00331 newPhi = oldPhi;
00332 newValue = oldValue;
00333 }
00334 else
00335 {
00336 updateGrp(newGrp, oldGrp, dir, step);
00337 newPhi = meritFuncPtr->computef(newGrp);
00338 newValue = computeValue(newGrp, newPhi);
00339 }
00340 }
00341
00342 string message = (isFailed) ? "(USING RECOVERY STEP!)" : "(STEP ACCEPTED!)";
00343 print.printStep(nIters, step, oldValue, newValue, message, (suffDecrCond != AredPred));
00344
00345 paramsPtr->set("Adjusted Tolerance", 1.0 - step * (1.0 - eta));
00346
00347 if (useCounter)
00348 counter.setValues(*paramsPtr);
00349
00350 return (!isFailed);
00351 }
00352
00353 bool NOX::LineSearch::Polynomial::
00354 checkConvergence(double newValue, double oldValue,
00355 double oldSlope,
00356 double step, double eta,
00357 int nIters,
00358 int nNonlinearIters) const
00359 {
00360 NOX::StatusTest::FiniteValue checkNAN;
00361 if (checkNAN.finiteNumberTest(newValue) !=0)
00362 return false;
00363
00364 if ((nIters == 1) && (doForceInterpolation))
00365 return false;
00366
00367 if ((doAllowIncrease) && (nNonlinearIters <= maxIncreaseIter))
00368 {
00369 double relativeIncrease = newValue / oldValue;
00370 if (relativeIncrease < maxRelativeIncrease)
00371 return true;
00372 }
00373
00374 bool returnVal = false;
00375 switch (suffDecrCond)
00376 {
00377
00378 case ArmijoGoldstein:
00379 returnVal = (newValue <= oldValue + alpha * step * oldSlope);
00380 break;
00381 case AredPred:
00382 {
00383 double newEta = 1.0 - step * (1.0 - eta);
00384 returnVal = (newValue <= oldValue * (1.0 - alpha * (1.0 - newEta)));
00385 break;
00386 }
00387 case None:
00388 returnVal = true;
00389 break;
00390 default:
00391
00392 print.err() << "NOX::LineSearch::Polynomial::isSufficientDecrease - Unknown convergence criteria" << endl;
00393 throw "NOX Error";
00394
00395 }
00396 return returnVal;
00397 }
00398
00399 bool NOX::LineSearch::Polynomial::
00400 updateGrp(NOX::Abstract::Group& newGrp,
00401 const NOX::Abstract::Group& oldGrp,
00402 const NOX::Abstract::Vector& dir,
00403 double step) const
00404 {
00405 newGrp.computeX(oldGrp, dir, step);
00406
00407 NOX::Abstract::Group::ReturnType status = newGrp.computeF();
00408 if (status != NOX::Abstract::Group::Ok)
00409 return false;
00410
00411 return true;
00412 }
00413
00414 double NOX::LineSearch::Polynomial::
00415 computeValue(const NOX::Abstract::Group& grp, double phi)
00416 {
00417 double value = phi;
00418
00419 if (suffDecrCond == AredPred)
00420 {
00421 value = grp.getNormF();
00422 }
00423
00424 return value;
00425 }
00426
00427 void NOX::LineSearch::Polynomial::printOpeningRemarks() const
00428 {
00429 if (print.isPrintType(NOX::Utils::InnerIteration))
00430 {
00431 print.out() << "\n" << NOX::Utils::fill(72) << "\n"
00432 << "-- Polynomial Line Search -- \n";
00433 }
00434
00435 if (print.isPrintType(NOX::Utils::Details))
00436 {
00437 if (!Teuchos::is_null(meritFuncPtr))
00438 print.out() << " Merit Function = " << meritFuncPtr->name()
00439 << endl;
00440 }
00441 }
00442
00443 void NOX::LineSearch::Polynomial::printBadSlopeWarning(double slope) const
00444 {
00445 print.out(NOX::Utils::Warning)
00446 << "WARNING: Computed slope is positive (slope = "
00447 << slope << ").\n" << "Using recovery step!" << endl;
00448 }