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 #include "NOX_LineSearch_NonlinearCG.H"
00040
00041 #include "NOX_Common.H"
00042 #include "NOX_Abstract_Vector.H"
00043 #include "NOX_Abstract_Group.H"
00044 #include "NOX_Solver_Generic.H"
00045 #include "Teuchos_ParameterList.hpp"
00046 #include "NOX_Utils.H"
00047 #include "NOX_GlobalData.H"
00048 #include "NOX_StatusTest_FiniteValue.H"
00049
00050 NOX::LineSearch::NonlinearCG::
00051 NonlinearCG(const Teuchos::RCP<NOX::GlobalData>& gd,
00052 Teuchos::ParameterList& params) :
00053 finiteValueTester( Teuchos::rcp(new StatusTest::FiniteValue()) )
00054 {
00055 reset(gd, params);
00056 }
00057
00058 NOX::LineSearch::NonlinearCG::~NonlinearCG()
00059 {
00060
00061 }
00062
00063 bool NOX::LineSearch::NonlinearCG::
00064 reset(const Teuchos::RCP<NOX::GlobalData>& gd,
00065 Teuchos::ParameterList& params)
00066 {
00067 utils = gd->getUtils();
00068
00069 return true;
00070 }
00071
00072 bool NOX::LineSearch::NonlinearCG::compute(Abstract::Group& newgrp,
00073 double& step,
00074 const Abstract::Vector& dir,
00075 const Solver::Generic& s)
00076 {
00077 if (utils->isPrintType(NOX::Utils::InnerIteration))
00078 {
00079 utils->out() << "\n" << NOX::Utils::fill(72) << "\n"
00080 << "-- NonlinearCG Line Search -- \n";
00081 }
00082
00083 const Abstract::Group& oldgrp = s.getPreviousSolutionGroup();
00084
00085
00086
00087
00088
00089
00090 double numerator = oldgrp.getF().innerProduct(dir);
00091 double denominator = computeDirectionalDerivative(dir, oldgrp).innerProduct(dir);
00092
00093 if( finiteValueTester->finiteNumberTest(step) )
00094 {
00095 utils->out() << "NOX::LineSearch::NonlinearCG::compute "
00096 << "- step value is NaN or Inf. " << endl;
00097 throw "NOX Error";
00098 }
00099
00100 step = - numerator / denominator;
00101 newgrp.computeX(oldgrp, dir, step);
00102 newgrp.computeF();
00103
00104 double checkOrthogonality = fabs( newgrp.getF().innerProduct(dir) );
00105
00106 if (utils->isPrintType(Utils::InnerIteration)) {
00107 utils->out() << setw(3) << "1" << ":";
00108 utils->out() << " step = " << utils->sciformat(step);
00109 utils->out() << " orth = " << utils->sciformat(checkOrthogonality);
00110 utils->out() << "\n" << NOX::Utils::fill(72) << "\n" << endl;
00111 }
00112
00113 return true;
00114 }
00115
00116
00117 NOX::Abstract::Vector& NOX::LineSearch::NonlinearCG::
00118 computeDirectionalDerivative(const Abstract::Vector& dir,
00119 const Abstract::Group& grp)
00120 {
00121
00122 if (Teuchos::is_null(vecPtr))
00123 vecPtr = dir.clone(ShapeCopy);
00124 if (Teuchos::is_null(grpPtr))
00125 grpPtr = grp.clone(ShapeCopy);
00126
00127
00128 if (!grp.isF())
00129 {
00130 utils->out() << "NOX::LineSearch::NonlinearCG::computeDirectionalDerivative "
00131 << "- Invalid F" << endl;
00132 throw "NOX Error";
00133 }
00134
00135
00136 double lambda = 1.0e-6;
00137 double denominator = dir.norm();
00138
00139
00140 if (denominator == 0.0)
00141 denominator = 1.0;
00142
00143 double eta = lambda * (lambda + grp.getX().norm() / denominator);
00144
00145
00146 if (eta == 0.0)
00147 eta = 1.0e-6;
00148
00149
00150 vecPtr->update(eta, dir, 1.0, grp.getX(), 0.0);
00151
00152
00153 grpPtr->setX(*vecPtr);
00154 grpPtr->computeF();
00155
00156
00157 vecPtr->update(-1.0/eta, grp.getF(), 1.0/eta, grpPtr->getF(), 0.0);
00158
00159 return(*vecPtr);
00160 }