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_Solver_TrustRegionBased.H"
00043 #include "NOX_Abstract_Vector.H"
00044 #include "NOX_Abstract_Group.H"
00045 #include "NOX_Common.H"
00046 #include "Teuchos_ParameterList.hpp"
00047 #include "NOX_MeritFunction_Generic.H"
00048 #include "NOX_Utils.H"
00049 #include "NOX_GlobalData.H"
00050 #include "NOX_Solver_SolverUtils.H"
00051 #include "NOX_Direction_Generic.H"
00052 #include "NOX_Direction_Factory.H"
00053
00054 using namespace NOX;
00055 using namespace NOX::Solver;
00056
00057 TrustRegionBased::
00058 TrustRegionBased(const Teuchos::RCP<NOX::Abstract::Group>& grp,
00059 const Teuchos::RCP<NOX::StatusTest::Generic>& t,
00060 const Teuchos::RCP<Teuchos::ParameterList>& p) :
00061 globalDataPtr(Teuchos::rcp(new NOX::GlobalData(p))),
00062 utilsPtr(globalDataPtr->getUtils()),
00063 solnPtr(grp),
00064 oldSolnPtr(grp->clone(DeepCopy)),
00065 newtonVecPtr(grp->getX().clone(ShapeCopy)),
00066 cauchyVecPtr(grp->getX().clone(ShapeCopy)),
00067 aVecPtr(grp->getX().clone(ShapeCopy)),
00068 bVecPtr(grp->getX().clone(ShapeCopy)),
00069 testPtr(t),
00070 paramsPtr(p),
00071 meritFuncPtr(globalDataPtr->getMeritFunction()),
00072 useAredPredRatio(false),
00073 prePostOperator(utilsPtr, paramsPtr->sublist("Solver Options"))
00074 {
00075 init();
00076 }
00077
00078
00079 void TrustRegionBased::init()
00080 {
00081
00082 nIter = 0;
00083 dx = 0;
00084 status = StatusTest::Unconverged;
00085
00086
00087 if (utilsPtr->isPrintType(NOX::Utils::Parameters)) {
00088 utilsPtr->out() << "\n" << NOX::Utils::fill(72) << "\n";
00089 utilsPtr->out() << "\n-- Parameters Passed to Nonlinear Solver --\n\n";
00090 paramsPtr->print(utilsPtr->out(),5);
00091 }
00092
00093
00094 paramsPtr->sublist("Direction").get("Method", "Newton");
00095 paramsPtr->sublist("Cauchy Direction").
00096 get("Method", "Steepest Descent");
00097 paramsPtr->sublist("Cauchy Direction").sublist("Steepest Descent").
00098 get("Scaling Type", "Quadratic Model Min");
00099
00100 newtonPtr = NOX::Direction::
00101 buildDirection(globalDataPtr, paramsPtr->sublist("Direction"));
00102
00103 cauchyPtr = NOX::Direction::
00104 buildDirection(globalDataPtr, paramsPtr->sublist("Cauchy Direction"));
00105
00106 minRadius = paramsPtr->sublist("Trust Region").
00107 get("Minimum Trust Region Radius", 1.0e-6);
00108 if (minRadius <= 0)
00109 invalid("Minimum Trust Region Radius", minRadius);
00110
00111 maxRadius = paramsPtr->sublist("Trust Region").get("Maximum Trust Region Radius", 1.0e+10);
00112 if (maxRadius <= minRadius)
00113 invalid("Maximum Trust Region Radius", maxRadius);
00114
00115 minRatio = paramsPtr->sublist("Trust Region").get("Minimum Improvement Ratio", 1.0e-4);
00116 if (minRatio <= 0)
00117 invalid("Minimum Improvement Ratio", minRatio);
00118
00119 contractTriggerRatio = paramsPtr->sublist("Trust Region").get("Contraction Trigger Ratio", 0.1);
00120 if (contractTriggerRatio < minRatio)
00121 invalid("Contraction Trigger Ratio", contractTriggerRatio);
00122
00123 expandTriggerRatio = paramsPtr->sublist("Trust Region").get("Expansion Trigger Ratio", 0.75);
00124 if (expandTriggerRatio <= contractTriggerRatio)
00125 invalid("Expansion Trigger Ratio", expandTriggerRatio);
00126
00127 contractFactor = paramsPtr->sublist("Trust Region").get("Contraction Factor", 0.25);
00128 if ((contractFactor <= 0) || (contractFactor >= 1))
00129 invalid("Contraction Factor", contractFactor);
00130
00131 expandFactor = paramsPtr->sublist("Trust Region").get("Expansion Factor", 4.0);
00132 if (expandFactor <= 1)
00133 invalid("Expansion Factor", expandFactor);
00134
00135 recoveryStep = paramsPtr->sublist("Trust Region").get("Recovery Step", 1.0);
00136 if (recoveryStep < 0)
00137 invalid("Recovery Step", recoveryStep);
00138
00139 checkType = parseStatusTestCheckType(paramsPtr->sublist("Solver Options"));
00140
00141
00142 useAredPredRatio =
00143 paramsPtr->sublist("Trust Region").
00144 get("Use Ared/Pred Ratio Calculation", false);
00145
00146 }
00147
00148
00149 void NOX::Solver::TrustRegionBased::
00150 invalid(const string& name, double value) const
00151 {
00152 utilsPtr->err() << "NOX::Solver::TrustRegionBased::init - "
00153 << "Invalid \"" << name << "\" (" << value << ")"
00154 << endl;
00155 throw "NOX Error";
00156 }
00157
00158 void TrustRegionBased::
00159 reset(const NOX::Abstract::Vector& initialGuess,
00160 const Teuchos::RCP<NOX::StatusTest::Generic>& t)
00161 {
00162
00163 solnPtr->setX(initialGuess);
00164 testPtr = t;
00165
00166
00167 nIter = 0;
00168 dx = 0;
00169 status = StatusTest::Unconverged;
00170
00171
00172 if (utilsPtr->isPrintType(NOX::Utils::Parameters)) {
00173 utilsPtr->out() << "\n" << NOX::Utils::fill(72) << "\n";
00174 utilsPtr->out() << "\n-- Parameters Passed to Nonlinear Solver --\n\n";
00175 paramsPtr->print(utilsPtr->out(),5);
00176 }
00177 }
00178
00179 void TrustRegionBased::
00180 reset(const NOX::Abstract::Vector& initialGuess)
00181 {
00182
00183 solnPtr->setX(initialGuess);
00184
00185
00186 nIter = 0;
00187 dx = 0;
00188 status = StatusTest::Unconverged;
00189
00190
00191 if (utilsPtr->isPrintType(NOX::Utils::Parameters)) {
00192 utilsPtr->out() << "\n" << NOX::Utils::fill(72) << "\n";
00193 utilsPtr->out() << "\n-- Parameters Passed to Nonlinear Solver --\n\n";
00194 paramsPtr->print(utilsPtr->out(),5);
00195 }
00196 }
00197
00198 TrustRegionBased::~TrustRegionBased()
00199 {
00200
00201 }
00202
00203
00204 NOX::StatusTest::StatusType TrustRegionBased::getStatus()
00205 {
00206 return status;
00207 }
00208
00209 NOX::StatusTest::StatusType TrustRegionBased::step()
00210 {
00211 prePostOperator.runPreIterate(*this);
00212
00213 if (nIter == 0) {
00214
00215 solnPtr->computeF();
00216 newF = meritFuncPtr->computef(*solnPtr);
00217
00218
00219 status = testPtr->checkStatus(*this, checkType);
00220
00221 printUpdate();
00222 }
00223
00224
00225 if (status != StatusTest::Unconverged)
00226 return status;
00227
00228
00229 Abstract::Group& soln = *solnPtr;
00230 StatusTest::Generic& test = *testPtr;
00231
00232
00233 bool ok;
00234 ok = newtonPtr->compute(*newtonVecPtr, soln, *this);
00235 if (!ok)
00236 {
00237 utilsPtr->out() << "NOX::Solver::TrustRegionBased::iterate - unable to calculate Newton direction" << endl;
00238 status = StatusTest::Failed;
00239 prePostOperator.runPostIterate(*this);
00240 printUpdate();
00241 return status;
00242 }
00243
00244 ok = cauchyPtr->compute(*cauchyVecPtr, soln, *this);
00245 if (!ok)
00246 {
00247 utilsPtr->err() << "NOX::Solver::TrustRegionBased::iterate - unable to calculate Cauchy direction" << endl;
00248 status = StatusTest::Failed;
00249 prePostOperator.runPostIterate(*this);
00250 printUpdate();
00251 return status;
00252 }
00253
00254 if (nIter == 0)
00255 {
00256 radius = newtonVecPtr->norm();
00257
00258 if (radius < minRadius)
00259 radius = 2 * minRadius;
00260 }
00261
00262
00263 nIter ++;
00264
00265
00266 *oldSolnPtr = *solnPtr;
00267
00268
00269 oldF = meritFuncPtr->computef(*oldSolnPtr);
00270
00271
00272 double ratio = -1;
00273
00274 if (utilsPtr->isPrintType(NOX::Utils::InnerIteration))
00275 {
00276 utilsPtr->out() << NOX::Utils::fill(72) << endl;
00277 utilsPtr->out() << "-- Trust Region Inner Iteration --" << endl;
00278 }
00279
00280
00281 while ((ratio < minRatio) && (radius > minRadius))
00282 {
00283
00284 Teuchos::RCP<NOX::Abstract::Vector> dirPtr;
00285 double step;
00286
00287
00288 double newtonVecNorm = newtonVecPtr->norm();
00289 double cauchyVecNorm = cauchyVecPtr->norm();
00290
00291 if (newtonVecNorm <= radius)
00292 {
00293 stepType = TrustRegionBased::Newton;
00294 step = 1.0;
00295 dirPtr = newtonVecPtr;
00296 }
00297 else if (cauchyVecNorm >= radius)
00298 {
00299 stepType = TrustRegionBased::Cauchy;
00300 step = radius / cauchyVecNorm;
00301 dirPtr = cauchyVecPtr;
00302 }
00303 else
00304 {
00305
00306
00307 aVecPtr->update(1.0, *newtonVecPtr, -1.0, *cauchyVecPtr, 0.0);
00308
00309
00310 double cta = cauchyVecPtr->innerProduct(*aVecPtr);
00311
00312 double ctc = cauchyVecPtr->innerProduct(*cauchyVecPtr);
00313
00314 double ata = aVecPtr->innerProduct(*aVecPtr);
00315
00316
00317 double tmp = (cta * cta) - ((ctc - (radius * radius)) * ata);
00318 if (tmp < 0) {
00319 utilsPtr->err() << "NOX::Solver::TrustRegionBased::iterate - invalid computation" << endl;
00320 throw "NOX Error";
00321 }
00322
00323
00324 double gamma = (sqrt(tmp) - cta) / ata;
00325 if ((gamma < 0) || (gamma > 1)) {
00326 utilsPtr->err() << "NOX::Solver::TrustRegionBased::iterate - invalid trust region step" << endl;
00327 throw "NOX Error";
00328 }
00329
00330
00331 aVecPtr->update(1.0 - gamma, *cauchyVecPtr, gamma, *newtonVecPtr, 0.0);
00332
00333
00334 stepType = TrustRegionBased::Dogleg;
00335 dirPtr = aVecPtr;
00336 step = 1.0;
00337 }
00338
00339
00340 const Abstract::Vector& dir = *dirPtr;
00341
00342
00343 dx = step * dir.norm();
00344
00345
00346 soln.computeX(*oldSolnPtr, dir, step);
00347
00348
00349 NOX::Abstract::Group::ReturnType rtype = soln.computeF();
00350 if (rtype != NOX::Abstract::Group::Ok)
00351 {
00352 utilsPtr->err() << "NOX::Solver::TrustRegionBased::iterate - unable to compute F" << endl;
00353 throw "NOX Error";
00354 }
00355
00356
00357
00358
00359 if (useAredPredRatio) {
00360
00361
00362 rtype = oldSolnPtr->applyJacobian(*dirPtr, *bVecPtr);
00363 if (rtype != NOX::Abstract::Group::Ok)
00364 {
00365 utilsPtr->out() << "NOX::Solver::TrustRegionBased::iterate - "
00366 << "unable to compute F" << endl;
00367 throw "NOX Error";
00368 }
00369 bVecPtr->update(1.0, oldSolnPtr->getF(), 1.0);
00370
00371
00372 double oldNormF = oldSolnPtr->getNormF();
00373 double newNormF = soln.getNormF();
00374 double normFLinear = bVecPtr->norm();
00375
00376 ratio = (oldNormF - newNormF) / (oldNormF - normFLinear);
00377
00378
00379 if (utilsPtr->isPrintType(NOX::Utils::InnerIteration)) {
00380 double numerator = oldNormF - newNormF;
00381 double denominator = oldNormF - normFLinear;
00382 utilsPtr->out() << "Ratio computation: "
00383 << utilsPtr->sciformat(numerator) << "/"
00384 << utilsPtr->sciformat(denominator) << "="
00385 << ratio << endl;
00386 }
00387
00388
00389 newF = meritFuncPtr->computef(*solnPtr);
00390
00391 }
00392 else {
00393
00394 newF = meritFuncPtr->computef(*solnPtr);
00395
00396 if (newF >= oldF)
00397 {
00398 ratio = -1;
00399 }
00400 else
00401 {
00402
00403 rtype = oldSolnPtr->applyJacobian(*dirPtr, *bVecPtr);
00404 if (rtype != NOX::Abstract::Group::Ok)
00405 {
00406 utilsPtr->err() << "NOX::Solver::TrustRegionBased::iterate - unable to compute F" << endl;
00407 throw "NOX Error";
00408 }
00409 double numerator = oldF - newF;
00410 double denominator = 0.0;
00411
00412 denominator = fabs(oldF - meritFuncPtr->
00413 computeQuadraticModel(dir,*oldSolnPtr));
00414
00415 ratio = numerator / denominator;
00416 if (utilsPtr->isPrintType(NOX::Utils::Debug))
00417 utilsPtr->out() << "Ratio computation: "
00418 << utilsPtr->sciformat(numerator) << "/"
00419 << utilsPtr->sciformat(denominator) << "="
00420 << utilsPtr->sciformat(ratio) << endl;
00421
00422
00423 if ((denominator < 1.0e-12) && ((newF / oldF) >= 0.5))
00424 ratio = -1;
00425 }
00426 }
00427
00428 if (utilsPtr->isPrintType(Utils::InnerIteration)) {
00429 utilsPtr->out() << "radius = " << utilsPtr->sciformat(radius, 1);
00430 utilsPtr->out() << " ratio = " << setprecision(2) << setw(4) << ratio;
00431 utilsPtr->out() << " f = " << utilsPtr->sciformat(sqrt(2*newF));
00432 utilsPtr->out() << " old f = " << utilsPtr->sciformat(sqrt(2*oldF));
00433 utilsPtr->out() << " ";
00434
00435 switch(stepType) {
00436 case TrustRegionBased::Newton:
00437 utilsPtr->out() << "Newton";
00438 break;
00439 case TrustRegionBased::Cauchy:
00440 utilsPtr->out() << "Cauchy";
00441 break;
00442 case TrustRegionBased::Dogleg:
00443 utilsPtr->out() << "Dogleg";
00444 break;
00445 }
00446
00447 utilsPtr->out() << endl;
00448 }
00449
00450
00451 if (ratio < contractTriggerRatio)
00452 {
00453 if (stepType == TrustRegionBased::Newton) {
00454 radius = newtonVecPtr->norm();
00455 }
00456 radius = NOX_MAX(contractFactor * radius, minRadius);
00457 }
00458 else if ((ratio > expandTriggerRatio) && (dx == radius))
00459 {
00460 radius = NOX_MIN(expandFactor * radius, maxRadius);
00461 }
00462
00463 }
00464
00465
00466
00467 if ((radius <= minRadius) && (ratio < minRatio))
00468 {
00469 if (utilsPtr->isPrintType(Utils::InnerIteration))
00470 utilsPtr->out() << "Using recovery step and resetting trust region." << endl;
00471 soln.computeX(*oldSolnPtr, *newtonVecPtr, recoveryStep);
00472 soln.computeF();
00473 radius = newtonVecPtr->norm();
00474
00475
00476 }
00477
00478 status = test.checkStatus(*this, checkType);
00479
00480 if (utilsPtr->isPrintType(Utils::InnerIteration))
00481 utilsPtr->out() << NOX::Utils::fill(72) << endl;
00482
00483 prePostOperator.runPostIterate(*this);
00484
00485 printUpdate();
00486
00487 return status;
00488 }
00489
00490 NOX::StatusTest::StatusType TrustRegionBased::solve()
00491 {
00492 prePostOperator.runPreSolve(*this);
00493
00494
00495 while (status == StatusTest::Unconverged) {
00496 status = step();
00497 }
00498
00499 Teuchos::ParameterList& outputParams = paramsPtr->sublist("Output");
00500 outputParams.set("Nonlinear Iterations", nIter);
00501 outputParams.set("2-Norm of Residual", solnPtr->getNormF());
00502
00503 prePostOperator.runPostSolve(*this);
00504
00505 return status;
00506 }
00507
00508 const Abstract::Group& TrustRegionBased::getSolutionGroup() const
00509 {
00510 return *solnPtr;
00511 }
00512
00513 const Abstract::Group& TrustRegionBased::getPreviousSolutionGroup() const
00514 {
00515 return *oldSolnPtr;
00516 }
00517
00518 int TrustRegionBased::getNumIterations() const
00519 {
00520 return nIter;
00521 }
00522
00523 const Teuchos::ParameterList& TrustRegionBased::getList() const
00524 {
00525 return *paramsPtr;
00526 }
00527
00528
00529 void TrustRegionBased::printUpdate()
00530 {
00531
00532 if ((status == StatusTest::Unconverged) &&
00533 (utilsPtr->isPrintType(NOX::Utils::OuterIterationStatusTest))) {
00534 utilsPtr->out() << NOX::Utils::fill(72) << "\n";
00535 utilsPtr->out() << "-- Status Test Results --\n";
00536 testPtr->print(utilsPtr->out());
00537 utilsPtr->out() << NOX::Utils::fill(72) << "\n";
00538 }
00539
00540 double fmax = solnPtr->getF().norm(Abstract::Vector::MaxNorm);
00541 if (utilsPtr->isPrintType(NOX::Utils::OuterIteration)) {
00542 utilsPtr->out() << "\n" << NOX::Utils::fill(72) << "\n";
00543 utilsPtr->out() << "-- Newton Trust-Region Step " << nIter << " -- \n";
00544 utilsPtr->out() << "f = " << utilsPtr->sciformat(sqrt(2*newF));
00545 utilsPtr->out() << " fmax = " << utilsPtr->sciformat(fmax);
00546 utilsPtr->out() << " dx = " << utilsPtr->sciformat(dx);
00547 utilsPtr->out() << " radius = " << utilsPtr->sciformat(radius);
00548 if (status == StatusTest::Converged)
00549 utilsPtr->out() << " (Converged!)";
00550 if (status == StatusTest::Failed)
00551 utilsPtr->out() << " (Failed!)";
00552 utilsPtr->out() << "\n" << NOX::Utils::fill(72) << "\n" << endl;
00553 }
00554
00555 if ((status != StatusTest::Unconverged) &&
00556 (utilsPtr->isPrintType(NOX::Utils::OuterIteration))) {
00557 utilsPtr->out() << NOX::Utils::fill(72) << "\n";
00558 utilsPtr->out() << "-- Final Status Test Results --\n";
00559 testPtr->print(utilsPtr->out());
00560 utilsPtr->out() << NOX::Utils::fill(72) << "\n";
00561 }
00562 }
00563