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_Common.H"
00043 #include "NOX_Solver_InexactTrustRegionBased.H"
00044 #include "NOX_Abstract_Vector.H"
00045 #include "NOX_Abstract_Group.H"
00046 #include "Teuchos_ParameterList.hpp"
00047 #include "NOX_Utils.H"
00048 #include "NOX_GlobalData.H"
00049 #include "NOX_MeritFunction_Generic.H"
00050 #include "NOX_Solver_PrePostOperator.H"
00051 #include "NOX_Solver_SolverUtils.H"
00052 #include "NOX_Direction_Generic.H"
00053 #include "NOX_Direction_Factory.H"
00054
00055 using namespace NOX;
00056 using namespace NOX::Solver;
00057
00058
00059
00060
00061 NOX::Solver::InexactTrustRegionBased::
00062 InexactTrustRegionBased(const Teuchos::RCP<NOX::Abstract::Group>& grp,
00063 const Teuchos::RCP<NOX::StatusTest::Generic>& t,
00064 const Teuchos::RCP<Teuchos::ParameterList>& p) :
00065 globalDataPtr(Teuchos::rcp(new NOX::GlobalData(p))),
00066 utils(globalDataPtr->getUtils()),
00067 solnPtr(grp),
00068 oldSolnPtr(grp->clone(DeepCopy)),
00069 newtonVecPtr(grp->getX().clone(ShapeCopy)),
00070 cauchyVecPtr(grp->getX().clone(ShapeCopy)),
00071 rCauchyVecPtr(grp->getX().clone(ShapeCopy)),
00072 residualVecPtr(grp->getX().clone(ShapeCopy)),
00073 aVecPtr(grp->getX().clone(ShapeCopy)),
00074 bVecPtr(grp->getX().clone(ShapeCopy)),
00075 testPtr(t),
00076 paramsPtr(p),
00077 inNewtonUtils(globalDataPtr, paramsPtr->sublist("Direction")),
00078 radius(0.0),
00079 meritFuncPtr(globalDataPtr->getMeritFunction()),
00080 useCauchyInNewtonDirection(false),
00081 writeOutputParamsToList(true),
00082 useCounters(true),
00083 numCauchySteps(0),
00084 numNewtonSteps(0),
00085 numDoglegSteps(0),
00086 numTrustRegionInnerIterations(0),
00087 sumDoglegFracCauchyToNewton(0.0),
00088 sumDoglegFracNewtonLength(0.0),
00089 useAredPredRatio(false),
00090 useDoglegMinimization(false),
00091 prePostOperator(utils, paramsPtr->sublist("Solver Options"))
00092 {
00093 init();
00094 }
00095
00096
00097
00098
00099 NOX::Solver::InexactTrustRegionBased::~InexactTrustRegionBased()
00100 {
00101
00102 }
00103
00104
00105
00106
00107
00108 void NOX::Solver::InexactTrustRegionBased::init()
00109 {
00110
00111 nIter = 0;
00112 dx = 0.0;
00113 status = StatusTest::Unconverged;
00114 if (useCounters)
00115 resetCounters();
00116
00117 checkType = parseStatusTestCheckType(paramsPtr->sublist("Solver Options"));
00118
00119
00120 if (utils->isPrintType(NOX::Utils::Parameters)) {
00121 utils->out() << "\n" << NOX::Utils::fill(72) << "\n";
00122 utils->out() << "\n-- Parameters Passed to Nonlinear Solver --\n\n";
00123 paramsPtr->print(utils->out(),5);
00124 }
00125
00126
00127 string methodChoice =
00128 paramsPtr->sublist("Trust Region").
00129 get("Inner Iteration Method", "Inexact Trust Region");
00130 if (methodChoice == "Standard Trust Region")
00131 method = Standard;
00132 else if (methodChoice == "Inexact Trust Region")
00133 method = Inexact;
00134 else {
00135 utils->err() << "NOX::Solver::InexactTrustRegionBased::init - \"" << methodChoice
00136 << "\" is an invalid choice for \"Method\" key!"
00137 << endl;
00138 throw "NOX Error";
00139 }
00140
00141
00142
00143 paramsPtr->sublist("Direction").get("Method", "Newton");
00144 paramsPtr->sublist("Cauchy Direction")
00145 .get("Method", "Steepest Descent");
00146 paramsPtr->sublist("Cauchy Direction").sublist("Steepest Descent")
00147 .get("Scaling Type", "Quadratic Model Min");
00148
00149 newtonPtr = NOX::Direction::
00150 buildDirection(globalDataPtr, paramsPtr->sublist("Direction"));
00151 cauchyPtr = NOX::Direction::
00152 buildDirection(globalDataPtr, paramsPtr->sublist("Cauchy Direction"));
00153 inNewtonUtils.reset(globalDataPtr, paramsPtr->sublist("Direction"));
00154
00155 minRadius = paramsPtr->sublist("Trust Region")
00156 .get("Minimum Trust Region Radius", 1.0e-6);
00157 if (minRadius <= 0.0)
00158 invalid("Minimum Trust Region Radius", minRadius);
00159
00160 maxRadius = paramsPtr->sublist("Trust Region")
00161 .get("Maximum Trust Region Radius", 1.0e+10);
00162 if (maxRadius <= minRadius)
00163 invalid("Maximum Trust Region Radius", maxRadius);
00164
00165 minRatio = paramsPtr->sublist("Trust Region")
00166 .get("Minimum Improvement Ratio", 1.0e-4);
00167 if (minRatio <= 0.0)
00168 invalid("Minimum Improvement Ratio", minRatio);
00169
00170 contractTriggerRatio = paramsPtr->sublist("Trust Region")
00171 .get("Contraction Trigger Ratio", 0.1);
00172 if (contractTriggerRatio < minRatio)
00173 invalid("Contraction Trigger Ratio", contractTriggerRatio);
00174
00175 expandTriggerRatio = paramsPtr->sublist("Trust Region")
00176 .get("Expansion Trigger Ratio", 0.75);
00177 if (expandTriggerRatio <= contractTriggerRatio)
00178 invalid("Expansion Trigger Ratio", expandTriggerRatio);
00179
00180 contractFactor = paramsPtr->sublist("Trust Region")
00181 .get("Contraction Factor", 0.25);
00182 if ((contractFactor <= 0.0) || (contractFactor >= 1))
00183 invalid("Contraction Factor", contractFactor);
00184
00185 expandFactor = paramsPtr->sublist("Trust Region")
00186 .get("Expansion Factor", 4.0);
00187 if (expandFactor <= 1.0)
00188 invalid("Expansion Factor", expandFactor);
00189
00190 recoveryStep = paramsPtr->sublist("Trust Region")
00191 .get("Recovery Step", 1.0);
00192 if (recoveryStep < 0.0)
00193 invalid("Recovery Step", recoveryStep);
00194
00195 useCauchyInNewtonDirection = paramsPtr->sublist("Trust Region")
00196 .get("Use Cauchy in Newton Direction", false);
00197
00198
00199 useAredPredRatio = paramsPtr->sublist("Trust Region")
00200 .get("Use Ared/Pred Ratio Calculation", false);
00201
00202
00203 useDoglegMinimization = paramsPtr->sublist("Trust Region")
00204 .get("Use Dogleg Segment Minimization", false);
00205
00206
00207 useCounters = paramsPtr->sublist("Trust Region")
00208 .get("Use Counters", true);
00209
00210
00211 useCounters = paramsPtr->sublist("Trust Region")
00212 .get("Write Output Parameters", true);
00213
00214 }
00215
00216
00217
00218
00219 void NOX::Solver::InexactTrustRegionBased::invalid(const string& name,
00220 double value) const
00221 {
00222 utils->out() << "NOX::Solver::InexactTrustRegionBased::init - "
00223 << "Invalid \"" << name << "\" (" << value << ")"
00224 << endl;
00225 throw "NOX Error";
00226 }
00227
00228
00229
00230
00231 void NOX::Solver::InexactTrustRegionBased::throwError(const string& method,
00232 const string& message) const
00233 {
00234 utils->out() << "NOX::Solver::InexactTrustRegionBased::" << method << " - "
00235 << message << endl;
00236 throw "NOX Error";
00237 }
00238
00239
00240
00241
00242 void NOX::Solver::InexactTrustRegionBased::
00243 reset(const NOX::Abstract::Vector& initialGuess,
00244 const Teuchos::RCP<NOX::StatusTest::Generic>& t)
00245 {
00246 solnPtr->setX(initialGuess);
00247 testPtr = t;
00248
00249
00250 nIter = 0;
00251 dx = 0.0;
00252 status = StatusTest::Unconverged;
00253 if (useCounters)
00254 resetCounters();
00255
00256
00257 if (utils->isPrintType(NOX::Utils::Parameters)) {
00258 utils->out() << "\n" << NOX::Utils::fill(72) << "\n";
00259 utils->out() << "\n-- Parameters Passed to Nonlinear Solver --\n\n";
00260 paramsPtr->print(utils->out(),5);
00261 }
00262
00263
00264 solnPtr->computeF();
00265 newF = meritFuncPtr->computef(*solnPtr);
00266
00267
00268 status = testPtr->checkStatus(*this, checkType);
00269
00270 if (utils->isPrintType(NOX::Utils::Parameters)) {
00271 utils->out() << "\n-- Status Tests Passed to Nonlinear Solver --\n\n";
00272 testPtr->print(utils->out(), 5);
00273 utils->out() <<"\n" << NOX::Utils::fill(72) << "\n";
00274 }
00275 }
00276
00277
00278
00279
00280 void NOX::Solver::InexactTrustRegionBased::
00281 reset(const NOX::Abstract::Vector& initialGuess)
00282 {
00283 solnPtr->setX(initialGuess);
00284
00285
00286 nIter = 0;
00287 dx = 0.0;
00288 status = StatusTest::Unconverged;
00289 if (useCounters)
00290 resetCounters();
00291
00292
00293 if (utils->isPrintType(NOX::Utils::Parameters)) {
00294 utils->out() << "\n" << NOX::Utils::fill(72) << "\n";
00295 utils->out() << "\n-- Parameters Passed to Nonlinear Solver --\n\n";
00296 paramsPtr->print(utils->out(),5);
00297 }
00298
00299
00300 solnPtr->computeF();
00301 newF = meritFuncPtr->computef(*solnPtr);
00302
00303
00304 status = testPtr->checkStatus(*this, checkType);
00305
00306 if (utils->isPrintType(NOX::Utils::Parameters)) {
00307 utils->out() << "\n-- Status Tests Passed to Nonlinear Solver --\n\n";
00308 testPtr->print(utils->out(), 5);
00309 utils->out() <<"\n" << NOX::Utils::fill(72) << "\n";
00310 }
00311 }
00312
00313
00314
00315
00316 NOX::StatusTest::StatusType NOX::Solver::InexactTrustRegionBased::getStatus()
00317 {
00318 return status;
00319 }
00320
00321
00322
00323
00324 NOX::StatusTest::StatusType NOX::Solver::InexactTrustRegionBased::step()
00325 {
00326 prePostOperator.runPreIterate(*this);
00327
00328 if (nIter == 0) {
00329
00330 solnPtr->computeF();
00331 newF = meritFuncPtr->computef(*solnPtr);
00332
00333
00334 status = testPtr->checkStatus(*this, checkType);
00335
00336 printUpdate();
00337 }
00338
00339 NOX::StatusTest::StatusType status = NOX::StatusTest::Unconverged;
00340
00341 switch(method) {
00342 case Standard:
00343 status = iterateStandard();
00344 break;
00345
00346 case Inexact:
00347 status = iterateInexact();
00348 break;
00349
00350 default:
00351 status = iterateInexact();
00352 break;
00353 }
00354
00355 prePostOperator.runPostIterate(*this);
00356
00357 printUpdate();
00358
00359 return status;
00360 }
00361
00362
00363
00364
00365 NOX::StatusTest::StatusType
00366 NOX::Solver::InexactTrustRegionBased::iterateStandard()
00367 {
00368
00369 if (status != StatusTest::Unconverged)
00370 return status;
00371
00372
00373 Abstract::Group& soln = *solnPtr;
00374 StatusTest::Generic& test = *testPtr;
00375
00376
00377 bool ok;
00378 ok = newtonPtr->compute(*newtonVecPtr, soln, *this);
00379 if (!ok)
00380 {
00381 utils->out() << "NOX::Solver::InexactTrustRegionBased::iterate - "
00382 << "unable to calculate Newton direction" << endl;
00383 status = StatusTest::Failed;
00384 return status;
00385 }
00386
00387 ok = cauchyPtr->compute(*cauchyVecPtr, soln, *this);
00388 if (!ok)
00389 {
00390 utils->out() << "NOX::Solver::InexactTrustRegionBased::iterate - "
00391 << "unable to calculate Cauchy direction" << endl;
00392 status = StatusTest::Failed;
00393 return status;
00394 }
00395
00396 if (nIter == 0)
00397 {
00398 radius = computeNorm(*newtonVecPtr);
00399
00400 if (radius < minRadius)
00401 radius = 2 * minRadius;
00402 }
00403
00404
00405 nIter ++;
00406
00407
00408 *oldSolnPtr = *solnPtr;
00409
00410
00411
00412 oldF = meritFuncPtr->computef(*oldSolnPtr);
00413
00414
00415 double ratio = -1;
00416
00417 if (utils->isPrintType(NOX::Utils::InnerIteration))
00418 {
00419 utils->out() << NOX::Utils::fill(72) << endl;
00420 utils->out() << "-- Trust Region Inner Iteration --" << endl;
00421 }
00422
00423
00424 double gamma = 0.0;
00425
00426
00427 while ((ratio < minRatio) && (radius > minRadius))
00428 {
00429 if (useCounters)
00430 numTrustRegionInnerIterations += 1;
00431
00432 Teuchos::RCP<NOX::Abstract::Vector> dirPtr;
00433 double step;
00434
00435
00436 double newtonVecNorm = computeNorm(*newtonVecPtr);
00437 double cauchyVecNorm = computeNorm(*cauchyVecPtr);
00438
00439 if (newtonVecNorm <= radius)
00440 {
00441 stepType = InexactTrustRegionBased::Newton;
00442 step = 1.0;
00443 dirPtr = newtonVecPtr;
00444 }
00445 else if (cauchyVecNorm >= radius)
00446 {
00447 stepType = InexactTrustRegionBased::Cauchy;
00448 step = radius / cauchyVecNorm;
00449 dirPtr = cauchyVecPtr;
00450 }
00451 else
00452 {
00453
00454
00455 aVecPtr->update(1.0, *newtonVecPtr, -1.0, *cauchyVecPtr, 0.0);
00456
00457
00458 double cta = cauchyVecPtr->innerProduct(*aVecPtr);
00459
00460 double ctc = cauchyVecPtr->innerProduct(*cauchyVecPtr);
00461
00462 double ata = aVecPtr->innerProduct(*aVecPtr);
00463
00464
00465 double tmp = (cta * cta) - ((ctc - (radius * radius)) * ata);
00466 if (tmp < 0) {
00467 utils->err() << "NOX::Solver::InexactTrustRegionBased::iterate - invalid computation" << endl;
00468 throw "NOX Error";
00469 }
00470
00471
00472 gamma = (sqrt(tmp) - cta) / ata;
00473 if ((gamma < 0) || (gamma > 1)) {
00474 utils->err() << "NOX::Solver::InexactTrustRegionBased::iterate - invalid trust region step" << endl;
00475 throw "NOX Error";
00476 }
00477
00478
00479 aVecPtr->update(1.0 - gamma, *cauchyVecPtr, gamma, *newtonVecPtr, 0.0);
00480
00481
00482 stepType = InexactTrustRegionBased::Dogleg;
00483 dirPtr = aVecPtr;
00484 step = 1.0;
00485 }
00486
00487
00488 const NOX::Abstract::Vector& dir = *dirPtr;
00489
00490
00491 dx = step * (computeNorm(dir));
00492
00493
00494 soln.computeX(*oldSolnPtr, dir, step);
00495
00496
00497 NOX::Abstract::Group::ReturnType rtype = soln.computeF();
00498 if (rtype != NOX::Abstract::Group::Ok)
00499 {
00500 utils->err() << "NOX::Solver::InexactTrustRegionBased::iterate - "
00501 << "unable to compute F" << endl;
00502 throw "NOX Error";
00503 }
00504
00505
00506
00507
00508 if (useAredPredRatio) {
00509
00510
00511 rtype = oldSolnPtr->applyJacobian(*dirPtr, *bVecPtr);
00512 if (rtype != NOX::Abstract::Group::Ok)
00513 {
00514 utils->out() << "NOX::Solver::TrustRegionBased::iterate - "
00515 << "unable to compute F" << endl;
00516 throw "NOX Error";
00517 }
00518 bVecPtr->update(1.0, oldSolnPtr->getF(), 1.0);
00519
00520
00521 double oldNormF = computeNorm(oldSolnPtr->getF());
00522 double newNormF = computeNorm(solnPtr->getF());
00523 double normFLinear = computeNorm(*bVecPtr);;
00524
00525 ratio = (oldNormF - newNormF) / (oldNormF - normFLinear);
00526
00527
00528 if (utils->isPrintType(NOX::Utils::InnerIteration)) {
00529 double numerator = oldNormF - newNormF;
00530 double denominator = oldNormF - normFLinear;
00531 utils->out() << "Ratio computation: " << utils->sciformat(numerator)
00532 << "/" << utils->sciformat(denominator) << "="
00533 << ratio << endl;
00534 }
00535
00536
00537 newF = meritFuncPtr->computef(*solnPtr);
00538
00539 }
00540 else {
00541
00542 newF = meritFuncPtr->computef(*solnPtr);
00543
00544 if (newF >= oldF)
00545 {
00546 ratio = -1;
00547 }
00548 else
00549 {
00550
00551 rtype = oldSolnPtr->applyJacobian(*dirPtr, *bVecPtr);
00552 if (rtype != NOX::Abstract::Group::Ok)
00553 {
00554 utils->err() << "NOX::Solver::InexactTrustRegionBased::iterate - "
00555 << "unable to compute F" << endl;
00556 throw "NOX Error";
00557 }
00558
00559 double numerator = oldF - newF;
00560 double denominator = 0.0;
00561
00562 if (!Teuchos::is_null(meritFuncPtr)) {
00563 denominator = fabs(oldF - meritFuncPtr->
00564 computeQuadraticModel(dir, *oldSolnPtr));
00565 }
00566 else
00567 denominator = fabs(dir.innerProduct(oldSolnPtr->getGradient()) +
00568 0.5 * bVecPtr->innerProduct(*bVecPtr));
00569
00570 ratio = numerator / denominator;
00571 if (utils->isPrintType(NOX::Utils::InnerIteration))
00572 utils->out() << "Ratio computation: " << utils->sciformat(numerator) << "/"
00573 << utils->sciformat(denominator) << "=" << ratio << endl;
00574
00575
00576 if ((denominator < 1.0e-12) && ((newF / oldF) >= 0.5))
00577 ratio = -1;
00578 }
00579 }
00580
00581 if (utils->isPrintType(Utils::InnerIteration)) {
00582 utils->out() << "radius = " << utils->sciformat(radius, 1);
00583 utils->out() << " ratio = " << setprecision(1) << setw(3) << ratio;
00584 utils->out() << " f = " << utils->sciformat(sqrt(2*newF));
00585 utils->out() << " oldF = " << utils->sciformat(sqrt(2*oldF));
00586 utils->out() << " ";
00587
00588 switch(stepType) {
00589 case InexactTrustRegionBased::Newton:
00590 utils->out() << "Newton";
00591 break;
00592 case InexactTrustRegionBased::Cauchy:
00593 utils->out() << "Cauchy";
00594 break;
00595 case InexactTrustRegionBased::Dogleg:
00596 utils->out() << "Dogleg";
00597 break;
00598 }
00599
00600 utils->out() << endl;
00601 }
00602
00603
00604 if (ratio < contractTriggerRatio)
00605 {
00606 if (stepType == InexactTrustRegionBased::Newton) {
00607 radius = computeNorm(*newtonVecPtr);
00608 }
00609 radius = NOX_MAX(contractFactor * radius, minRadius);
00610 }
00611 else if ((ratio > expandTriggerRatio) && (dx == radius))
00612
00613
00614 {
00615 radius = NOX_MIN(expandFactor * radius, maxRadius);
00616 }
00617
00618 }
00619
00620
00621 if (useCounters) {
00622 if (stepType == InexactTrustRegionBased::Newton)
00623 numNewtonSteps += 1;
00624 else if (stepType == InexactTrustRegionBased::Cauchy)
00625 numCauchySteps += 1;
00626 else if (stepType == InexactTrustRegionBased::Dogleg) {
00627 numDoglegSteps += 1;
00628 sumDoglegFracCauchyToNewton += gamma;
00629 double tmp = radius/computeNorm(*newtonVecPtr);
00630 sumDoglegFracNewtonLength += tmp;
00631
00632 if (utils->isPrintType(Utils::Details)) {
00633 utils->out() << " Fraction of Newton Step Length = " << tmp << endl;
00634 utils->out() << " Fraction Between Cauchy and Newton Direction = "
00635 << gamma << endl;
00636 }
00637
00638 }
00639 }
00640
00641
00642 if ((radius <= minRadius) && (ratio < minRatio))
00643 {
00644 if (utils->isPrintType(Utils::InnerIteration))
00645 utils->out() << "Using recovery step and resetting trust region." << endl;
00646 solnPtr->computeX(*oldSolnPtr, *newtonVecPtr, recoveryStep);
00647 solnPtr->computeF();
00648 radius = computeNorm(*newtonVecPtr);
00649
00650
00651 }
00652
00653 status = test.checkStatus(*this, checkType);
00654
00655 if (utils->isPrintType(Utils::InnerIteration))
00656 utils->out() << NOX::Utils::fill(72) << endl;
00657
00658
00659 return status;
00660 }
00661
00662
00663
00664
00665 NOX::StatusTest::StatusType
00666 NOX::Solver::InexactTrustRegionBased::iterateInexact()
00667 {
00668
00669 if (status != StatusTest::Unconverged)
00670 return status;
00671
00672
00673 NOX::Abstract::Group& soln = *solnPtr;
00674 NOX::Abstract::Group& oldSoln = *oldSolnPtr;
00675 NOX::Abstract::Vector& newtonVec = *newtonVecPtr;
00676 NOX::Abstract::Vector& cauchyVec = *cauchyVecPtr;
00677 NOX::Abstract::Vector& rCauchyVec = *rCauchyVecPtr;
00678 NOX::Abstract::Vector& residualVec = *residualVecPtr;
00679 NOX::Abstract::Vector& aVec = *aVecPtr;
00680 NOX::Abstract::Vector& bVec = *bVecPtr;
00681 NOX::StatusTest::Generic& test = *testPtr;
00682
00683
00684 eta_last = eta;
00685 eta = inNewtonUtils.computeForcingTerm(soln, oldSoln, nIter,
00686 *this, eta_last);
00687
00688
00689 oldSoln = soln;
00690
00691
00692 bool ok = cauchyPtr->compute(cauchyVec, oldSoln, *this);
00693 if (!ok) {
00694 utils->out() << "NOX::Solver::InexactTrustRegionBased::iterate - "
00695 << "unable to calculate Cauchy direction" << endl;
00696 status = StatusTest::Failed;
00697 return status;
00698 }
00699
00700
00701
00702 oldF = meritFuncPtr->computef(oldSoln);
00703
00704
00705 oldSoln.applyJacobian(cauchyVec, rCauchyVec);
00706 rCauchyVec.update(1.0, oldSoln.getF(), 1.0);
00707 double rCauchy = computeNorm(rCauchyVec);
00708 double normF = computeNorm(oldSoln.getF());
00709 double etaCauchy = rCauchy/normF;
00710
00711
00712 double ratio = -1.0;
00713
00714
00715 if (nIter == 0)
00716 radius = maxRadius;
00717
00718 if (utils->isPrintType(NOX::Utils::InnerIteration))
00719 {
00720 utils->out() << NOX::Utils::fill(72) << endl;
00721 utils->out() << "-- Trust Region Inner Iteration --" << endl;
00722 }
00723
00724
00725 NOX::Abstract::Vector* dirPtr = 0;
00726 NOX::Abstract::Vector& doglegVec = bVec;
00727 NOX::Abstract::Vector& zVec = aVec;
00728 double step = 0.0;
00729 double tau = 0.0;
00730 double normS = computeNorm(cauchyVec);
00731 bool computedNewtonDir = false;
00732 innerIterationStatus = Unconverged;
00733
00734
00735 nIter ++;
00736
00737
00738 while (innerIterationStatus == Unconverged) {
00739
00740 if (useCounters)
00741 numTrustRegionInnerIterations += 1;
00742
00743
00744 if (computeNorm(cauchyVec) >= radius) {
00745 stepType = InexactTrustRegionBased::Cauchy;
00746 step = radius / computeNorm(cauchyVec);
00747 dirPtr = &cauchyVec;
00748 }
00749 else {
00750 if (etaCauchy <= eta) {
00751 stepType = InexactTrustRegionBased::Cauchy;
00752 step = 1.0;
00753 dirPtr = &cauchyVec;
00754 }
00755 else {
00756
00757 if (!computedNewtonDir) {
00758 string directionMethod = paramsPtr->sublist("Direction").
00759 get("Method", "Newton");
00760 Teuchos::ParameterList& lsParams = paramsPtr->sublist("Direction").
00761 sublist(directionMethod).sublist("Linear Solver");
00762 lsParams.set("Tolerance", eta);
00763 if (useCauchyInNewtonDirection)
00764 newtonVec.update(1.0, cauchyVec, 0.0);
00765 else
00766 newtonVec.init(0.0);
00767 if (!(oldSoln.isJacobian()))
00768 oldSoln.computeJacobian();
00769 oldSoln.applyJacobianInverse(lsParams, oldSoln.getF(), newtonVec);
00770 newtonVec.scale(-1.0);
00771 computedNewtonDir = true;
00772 }
00773
00774
00775 if ((computeNorm(newtonVec) <= radius) && (!useDoglegMinimization)) {
00776 stepType = InexactTrustRegionBased::Newton;
00777 step = 1.0;
00778 dirPtr = &newtonVec;
00779 }
00780 else {
00781
00782
00783 zVec.update(1.0, newtonVec, -1.0, cauchyVec, 0.0);
00784
00785
00786 double sDotZ = 0.0;
00787 sDotZ = cauchyVec.innerProduct(zVec);
00788
00789
00790 double normZ = computeNorm(zVec);
00791 if (sDotZ <= 0.0)
00792 tau = (-1.0 * sDotZ + sqrt(sDotZ * sDotZ + normZ * normZ *
00793 (radius * radius - normS * normS)))/(normZ * normZ);
00794 else
00795 tau = (radius * radius - normS * normS)/
00796 (sDotZ + sqrt(sDotZ * sDotZ + normZ * normZ *
00797 (radius * radius - normS * normS)));
00798
00799
00800 if (useDoglegMinimization) {
00801 double tauMin = 1.0;
00802
00803 oldSoln.applyJacobian(newtonVec, residualVec);
00804 residualVec.update(1.0, oldSoln.getF(), 1.0);
00805 residualVec.update(1.0, rCauchyVec, -1.0);
00806 double norm = computeNorm(residualVec);
00807 tauMin = rCauchyVec.innerProduct(residualVec) / (norm * norm);
00808
00809 tau = NOX_MIN(tauMin, tau);
00810 }
00811
00812
00813 doglegVec.update(1.0, cauchyVec, tau, zVec, 0.0);
00814 stepType = InexactTrustRegionBased::Dogleg;
00815 step = 1.0;
00816 dirPtr = &doglegVec;
00817 }
00818 }
00819 }
00820
00821
00822
00823 NOX::Abstract::Vector& dir = *dirPtr;
00824
00825
00826 soln.computeX(oldSoln, dir, step);
00827 NOX::Abstract::Group::ReturnType rtype = soln.computeF();
00828 if (rtype != NOX::Abstract::Group::Ok)
00829 throwError("iterateInexact", "unable to compute F!");
00830
00831
00832 newF = meritFuncPtr->computef(*solnPtr);
00833
00834 if (newF >= oldF)
00835 ratio = -1;
00836 else {
00837
00838
00839
00840 if (useAredPredRatio) {
00841
00842
00843 rtype = oldSoln.applyJacobian(*dirPtr, bVec);
00844 if (rtype != NOX::Abstract::Group::Ok)
00845 {
00846 utils->out() << "NOX::Solver::TrustRegionBased::iterate - "
00847 << "unable to compute F" << endl;
00848 throw "NOX Error";
00849 }
00850 bVec.update(1.0, oldSoln.getF(), 1.0);
00851
00852
00853 double oldNormF = computeNorm(oldSoln.getF());
00854 double newNormF = computeNorm(soln.getF());
00855 double normFLinear = computeNorm(bVec);;
00856
00857 ratio = (oldNormF - newNormF) / (oldNormF - normFLinear);
00858
00859
00860 if (utils->isPrintType(NOX::Utils::InnerIteration)) {
00861 double numerator = oldNormF - newNormF;
00862 double denominator = oldNormF - normFLinear;
00863 utils->out() << "Ratio computation: " << utils->sciformat(numerator) << "/"
00864 << utils->sciformat(denominator) << "=" << ratio << endl;
00865 }
00866
00867
00868 newF = meritFuncPtr->computef(*solnPtr);
00869
00870 }
00871 else {
00872
00873 rtype = oldSoln.applyJacobian(*dirPtr, zVec);
00874 if (rtype != NOX::Abstract::Group::Ok)
00875 throwError("iterateInexact", "unable to applyJacobian!");
00876
00877 double numerator = oldF - newF;
00878 double denominator = 0.0;
00879
00880 denominator = fabs(oldF - meritFuncPtr->
00881 computeQuadraticModel(dir,oldSoln));
00882
00883 ratio = numerator / denominator;
00884 if (utils->isPrintType(NOX::Utils::InnerIteration))
00885 utils->out() << "Ratio computation: " << utils->sciformat(numerator) << "/"
00886 << utils->sciformat(denominator) << "=" << ratio << endl;
00887
00888
00889 if ((denominator < 1.0e-12) && ((newF / oldF) >= 0.5))
00890 ratio = -1;
00891 }
00892 }
00893
00894 if (utils->isPrintType(Utils::InnerIteration)) {
00895 utils->out() << "radius = " << utils->sciformat(radius, 1);
00896 utils->out() << " ratio = " << setprecision(1) << setw(3) << ratio;
00897 utils->out() << " f = " << utils->sciformat(sqrt(2*newF));
00898 utils->out() << " oldF = " << utils->sciformat(sqrt(2*oldF));
00899 utils->out() << " ";
00900
00901 switch(stepType) {
00902 case InexactTrustRegionBased::Newton:
00903 utils->out() << "Newton";
00904 break;
00905 case InexactTrustRegionBased::Cauchy:
00906 utils->out() << "Cauchy";
00907 break;
00908 case InexactTrustRegionBased::Dogleg:
00909 utils->out() << "Dogleg";
00910 break;
00911 }
00912
00913 utils->out() << endl;
00914 }
00915
00916
00917 if (ratio >= minRatio) {
00918 innerIterationStatus = Converged;
00919 }
00920 else if ((radius <= minRadius) && (ratio < minRatio)) {
00921 innerIterationStatus = Failed;
00922 }
00923
00924
00925 if (ratio < contractTriggerRatio) {
00926
00927 if (stepType == InexactTrustRegionBased::Newton)
00928 radius = computeNorm(newtonVec);
00929
00930 radius = NOX_MAX(contractFactor * radius, minRadius);
00931 }
00932 else if (ratio > expandTriggerRatio)
00933 {
00934 radius = NOX_MIN(expandFactor * radius, maxRadius);
00935 }
00936
00937 }
00938
00939
00940 if (useCounters) {
00941 if (stepType == InexactTrustRegionBased::Newton)
00942 numNewtonSteps += 1;
00943 else if (stepType == InexactTrustRegionBased::Cauchy)
00944 numCauchySteps += 1;
00945 else if (stepType == InexactTrustRegionBased::Dogleg) {
00946 numDoglegSteps += 1;
00947 sumDoglegFracCauchyToNewton += tau;
00948 double tmp = radius/computeNorm(newtonVec);
00949 sumDoglegFracNewtonLength += tmp;
00950
00951 if (utils->isPrintType(Utils::Details)) {
00952 utils->out() << " Fraction of Newton Step Length = " << tmp << endl;
00953 utils->out() << " Fraction Between Cauchy and Newton Direction = "
00954 << tau << endl;
00955 }
00956
00957 }
00958 }
00959
00960
00961
00962
00963
00964 if (stepType == InexactTrustRegionBased::Cauchy)
00965 eta = 1.0 - (computeNorm(*dirPtr) / computeNorm(cauchyVec))*(1.0-etaCauchy);
00966 else if (stepType == InexactTrustRegionBased::Dogleg)
00967 eta = etaCauchy - tau * (etaCauchy - eta);
00968
00969
00970 if (innerIterationStatus == Failed) {
00971
00972 if (utils->isPrintType(Utils::InnerIteration))
00973 utils->out() << "Inner Iteration Failed!\n ";
00974 if (computedNewtonDir) {
00975 soln.computeX(oldSoln, cauchyVec, recoveryStep);
00976 utils->out() << "Using Newton recovery step and resetting trust region!" << endl;
00977 }
00978 else {
00979 soln.computeX(oldSoln, cauchyVec, recoveryStep);
00980 utils->out() << "Using Cauchy recovery step and resetting trust region!" << endl;
00981 }
00982 soln.computeF();
00983 radius = computeNorm(newtonVec);
00984 }
00985
00986
00987 status = test.checkStatus(*this, checkType);
00988
00989 if (utils->isPrintType(Utils::InnerIteration))
00990 utils->out() << NOX::Utils::fill(72) << endl;
00991
00992 return status;
00993 }
00994
00995
00996
00997
00998 NOX::StatusTest::StatusType
00999 NOX::Solver::InexactTrustRegionBased::
01000 checkStep(const NOX::Abstract::Vector& step,
01001 double& radius)
01002 {
01003 return NOX::StatusTest::Converged;
01004 }
01005
01006
01007
01008
01009 NOX::StatusTest::StatusType NOX::Solver::InexactTrustRegionBased::solve()
01010 {
01011 prePostOperator.runPreSolve(*this);
01012
01013
01014 while (status == StatusTest::Unconverged) {
01015 status = step();
01016 }
01017
01018 if (writeOutputParamsToList) {
01019 Teuchos::ParameterList& outputParams = paramsPtr->sublist("Output");
01020 outputParams.set("Nonlinear Iterations", nIter);
01021 outputParams.set("2-Norm of Residual", solnPtr->getNormF());
01022 if (useCounters) {
01023 Teuchos::ParameterList& trOutputParams = paramsPtr->
01024 sublist("Trust Region").sublist("Output");
01025 trOutputParams.set("Number of Cauchy Steps", numCauchySteps);
01026 trOutputParams.set("Number of Newton Steps", numNewtonSteps);
01027 trOutputParams.set("Number of Dogleg Steps", numDoglegSteps);
01028 trOutputParams.set("Number of Trust Region Inner Iterations",
01029 numTrustRegionInnerIterations);
01030 if (numDoglegSteps != 0) {
01031 trOutputParams.set("Dogleg Steps: Average Fraction of Newton Step Length", (sumDoglegFracNewtonLength/((double)numDoglegSteps)));
01032 trOutputParams.set("Dogleg Steps: Average Fraction Between Cauchy and Newton Direction", (sumDoglegFracCauchyToNewton/((double)numDoglegSteps)));
01033 }
01034
01035 }
01036 }
01037
01038 prePostOperator.runPostSolve(*this);
01039
01040 return status;
01041 }
01042
01043
01044
01045
01046 const Abstract::Group&
01047 NOX::Solver::InexactTrustRegionBased::getSolutionGroup() const
01048 {
01049 return *solnPtr;
01050 }
01051
01052
01053
01054
01055 const Abstract::Group&
01056 NOX::Solver::InexactTrustRegionBased::getPreviousSolutionGroup() const
01057 {
01058 return *oldSolnPtr;
01059 }
01060
01061
01062
01063
01064 int NOX::Solver::InexactTrustRegionBased::getNumIterations() const
01065 {
01066 return nIter;
01067 }
01068
01069
01070
01071
01072 const Teuchos::ParameterList&
01073 NOX::Solver::InexactTrustRegionBased::getList() const
01074 {
01075 return *paramsPtr;
01076 }
01077
01078
01079
01080
01081
01082 void NOX::Solver::InexactTrustRegionBased::printUpdate()
01083 {
01084
01085 if ((status == StatusTest::Unconverged) &&
01086 (utils->isPrintType(NOX::Utils::OuterIterationStatusTest))) {
01087 utils->out() << NOX::Utils::fill(72) << "\n";
01088 utils->out() << "-- Status Test Results --\n";
01089 testPtr->print(utils->out());
01090 utils->out() << NOX::Utils::fill(72) << "\n";
01091 }
01092
01093 double fmax = solnPtr->getF().norm(Abstract::Vector::MaxNorm);
01094 if (utils->isPrintType(NOX::Utils::OuterIteration)) {
01095 utils->out() << "\n" << NOX::Utils::fill(72) << "\n";
01096 utils->out() << "-- Newton Trust-Region Step " << nIter << " -- \n";
01097 utils->out() << "f = " << utils->sciformat(sqrt(2*newF));
01098 utils->out() << " fmax = " << utils->sciformat(fmax);
01099 utils->out() << " dx = " << utils->sciformat(dx);
01100 utils->out() << " radius = " << utils->sciformat(radius);
01101 if (status == StatusTest::Converged)
01102 utils->out() << " (Converged!)";
01103 if (status == StatusTest::Failed)
01104 utils->out() << " (Failed!)";
01105 utils->out() << "\n" << NOX::Utils::fill(72) << "\n" << endl;
01106 }
01107
01108 if ((status != StatusTest::Unconverged) &&
01109 (utils->isPrintType(NOX::Utils::OuterIteration))) {
01110 utils->out() << NOX::Utils::fill(72) << "\n";
01111 utils->out() << "-- Final Status Test Results --\n";
01112 testPtr->print(utils->out());
01113 utils->out() << NOX::Utils::fill(72) << "\n";
01114 }
01115 }
01116
01117
01118
01119
01120 void NOX::Solver::InexactTrustRegionBased::resetCounters()
01121 {
01122 numCauchySteps = 0;
01123 numNewtonSteps = 0;
01124 numDoglegSteps = 0;
01125 numTrustRegionInnerIterations = 0;
01126 sumDoglegFracCauchyToNewton = 0.0;
01127 sumDoglegFracNewtonLength = 0.0;
01128 return;
01129 }
01130
01131
01132
01133
01134 double NOX::Solver::InexactTrustRegionBased::
01135 computeNorm(const NOX::Abstract::Vector& v)
01136 {
01137 double norm = 0.0;
01138 norm = v.norm();
01139 return norm;
01140 }