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_TensorBased.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_Utils.H"
00048 #include "NOX_GlobalData.H"
00049 #include "NOX_Solver_SolverUtils.H"
00050
00051 #include "NOX_LineSearch_Utils_Printing.H"
00052 #include "NOX_LineSearch_Utils_Counters.H"
00053 #include "NOX_LineSearch_Utils_Slope.H"
00054
00055
00056 #define CHECK_RESIDUALS
00057 #define DEBUG_LEVEL 0
00058 #define DEVELOPER_CODE
00059 #define USE_INITIAL_GUESS_LOGIC
00060
00061 NOX::Solver::TensorBased::
00062 TensorBased(const Teuchos::RCP<NOX::Abstract::Group>& xGrp,
00063 const Teuchos::RCP<NOX::StatusTest::Generic>& t,
00064 const Teuchos::RCP<Teuchos::ParameterList>& p) :
00065 globalDataPtr(Teuchos::rcp(new NOX::GlobalData(p))),
00066 utilsPtr(globalDataPtr->getUtils()),
00067 solnPtr(xGrp),
00068 oldSolnPtr(xGrp->clone(DeepCopy)),
00069 newtonVecPtr(xGrp->getX().clone(ShapeCopy)),
00070 tensorVecPtr(xGrp->getX().clone(ShapeCopy)),
00071 aVecPtr(xGrp->getX().clone(ShapeCopy)),
00072 sVecPtr(xGrp->getX().clone(ShapeCopy)),
00073 tmpVecPtr(xGrp->getX().clone(ShapeCopy)),
00074 residualVecPtr(xGrp->getX().clone(ShapeCopy)),
00075 testPtr(t),
00076 paramsPtr(p),
00077 print(utilsPtr),
00078 slopeObj(globalDataPtr),
00079 prePostOperator(utilsPtr, paramsPtr->sublist("Solver Options"))
00080 {
00081 reset(xGrp, t, p);
00082 }
00083
00084
00085 void NOX::Solver::TensorBased::init()
00086 {
00087
00088 stepSize = 0;
00089 nIter = 0;
00090 status = NOX::StatusTest::Unconverged;
00091
00092
00093 counter.reset();
00094 numJvMults = 0;
00095 numJ2vMults = 0;
00096
00097
00098 if (utilsPtr->isPrintType(NOX::Utils::Parameters))
00099 {
00100 utilsPtr->out() << "\n" << NOX::Utils::fill(72) << "\n";
00101 utilsPtr->out() << "\n-- Parameters Passed to Nonlinear Solver --\n\n";
00102 paramsPtr->print(utilsPtr->out(),5);
00103 utilsPtr->out() << "\n" << NOX::Utils::fill(72) << "\n";
00104 }
00105
00106 }
00107
00108
00109 bool NOX::Solver::TensorBased::
00110 reset(const Teuchos::RCP<NOX::Abstract::Group>& xGrp,
00111 const Teuchos::RCP<NOX::StatusTest::Generic>& t,
00112 const Teuchos::RCP<Teuchos::ParameterList>& p)
00113 {
00114 solnPtr = xGrp;
00115 testPtr = t;
00116 paramsPtr = p;
00117
00118 globalDataPtr = Teuchos::rcp(new NOX::GlobalData(p));
00119 utilsPtr->reset(paramsPtr->sublist("Printing"));
00120 print.reset(utilsPtr);
00121 slopeObj.reset(globalDataPtr);
00122 prePostOperator.reset(utilsPtr, paramsPtr->sublist("Solver Options"));
00123
00124
00125 Teuchos::ParameterList& dirParams = paramsPtr->sublist("Direction");
00126
00127
00128 string choice = dirParams.get("Method", "Tensor");
00129 if (choice == "Tensor")
00130 requestedBaseStep = TensorStep;
00131 else if (choice == "Newton")
00132 requestedBaseStep = NewtonStep;
00133 else
00134 {
00135 if (utilsPtr->isPrintType(NOX::Utils::Error))
00136 utilsPtr->err() << "NOX::Direction::Tensor::reset() - The choice of "
00137 << "\"Method\" parameter \"" << choice
00138 << "\" is invalid." << endl;
00139 throw "NOX error";
00140 }
00141
00142
00143 Teuchos::ParameterList& teParams = dirParams.sublist(choice);
00144
00145
00146 dirParams.set("Compute Step", choice);
00147
00148
00149 doRescue = teParams.get("Rescue Bad Newton Solve", true);
00150
00151 checkType = parseStatusTestCheckType(paramsPtr->sublist("Solver Options"));
00152
00153
00154 useModifiedMethod = false;
00155 if (requestedBaseStep == TensorStep)
00156 {
00157 useModifiedMethod =
00158 dirParams.get("Use Modified Bouaricha", true);
00159 if (useModifiedMethod &&
00160 utilsPtr->isPrintType(NOX::Utils::Parameters))
00161 utilsPtr->out() << "Using Modifed Bouaricha method" << endl;
00162 }
00163
00164
00165
00166 Teuchos::ParameterList& lsParams = paramsPtr->sublist("Line Search");
00167
00168
00169 choice = lsParams.get("Method", "Curvilinear");
00170
00171 if (choice == "Curvilinear")
00172 lsType = Curvilinear;
00173 else if (choice == "Dual")
00174 lsType = Dual;
00175 else if (choice == "Standard")
00176 lsType = Standard;
00177 else if (choice == "Full Step")
00178 lsType = FullStep;
00179 else if (choice == "Newton")
00180 lsType = Newton;
00181 else
00182 {
00183 if (utilsPtr->isPrintType(NOX::Utils::Error))
00184 utilsPtr->err() << "NOX::Direction::Tensor::reset() - The choice of "
00185 << "\"Line Search\" parameter " << choice
00186 << " is invalid." << endl;
00187 throw "NOX Error";
00188 }
00189
00190 lsParams.set("Submethod", choice);
00191
00192
00193 Teuchos::ParameterList& gsParams = lsParams.sublist(choice);
00194
00195
00196 choice = gsParams.get("Recovery Step Type", "Constant");
00197 if (choice == "Constant")
00198 recoveryStepType = Constant;
00199 else if (choice == "Last Computed Step")
00200 recoveryStepType = LastComputedStep;
00201 else
00202 {
00203 utilsPtr->err() << "NOX::Solver::TensorBased::reset() - "
00204 << "Invalid \"Recovery Step Type\"" << endl;
00205 throw "NOX Error";
00206 }
00207
00208
00209 minStep = gsParams.get("Minimum Step", 1.0e-12);
00210 defaultStep = gsParams.get("Default Step", 1.0);
00211 recoveryStep = gsParams.get("Recovery Step", 0.0);
00212 maxIters = gsParams.get("Max Iters", 40);
00213 alpha = gsParams.get("Alpha Factor", 1.0e-4);
00214
00215 choice = gsParams.get("Lambda Selection", "Halving");
00216 if (choice == "Halving")
00217 lambdaSelection = Halving;
00218 else if (choice == "Quadratic")
00219 lambdaSelection = Quadratic;
00220 else
00221 {
00222 if (utilsPtr->isPrintType(NOX::Utils::Error))
00223 utilsPtr->err() << "NOX::Solver::TensorBased::reset() - The choice of "
00224 << "\"Lambda Selection\" parameter " << choice
00225 << " is invalid." << endl;
00226 throw "NOX Error";
00227 }
00228
00229 choice = gsParams.get("Sufficient Decrease Condition",
00230 "Armijo-Goldstein");
00231 if (choice == "Armijo-Goldstein")
00232 convCriteria = ArmijoGoldstein;
00233
00234
00235
00236
00237 else
00238 {
00239 if (utilsPtr->isPrintType(NOX::Utils::Error))
00240 utilsPtr->err() << "NOX::Solver::TensorBased::reset() - The choice of "
00241 << "\"Sufficient Decrease Condition\" parameter " << choice
00242 << " is invalid." << endl;
00243 throw "NOX Error";
00244 }
00245
00246 init();
00247 return true;
00248 }
00249
00250 void NOX::Solver::TensorBased::
00251 reset(const NOX::Abstract::Vector& initialGuess,
00252 const Teuchos::RCP<NOX::StatusTest::Generic>& t)
00253 {
00254 solnPtr->setX(initialGuess);
00255 testPtr = t;
00256 init();
00257 }
00258
00259 void NOX::Solver::TensorBased::
00260 reset(const NOX::Abstract::Vector& initialGuess)
00261 {
00262 solnPtr->setX(initialGuess);
00263 init();
00264 }
00265
00266 NOX::Solver::TensorBased::~TensorBased()
00267 {
00268 #ifdef DEVELOPER_CODE
00269 if (utilsPtr->isPrintType(NOX::Utils::Details))
00270 {
00271 utilsPtr->out() << "multsJv = " << numJvMults << " (linesearch)" << endl;
00272 utilsPtr->out() << "mults2Jv = " << numJ2vMults << endl;
00273 }
00274 #endif
00275 }
00276
00277
00278 NOX::StatusTest::StatusType NOX::Solver::TensorBased::getStatus()
00279 {
00280 return status;
00281 }
00282
00283 NOX::StatusTest::StatusType NOX::Solver::TensorBased::step()
00284 {
00285 prePostOperator.runPreIterate(*this);
00286
00287
00288 if (nIter ==0) {
00289
00290 NOX::Abstract::Group::ReturnType rtype = solnPtr->computeF();
00291 if (rtype != NOX::Abstract::Group::Ok) {
00292 utilsPtr->err() << "NOX::Solver::TensorBased::init - "
00293 << "Unable to compute F" << endl;
00294 throw "NOX Error";
00295 }
00296
00297
00298 status = testPtr->checkStatus(*this, checkType);
00299 if ((status == NOX::StatusTest::Converged) &&
00300 (utilsPtr->isPrintType(NOX::Utils::Warning))) {
00301 utilsPtr->out() << "Warning: NOX::Solver::TensorBased::init() - "
00302 << "The solution passed into the solver (either "
00303 << "through constructor or reset method) "
00304 << "is already converged! The solver will not "
00305 << "attempt to solve this system since status "
00306 << "is flagged as converged." << endl;
00307 }
00308
00309 printUpdate();
00310 }
00311
00312
00313 if (status != NOX::StatusTest::Unconverged)
00314 {
00315 prePostOperator.runPostIterate(*this);
00316 printUpdate();
00317 return status;
00318 }
00319
00320
00321 NOX::Abstract::Group& soln = *solnPtr;
00322 NOX::StatusTest::Generic& test = *testPtr;
00323
00324
00325 bool ok = computeTensorDirection(soln, *this);
00326 if (!ok)
00327 {
00328 if (utilsPtr->isPrintType(NOX::Utils::Error))
00329 utilsPtr->out() << "NOX::Solver::TensorBased::iterate - "
00330 << "unable to calculate direction" << endl;
00331 status = NOX::StatusTest::Failed;
00332 prePostOperator.runPostIterate(*this);
00333 printUpdate();
00334 return status;
00335 }
00336
00337
00338 nIter ++;
00339
00340
00341 *oldSolnPtr = soln;
00342
00343
00344 ok = implementGlobalStrategy(soln, stepSize, *this);
00345 if (!ok)
00346 {
00347 if (stepSize == 0.0)
00348 {
00349 if (utilsPtr->isPrintType(NOX::Utils::Error))
00350 utilsPtr->out() << "NOX::Solver::TensorBased::iterate - line search failed"
00351 << endl;
00352 status = NOX::StatusTest::Failed;
00353 prePostOperator.runPostIterate(*this);
00354 printUpdate();
00355 return status;
00356 }
00357 else if (utilsPtr->isPrintType(NOX::Utils::Warning))
00358 utilsPtr->out() << "NOX::Solver::TensorBased::iterate - "
00359 << "using recovery step for line search" << endl;
00360 }
00361
00362
00363 NOX::Abstract::Group::ReturnType rtype = soln.computeF();
00364 if (rtype != NOX::Abstract::Group::Ok)
00365 {
00366 if (utilsPtr->isPrintType(NOX::Utils::Error))
00367 utilsPtr->out() << "NOX::Solver::TensorBased::iterate - "
00368 << "unable to compute F" << endl;
00369 status = NOX::StatusTest::Failed;
00370 prePostOperator.runPostIterate(*this);
00371 printUpdate();
00372 return status;
00373 }
00374
00375 status = test.checkStatus(*this, checkType);
00376
00377 prePostOperator.runPostIterate(*this);
00378
00379 printUpdate();
00380
00381 return status;
00382 }
00383
00384
00385 NOX::StatusTest::StatusType NOX::Solver::TensorBased::solve()
00386 {
00387 prePostOperator.runPreSolve(*this);
00388
00389
00390 while (status == NOX::StatusTest::Unconverged)
00391 {
00392 status = step();
00393 }
00394
00395 Teuchos::ParameterList& outputParams = paramsPtr->sublist("Output");
00396 outputParams.set("Nonlinear Iterations", nIter);
00397 outputParams.set("2-Norm of Residual", solnPtr->getNormF());
00398
00399 prePostOperator.runPostSolve(*this);
00400
00401 return status;
00402 }
00403
00404 const NOX::Abstract::Group&
00405 NOX::Solver::TensorBased::getSolutionGroup() const
00406 {
00407 return *solnPtr;
00408 }
00409
00410 const NOX::Abstract::Group&
00411 NOX::Solver::TensorBased::getPreviousSolutionGroup() const
00412 {
00413 return *oldSolnPtr;
00414 }
00415
00416 int NOX::Solver::TensorBased::getNumIterations() const
00417 {
00418 return nIter;
00419 }
00420
00421 const Teuchos::ParameterList&
00422 NOX::Solver::TensorBased::getList() const
00423 {
00424 return *paramsPtr;
00425 }
00426
00427
00428 void NOX::Solver::TensorBased::printUpdate()
00429 {
00430 double normSoln = 0;
00431 double normStep = 0;
00432
00433
00434 if ((status == NOX::StatusTest::Unconverged) &&
00435 (utilsPtr->isPrintType(NOX::Utils::OuterIterationStatusTest)))
00436 {
00437 utilsPtr->out() << NOX::Utils::fill(72) << "\n";
00438 utilsPtr->out() << "-- Status Test Results --\n";
00439 testPtr->print(utilsPtr->out());
00440 utilsPtr->out() << NOX::Utils::fill(72) << "\n";
00441 }
00442
00443
00444 if (utilsPtr->isPrintType(NOX::Utils::OuterIteration))
00445 {
00446 normSoln = solnPtr->getNormF();
00447 normStep = (nIter > 0) ? tensorVecPtr->norm() : 0;
00448 }
00449
00450
00451 if (utilsPtr->isPrintType(NOX::Utils::OuterIteration))
00452 {
00453 utilsPtr->out() << "\n" << NOX::Utils::fill(72) << "\n";
00454 utilsPtr->out() << "-- Nonlinear Solver Step " << nIter << " -- \n";
00455 utilsPtr->out() << "f = " << utilsPtr->sciformat(normSoln);
00456 utilsPtr->out() << " step = " << utilsPtr->sciformat(stepSize);
00457 utilsPtr->out() << " dx = " << utilsPtr->sciformat(normStep);
00458 if (status == NOX::StatusTest::Converged)
00459 utilsPtr->out() << " (Converged!)";
00460 if (status == NOX::StatusTest::Failed)
00461 utilsPtr->out() << " (Failed!)";
00462 utilsPtr->out() << "\n" << NOX::Utils::fill(72) << "\n" << endl;
00463 }
00464
00465
00466 if ((status != NOX::StatusTest::Unconverged) &&
00467 (utilsPtr->isPrintType(NOX::Utils::OuterIteration)))
00468 {
00469 utilsPtr->out() << NOX::Utils::fill(72) << "\n";
00470 utilsPtr->out() << "-- Final Status Test Results --\n";
00471 testPtr->print(utilsPtr->out());
00472 utilsPtr->out() << NOX::Utils::fill(72) << "\n";
00473 }
00474 }
00475
00476
00477 bool
00478 NOX::Solver::TensorBased::computeTensorDirection(NOX::Abstract::Group& soln,
00479 const NOX::Solver::Generic& solver)
00480 {
00481 NOX::Abstract::Group::ReturnType dir_status;
00482
00483 Teuchos::ParameterList& linearParams = paramsPtr->sublist("Direction").
00484 sublist(paramsPtr->sublist("Direction").
00485 get("Method","Tensor")).
00486 sublist("Linear Solver");
00487
00488
00489 dir_status = soln.computeF();
00490 if (dir_status != NOX::Abstract::Group::Ok)
00491 throwError("computeTensorDirection", "Unable to compute F");
00492
00493
00494 dir_status = soln.computeJacobian();
00495 if (dir_status != NOX::Abstract::Group::Ok)
00496 throwError("computeTensorDirection", "Unable to compute Jacobian");
00497
00498
00499 double sDotS = 0.0;
00500 int tempVal1 = 0;
00501 if ((nIter > 0) && (requestedBaseStep == TensorStep))
00502 {
00503
00504 *sVecPtr = soln.getX();
00505 sVecPtr->update(1.0, solver.getPreviousSolutionGroup().getX(), -1.0);
00506 double normS = sVecPtr->norm();
00507 sDotS = normS * normS;
00508
00509
00510 soln.applyJacobian(*sVecPtr, *aVecPtr);
00511 numJvMults++;
00512 aVecPtr->update(1.0, solver.getPreviousSolutionGroup().getF(), -1.0);
00513 aVecPtr->update(-1.0, soln.getF(), 1.0);
00514 if (sDotS != 0)
00515 aVecPtr->scale(1.0 / (sDotS * sDotS));
00516
00517
00518 *tmpVecPtr = *newtonVecPtr;
00519 tmpVecPtr->scale(-1.0);
00520
00521
00522 soln.applyJacobian(*tmpVecPtr, *residualVecPtr);
00523 numJvMults++;
00524 residualVecPtr->update(1.0, solver.getPreviousSolutionGroup().getF(),-1.0);
00525 double residualNorm = residualVecPtr->norm();
00526
00527 #if DEBUG_LEVEL > 0
00528 double tmpVecNorm = tmpVecPtr->norm();
00529 double residualNormRel = residualNorm /
00530 solver.getPreviousSolutionGroup().getNormF();
00531 if (utilsPtr->isPrintType(NOX::Utils::Details))
00532 {
00533 utilsPtr->out() << " Norm of initial guess: " << utilsPtr->sciformat(tmpVecNorm, 6)
00534 << endl;
00535 utilsPtr->out() << " initg norm of model residual = "
00536 << utilsPtr->sciformat(residualNorm, 6) << " (abs) "
00537 << utilsPtr->sciformat(residualNormRel, 6) << " (rel)" << endl;
00538 }
00539 #endif
00540
00541
00542 double tol = linearParams.get("Tolerance", 1e-4);
00543 double relativeResidual = residualNorm /
00544 solver.getPreviousSolutionGroup().getNormF();
00545
00546
00547 bool isInitialGuessGood = false;
00548 #ifdef USE_INITIAL_GUESS_LOGIC
00549 if (relativeResidual < 1.0)
00550 {
00551 if (utilsPtr->isPrintType(NOX::Utils::Details))
00552 utilsPtr->out() << " Initial guess is good..." << endl;
00553 isInitialGuessGood = true;
00554
00555 *tensorVecPtr = *tmpVecPtr;
00556 double newTol = tol / relativeResidual;
00557 if (newTol > 0.99)
00558 newTol = 0.99;
00559 linearParams.set("Tolerance", newTol);
00560 if (utilsPtr->isPrintType(NOX::Utils::Details))
00561 utilsPtr->out() << " Setting tolerance to " << utilsPtr->sciformat(newTol,6) << endl;
00562 }
00563 else
00564 #endif // USE_INITIAL_GUESS_LOGIC
00565 {
00566
00567 isInitialGuessGood = false;
00568 *residualVecPtr = solver.getPreviousSolutionGroup().getF();
00569 }
00570
00571
00572 tmpVecPtr->init(0.0);
00573 dir_status = soln.applyJacobianInverse(linearParams, *residualVecPtr,
00574 *tmpVecPtr);
00575
00576
00577 if (dir_status != NOX::Abstract::Group::Ok)
00578 {
00579 if (doRescue == false)
00580 throwError("computeTensorDirection", "Unable to apply Jacobian inverse");
00581 else if ((doRescue == true) &&
00582 (utilsPtr->isPrintType(NOX::Utils::Warning)))
00583 utilsPtr->out() << "WARNING: NOX::Solver::TensorBased::computeTensorDirection() - "
00584 << "Linear solve failed to achieve convergence - "
00585 << "using the step anyway "
00586 << "since \"Rescue Bad Newton Solve\" is true." << endl;
00587 }
00588
00589
00590 #ifdef USE_INITIAL_GUESS_LOGIC
00591 if (isInitialGuessGood)
00592 {
00593 tmpVecPtr->update(1.0, *tensorVecPtr, 1.0);
00594 linearParams.set("Tolerance", tol);
00595 }
00596 #endif
00597
00598
00599 if (linearParams.sublist("Output").
00600 isParameter("Number of Linear Iterations"))
00601 tempVal1 = linearParams.sublist("Output").
00602 get("Number of Linear Iterations",0);
00603
00604 #if DEBUG_LEVEL > 0
00605
00606 soln.applyJacobian(*tmpVecPtr, *residualVecPtr);
00607 numJvMults++;
00608 residualVec.update(-1.0, solver.getPreviousSolutionGroup().getF(),1.0);
00609 double residualNorm2 = residualVec.norm();
00610 double residualNorm2Rel = residualNorm2 /
00611 solver.getPreviousSolutionGroup().getNormF();
00612 if (utilsPtr->isPrintType(NOX::Utils::Details))
00613 utilsPtr->out() << " jifp norm of model residual = "
00614 << utilsPtr->sciformat(residualNorm2, 6) << " (abs) "
00615 << utilsPtr->sciformat(residualNorm2Rel, 6) << " (rel)" << endl;
00616 #endif
00617 }
00618
00619
00620 dir_status = soln.computeNewton(linearParams);
00621
00622
00623 if (dir_status != NOX::Abstract::Group::Ok)
00624 {
00625 if (doRescue == false)
00626 throwError("computeTensorDirection", "Unable to apply Jacobian inverse");
00627 else if ((doRescue == true) &&
00628 (utilsPtr->isPrintType(NOX::Utils::Warning)))
00629 utilsPtr->out() << "WARNING: NOX::Solver::TensorBased::computeTensorDirection() - "
00630 << "Linear solve failed to achieve convergence - "
00631 << "using the step anyway "
00632 << "since \"Rescue Bad Newton Solve\" is true." << endl;
00633 }
00634
00635
00636 *newtonVecPtr = soln.getNewton();
00637
00638
00639 int tempVal2 = 0;
00640 if (linearParams.sublist("Output").
00641 isParameter("Number of Linear Iterations"))
00642 tempVal2 = linearParams.sublist("Output").
00643 get("Number of Linear Iterations",0);
00644 numJ2vMults += (tempVal1 > tempVal2) ? tempVal1 : tempVal2;
00645
00646 #ifdef CHECK_RESIDUALS
00647 printDirectionInfo("newtonVec", *newtonVecPtr, soln, false);
00648 #endif // CHECK_RESIDUALS
00649
00650
00651 if ((nIter > 0) && (requestedBaseStep == TensorStep))
00652 {
00653
00654
00655
00656 tmpVecPtr->update(1.0, *newtonVecPtr, 1.0);
00657 tmpVecPtr->update(-1.0, *sVecPtr, 1.0);
00658 if (sDotS != 0.0)
00659 tmpVecPtr->scale( 1.0 / (sDotS * sDotS));
00660
00661
00662 sTinvJF = -sVecPtr->innerProduct(*newtonVecPtr);
00663 sTinvJa = sVecPtr->innerProduct(*tmpVecPtr);
00664 double qval = 0;
00665 double lambdaBar = 1;
00666 beta = calculateBeta(sTinvJa, 1.0, sTinvJF, qval, lambdaBar);
00667
00668 double sVecNorm = sVecPtr->norm();
00669 double aVecNorm = aVecPtr->norm();
00670 if (utilsPtr->isPrintType(NOX::Utils::Details))
00671 {
00672 utilsPtr->out() << " sTinvJF = " << utilsPtr->sciformat(sTinvJF, 6)
00673 << " sTinvJa = " << utilsPtr->sciformat(sTinvJa, 6) << endl;
00674 utilsPtr->out() << " norm(s) = " << utilsPtr->sciformat(sVecNorm, 6)
00675 << " norm(a) = " << utilsPtr->sciformat(aVecNorm, 6) << endl;
00676 }
00677
00678 if (useModifiedMethod)
00679 {
00680 double alpha2 = lambdaBar;
00681 if (utilsPtr->isPrintType(NOX::Utils::Details))
00682 utilsPtr->out() << " Beta = " << utilsPtr->sciformat(beta, 6)
00683 << " Alpha2 = " << utilsPtr->sciformat(alpha2, 6) << endl;
00684 if (alpha2 != 1.0)
00685 {
00686 if (utilsPtr->isPrintType(NOX::Utils::Details))
00687 utilsPtr->out() << " *** Scaling tensor term a ***" << endl;
00688 aVecPtr->scale(alpha2);
00689 tmpVecPtr->scale(alpha2);
00690 sTinvJa *= alpha2;
00691 beta /= alpha2;
00692 lambdaBar = 1.0;
00693 qval = 0;
00694 }
00695 }
00696
00697
00698 tensorVecPtr->update(1.0, *newtonVecPtr, -beta*beta, *tmpVecPtr, 0.0);
00699
00700 #ifdef CHECK_RESIDUALS
00701 printDirectionInfo("tensorVec", *tensorVecPtr, soln, true);
00702 #endif // CHECK_RESIDUALS
00703 #if DEBUG_LEVEL > 0
00704 double sDotT = tensorVecPtr->innerProduct(sVec);
00705 if (utilsPtr->isPrintType(NOX::Utils::Details))
00706 utilsPtr->out() << " Beta = " << utilsPtr->sciformat(beta, 6)
00707 << " std = " << utilsPtr->sciformat(sDotT, 6)
00708 << " qval = " << utilsPtr->sciformat(qval, 2)
00709 << " lambdaBar = " << lambdaBar << endl;
00710 #endif
00711 }
00712 else
00713 *tensorVecPtr = *newtonVecPtr;
00714
00715 return true;
00716 }
00717
00718
00719 double NOX::Solver::TensorBased::calculateBeta(double qa,
00720 double qb,
00721 double qc,
00722 double& qval,
00723 double& lambdaBar,
00724 double lambda) const
00725 {
00726 double beta_value = 0.0;
00727 double discriminant = qb*qb - 4*qa*qc*lambda;
00728
00729 if (discriminant < 0.0)
00730 {
00731
00732 beta_value = -qb / qa / 2.0;
00733 qval = (qa * beta_value * beta_value) + (qb * beta_value) + (lambda * qc);
00734 lambdaBar = qb*qb / (4*qa*qc);
00735 #if DEBUG_LEVEL > 0
00736 if (utilsPtr->isPrintType(NOX::Utils::Details))
00737 utilsPtr->out() << " #### LambdaBar = " << lambdaBar << " ####\n";
00738 #endif
00739 }
00740 else
00741 {
00742 qval = 0;
00743 lambdaBar = 1.0;
00744 if ( (fabs(qa / qb) < 1e-8) && (fabs(lambda * qc / qb) < 1) )
00745 {
00746 #if DEBUG_LEVEL > 0
00747 if (utilsPtr->isPrintType(NOX::Utils::Details))
00748 utilsPtr->out() << " qa is relatively small\n";
00749 #endif
00750 beta_value = -lambda * qc / qb;
00751 }
00752 else
00753 {
00754 double tmp1 = (-qb + sqrt(discriminant)) / (2*qa);
00755 double tmp2 = (-qb - sqrt(discriminant)) / (2*qa);
00756 beta_value = (fabs(tmp1) < fabs(tmp2)) ? tmp1 : tmp2;
00757 #if DEBUG_LEVEL > 1
00758 if (utilsPtr->isPrintType(NOX::Utils::Details))
00759 utilsPtr->out() << " tmp1 = " << utilsPtr->sciformat(tmp1, 6)
00760 << " tmp2 = " << utilsPtr->sciformat(tmp2, 6)
00761 << " dir0xsc = " << utilsPtr->sciformat(dir0xsc, 6)
00762 << " normS = " << utilsPtr->sciformat(normS, 6)
00763 << endl;
00764 #endif
00765 }
00766 }
00767 #if DEBUG_LEVEL > 1
00768 if (utilsPtr->isPrintType(NOX::Utils::Details))
00769 utilsPtr->out() << " qa,qb,qc = " << utilsPtr->sciformat(qa, 6)
00770 << utilsPtr->sciformat(qb, 6)
00771 << utilsPtr->sciformat(qc, 6)
00772 << " beta = " << utilsPtr->sciformat(beta_value, 6)
00773 << endl;
00774 #endif
00775
00776 return beta_value;
00777 }
00778
00779
00780 bool
00781 NOX::Solver::TensorBased::computeCurvilinearStep(NOX::Abstract::Vector& dir,
00782 const NOX::Abstract::Group& soln,
00783 const NOX::Solver::Generic& s,
00784 double& lambda)
00785 {
00786 double qval = 0;
00787 double lambdaBar = 1;
00788 double beta1 = calculateBeta(sTinvJa, 1, sTinvJF, qval, lambdaBar, lambda);
00789 double betaFactor = ( (beta == 0.0) ? 0.0 : beta1*beta1 / (beta*beta));
00790
00791 dir.update(lambda - betaFactor, *newtonVecPtr, betaFactor, *tensorVecPtr, 0.0);
00792
00793 #if DEBUG_LEVEL > 0
00794 double sDotD = dir.innerProduct(sVec);
00795 if (utilsPtr->isPrintType(NOX::Utils::Details))
00796 {
00797 utilsPtr->out() << " Beta = " << utilsPtr->sciformat(beta, 6)
00798 << " std = " << utilsPtr->sciformat(sDotD, 6)
00799 << " qval = " << qval
00800 << " lambdaBar = " << lambdaBar
00801 << endl;
00802 utilsPtr->out() << " betaFactor = " << utilsPtr->sciformat(betaFactor,6)
00803 << " beta1 = " << utilsPtr->sciformat(beta1, 6)
00804 << endl;
00805 }
00806 #endif
00807
00808 return true;
00809 }
00810
00811
00812 bool
00813 NOX::Solver::TensorBased::implementGlobalStrategy(NOX::Abstract::Group& newGrp,
00814 double& in_stepSize,
00815 const NOX::Solver::Generic& s)
00816 {
00817 bool ok;
00818 counter.incrementNumLineSearches();
00819 isNewtonDirection = false;
00820 NOX::Abstract::Vector& searchDirection = *tensorVecPtr;
00821
00822 if ((counter.getNumLineSearches() == 1) || (lsType == Newton))
00823 {
00824 isNewtonDirection = true;
00825 searchDirection = *newtonVecPtr;
00826 }
00827
00828
00829 if ((lsType != Dual) || (isNewtonDirection))
00830 ok = performLinesearch(newGrp, in_stepSize, searchDirection, s);
00831 else if (lsType == Dual)
00832 {
00833 double fTensor = 0.0;
00834 double fNew = 0.0;
00835 double tensorStep = 1.0;
00836 bool isTensorDescent = false;
00837
00838 const Abstract::Group& oldGrp = s.getPreviousSolutionGroup();
00839 double fprime = slopeObj.computeSlope(searchDirection, oldGrp);
00840
00841
00842 if (fprime < 0)
00843 {
00844 ok = performLinesearch(newGrp, in_stepSize, searchDirection, s);
00845 fTensor = 0.5 * newGrp.getNormF() * newGrp.getNormF();
00846 tensorStep = in_stepSize;
00847 isTensorDescent = true;
00848 }
00849
00850
00851 ok = performLinesearch(newGrp, in_stepSize, *newtonVecPtr, s);
00852 fNew = 0.5 * newGrp.getNormF() * newGrp.getNormF();
00853
00854
00855 if (isTensorDescent && (fTensor <= fNew))
00856 {
00857 newGrp.computeX(oldGrp, *tensorVecPtr, tensorStep);
00858 newGrp.computeF();
00859 }
00860 }
00861
00862 return ok;
00863 }
00864
00865
00866 bool
00867 NOX::Solver::TensorBased::performLinesearch(NOX::Abstract::Group& newSoln,
00868 double& in_stepSize,
00869 const NOX::Abstract::Vector& lsDir,
00870 const NOX::Solver::Generic& s)
00871 {
00872 if (print.isPrintType(NOX::Utils::InnerIteration))
00873 {
00874 utilsPtr->out() << "\n" << NOX::Utils::fill(72) << "\n";
00875 utilsPtr->out() << "-- Tensor Line Search (";
00876 if (lsType == Curvilinear)
00877 utilsPtr->out() << "Curvilinear";
00878 else if (lsType == Standard)
00879 utilsPtr->out() << "Standard";
00880 else if (lsType == FullStep)
00881 utilsPtr->out() << "Full Step";
00882 else if (lsType == Dual)
00883 utilsPtr->out() << "Dual";
00884 utilsPtr->out() << ") -- " << endl;
00885 }
00886
00887
00888 bool isFailed = false;
00889 bool isAcceptable = false;
00890 bool isFirstPass = true;
00891 string message = "(STEP ACCEPTED!)";
00892
00893
00894 int lsIterations = 1;
00895
00896
00897 const Abstract::Group& oldSoln = s.getPreviousSolutionGroup();
00898 double fOld = 0.5 * oldSoln.getNormF() * oldSoln.getNormF();
00899
00900
00901 in_stepSize = defaultStep;
00902 newSoln.computeX(oldSoln, lsDir, in_stepSize);
00903 newSoln.computeF();
00904 double fNew = 0.5 * newSoln.getNormF() * newSoln.getNormF();
00905
00906
00907 if (lsType == FullStep)
00908 {
00909 print.printStep(lsIterations, in_stepSize, fOld, fNew, message);
00910 return (!isFailed);
00911 }
00912
00913
00914 double fprime;
00915 if ((lsType == Curvilinear) && !(isNewtonDirection))
00916 fprime = slopeObj.computeSlope(*newtonVecPtr, oldSoln);
00917 else
00918 fprime = slopeObj.computeSlope(lsDir, oldSoln);
00919 numJvMults++;
00920
00921
00922 double threshold = fOld + alpha*in_stepSize*fprime;
00923 isAcceptable = (fNew < threshold);
00924
00925
00926 if (!isAcceptable)
00927 {
00928 counter.incrementNumNonTrivialLineSearches();
00929 *tmpVecPtr = lsDir;
00930 }
00931
00932
00933 while (!isAcceptable)
00934 {
00935
00936 if (lsIterations > maxIters)
00937 {
00938 isFailed = true;
00939 message = "(FAILED - Max Iters)";
00940 break;
00941 }
00942
00943 print.printStep(lsIterations, in_stepSize, fOld, fNew);
00944
00945
00946 if (isFirstPass &&
00947 (!isNewtonDirection) &&
00948 (fprime >= 0) &&
00949 (lsType != Curvilinear) )
00950 {
00951 *tmpVecPtr = *newtonVecPtr;
00952 fprime = slopeObj.computeSlope(*tmpVecPtr, oldSoln);
00953 numJvMults++;
00954
00955 if (utilsPtr->isPrintType(NOX::Utils::Details))
00956 utilsPtr->out() << " Switching to Newton step. New fprime = "
00957 << utilsPtr->sciformat(fprime, 6) << endl;
00958 }
00959 else
00960 {
00961 in_stepSize = selectLambda(fNew, fOld, fprime, in_stepSize);
00962 }
00963
00964 isFirstPass = false;
00965
00966
00967 if (in_stepSize < minStep)
00968 {
00969 isFailed = true;
00970 message = "(FAILED - Min Step)";
00971 break;
00972 }
00973
00974
00975 counter.incrementNumIterations();
00976 lsIterations ++;
00977
00978
00979 if ((lsType == Curvilinear) && !(isNewtonDirection))
00980 {
00981 computeCurvilinearStep(*tmpVecPtr, oldSoln, s, in_stepSize);
00982
00983 newSoln.computeX(oldSoln, *tmpVecPtr, 1.0);
00984 }
00985 else
00986 {
00987 newSoln.computeX(oldSoln, *tmpVecPtr, in_stepSize);
00988 }
00989 newSoln.computeF();
00990 fNew = 0.5 * newSoln.getNormF() * newSoln.getNormF();
00991
00992
00993 threshold = fOld + alpha*in_stepSize*fprime;
00994 isAcceptable = (fNew < threshold);
00995 }
00996
00997
00998 if (isFailed)
00999 {
01000 counter.incrementNumFailedLineSearches();
01001
01002 if (recoveryStepType == Constant)
01003 {
01004 in_stepSize = recoveryStep;
01005 if (in_stepSize == 0.0)
01006 {
01007 newSoln = oldSoln;
01008 newSoln.computeF();
01009 fNew = fOld;
01010 }
01011 else
01012 {
01013
01014 if ((lsType == Curvilinear) && !(isNewtonDirection))
01015 {
01016 computeCurvilinearStep(*tmpVecPtr, oldSoln, s, in_stepSize);
01017
01018 newSoln.computeX(oldSoln, *tmpVecPtr, 1.0);
01019 }
01020 else
01021 {
01022 newSoln.computeX(oldSoln, *tmpVecPtr, in_stepSize);
01023 }
01024
01025 newSoln.computeF();
01026 fNew = 0.5 * newSoln.getNormF() * newSoln.getNormF();
01027 message = "(USING RECOVERY STEP!)";
01028 }
01029 }
01030 else
01031 message = "(USING LAST STEP!)";
01032 }
01033
01034 print.printStep(lsIterations, in_stepSize, fOld, fNew, message);
01035 counter.setValues(paramsPtr->sublist("Line Search"));
01036
01037 return (!isFailed);
01038 }
01039
01040
01041 double
01042 NOX::Solver::TensorBased::getNormModelResidual(
01043 const NOX::Abstract::Vector& dir,
01044 const NOX::Abstract::Group& soln,
01045 bool isTensorModel) const
01046 {
01047
01048
01049 Teuchos::RCP<NOX::Abstract::Vector> residualPtr =
01050 soln.getF().clone(ShapeCopy);
01051 soln.applyJacobian(dir, *residualPtr);
01052 numJvMults++;
01053 residualPtr->update(1.0, soln.getF(), 1.0);
01054
01055
01056 if (isTensorModel)
01057 {
01058 double tmp = sVecPtr->innerProduct(dir);
01059 if (utilsPtr->isPrintType(NOX::Utils::Details))
01060 utilsPtr->out() << " sc'*dt = " << utilsPtr->sciformat(tmp, 6) << endl;
01061 residualPtr->update(tmp*tmp, *aVecPtr, 1.0);
01062 }
01063
01064 double modelNorm = residualPtr->norm();
01065 return modelNorm;
01066 }
01067
01068
01069 void
01070 NOX::Solver::TensorBased::printDirectionInfo(std::string dirName,
01071 const NOX::Abstract::Vector& dir,
01072 const NOX::Abstract::Group& soln,
01073 bool isTensorModel) const
01074 {
01075 double dirNorm = dir.norm();
01076
01077 double residual = getNormModelResidual(dir, soln, isTensorModel);
01078 double residualRel = residual / soln.getNormF();
01079
01080 double fprime = getDirectionalDerivative(dir, soln);
01081 double fprimeRel = fprime / dirNorm;
01082
01083 if (utilsPtr->isPrintType(NOX::Utils::Details))
01084 {
01085 utilsPtr->out() << " " << dirName << " norm of model residual = "
01086 << utilsPtr->sciformat(residual, 6) << " (abs) "
01087 << utilsPtr->sciformat(residualRel, 6) << " (rel)" << endl;
01088 utilsPtr->out() << " " << dirName << " directional derivative = "
01089 << utilsPtr->sciformat(fprime, 6) << " (abs) "
01090 << utilsPtr->sciformat(fprimeRel, 6) << " (rel)" << endl;
01091 utilsPtr->out() << " " << dirName << " norm = "
01092 << utilsPtr->sciformat(dirNorm, 6) << endl;
01093 }
01094 }
01095
01096
01097 double NOX::Solver::TensorBased::getDirectionalDerivative(
01098 const NOX::Abstract::Vector& dir,
01099 const NOX::Abstract::Group& soln) const
01100 {
01101 Teuchos::RCP<NOX::Abstract::Vector> tmpPtr =
01102 soln.getF().clone(ShapeCopy);
01103 soln.applyJacobian(dir, *tmpPtr);
01104 numJvMults++;
01105 double fprime = tmpPtr->innerProduct(soln.getF());
01106 return fprime;
01107 }
01108
01109
01110 double NOX::Solver::TensorBased::selectLambda(double fNew, double fOld,
01111 double fOldPrime,
01112 double lambda)
01113 {
01114 double lambdaRet;
01115 double temp;
01116
01117 if (lambdaSelection == Quadratic)
01118 {
01119 temp = -fOldPrime / (2.0*(fNew - fOld - fOldPrime));
01120 if (temp < 0.1)
01121 temp = 0.1;
01122 lambdaRet = temp * lambda;
01123 }
01124 else
01125 {
01126 lambdaRet = 0.5 * lambda;
01127 }
01128 return lambdaRet;
01129 }
01130
01131
01132 void NOX::Solver::TensorBased::throwError(const string& functionName,
01133 const string& errorMsg) const
01134 {
01135 if (utilsPtr->isPrintType(NOX::Utils::Error))
01136 utilsPtr->err() << "NOX::Solver::TensorBased::" << functionName
01137 << " - " << errorMsg << endl;
01138 throw "NOX Error";
01139 }
01140