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_MeritFunction_SumOfSquares.H"
00044 #include "NOX_Abstract_Vector.H"
00045 #include "NOX_Abstract_Group.H"
00046
00047 NOX::MeritFunction::SumOfSquares::
00048 SumOfSquares(const Teuchos::RCP<NOX::Utils>& u) :
00049 meritFunctionName("Sum Of Squares (default): 0.5 * ||F|| * ||F||")
00050 {
00051 utils = u;
00052 }
00053
00054 NOX::MeritFunction::SumOfSquares::~SumOfSquares()
00055 {
00056 }
00057
00058 double NOX::MeritFunction::SumOfSquares::
00059 computef(const NOX::Abstract::Group& grp) const
00060 {
00061 if ( !(grp.isF()) ) {
00062 utils->err()
00063 << "ERROR: NOX::MeritFunction::SumOfSquares::computef() - "
00064 << "F has not been computed yet!. Please call "
00065 << "computeF() on the group passed into this function."
00066 << endl;
00067 throw "NOX Error";
00068 }
00069
00070 return (0.5 * grp.getNormF() * grp.getNormF());
00071 }
00072
00073 void NOX::MeritFunction::SumOfSquares::
00074 computeGradient(const NOX::Abstract::Group& grp,
00075 NOX::Abstract::Vector& result) const
00076 {
00077 if ( !(grp.isF()) ) {
00078 utils->err()
00079 << "ERROR: NOX::MeritFunction::SumOfSquares::computeGradient() - "
00080 << "F has not been computed yet!. Please call "
00081 << "computeF() on the group passed into this function."
00082 << endl;
00083 throw "NOX Error";
00084 }
00085
00086 if ( !(grp.isJacobian()) ) {
00087 utils->err()
00088 << "ERROR: NOX::MeritFunction::SumOfSquares::computeGradient() - "
00089 << "Jacobian has not been computed yet!. Please call "
00090 << "computeJacobian() on the group passed into this function."
00091 << endl;
00092 throw "NOX Error";
00093 }
00094
00095 NOX::Abstract::Group::ReturnType status =
00096 grp.applyJacobianTranspose(grp.getF(), result);
00097
00098 if (status != NOX::Abstract::Group::Ok) {
00099 utils->err() << "ERROR: NOX::MeritFunction::SumOfSquares::compute"
00100 << "Gradient - applyJacobianTranspose failed!" << endl;
00101 throw "NOX Error";
00102 }
00103
00104 return;
00105 }
00106
00107 double NOX::MeritFunction::SumOfSquares::
00108 computeSlope(const NOX::Abstract::Vector& dir,
00109 const NOX::Abstract::Group& grp) const
00110 {
00111 if (Teuchos::is_null(tmpVecPtr))
00112 tmpVecPtr = grp.getF().clone();
00113
00114
00115
00116 if (!(grp.isJacobian()))
00117 return this->computeSlopeWithoutJacobian(dir, grp);
00118
00119 this->computeGradient(grp, *(tmpVecPtr.get()));
00120
00121 return dir.innerProduct(*(tmpVecPtr.get()));
00122 }
00123
00124 double NOX::MeritFunction::SumOfSquares::
00125 computeQuadraticModel(const NOX::Abstract::Vector& dir,
00126 const NOX::Abstract::Group& grp) const
00127 {
00128 if (Teuchos::is_null(tmpVecPtr))
00129 tmpVecPtr = grp.getF().clone();
00130
00131 double m = 0.0;
00132
00133 m = this->computef(grp);
00134
00135 m += this->computeSlope(dir, grp);
00136
00137 grp.applyJacobian(dir, *(tmpVecPtr.get()));
00138
00139 m += 0.5 * tmpVecPtr->innerProduct(*(tmpVecPtr.get()));
00140
00141 return m;
00142 }
00143
00144 double NOX::MeritFunction::SumOfSquares::
00145 computeSlopeWithoutJacobian(const NOX::Abstract::Vector& dir,
00146 const NOX::Abstract::Group& grp) const
00147 {
00148 if (Teuchos::is_null(tmpVecPtr))
00149 tmpVecPtr = grp.getF().clone(NOX::ShapeCopy);
00150
00151 if (Teuchos::is_null(tmpGrpPtr))
00152 tmpGrpPtr = grp.clone(NOX::ShapeCopy);
00153
00154
00155
00156 double lambda = 1.0e-6;
00157 double denominator = dir.norm();
00158
00159
00160 if (denominator == 0.0)
00161 denominator = 1.0;
00162
00163 double eta = lambda * (lambda + grp.getX().norm() / denominator);
00164
00165
00166 if (eta == 0.0)
00167 eta = 1.0e-6;
00168
00169
00170 tmpVecPtr->update(eta, dir, 1.0, grp.getX(), 0.0);
00171
00172
00173 tmpGrpPtr->setX(*(tmpVecPtr.get()));
00174 tmpGrpPtr->computeF();
00175
00176
00177 tmpVecPtr->update(-1.0/eta, grp.getF(), 1.0/eta, tmpGrpPtr->getF(), 0.0);
00178
00179 return(tmpVecPtr->innerProduct(grp.getF()));
00180 }
00181
00182 void NOX::MeritFunction::SumOfSquares::
00183 computeQuadraticMinimizer(const NOX::Abstract::Group& grp,
00184 NOX::Abstract::Vector& result) const
00185 {
00186
00187 if (Teuchos::is_null(tmpVecPtr))
00188 tmpVecPtr = grp.getF().clone(NOX::ShapeCopy);
00189
00190
00191 if ( !(grp.isF()) ) {
00192 utils->err()
00193 << "ERROR: NOX::MeritFunction::SumOfSquares::"
00194 << "computeQuadraticMinimizer() - "
00195 << "F has not been computed yet!. Please call "
00196 << "computeF() on the group passed into this function."
00197 << endl;
00198 throw "NOX Error";
00199 }
00200
00201 if ( !(grp.isJacobian()) ) {
00202 utils->err()
00203 << "ERROR: NOX::MeritFunction::SumOfSquares::"
00204 << "computeQuadraticMinimizer() - "
00205 << "Jacobian has not been computed yet!. Please call "
00206 << "computeJacobian() on the group passed into this function."
00207 << endl;
00208 throw "NOX Error";
00209 }
00210
00211
00212 this->computeGradient(grp, result);
00213
00214
00215 NOX::Abstract::Group::ReturnType status =
00216 grp.applyJacobian(result, *tmpVecPtr);
00217 if (status != NOX::Abstract::Group::Ok) {
00218 utils->err()
00219 << "ERROR: NOX::MeritFunction::SumOfSquares::"
00220 << "computeQuadraticMinimizer() - grp->applyJacobian() has failed!"
00221 << endl;
00222 throw "NOX Error";
00223 }
00224
00225 result.scale( -1.0 * result.innerProduct(result) /
00226 tmpVecPtr->innerProduct(*tmpVecPtr) );
00227
00228 }
00229
00230 const string& NOX::MeritFunction::SumOfSquares::name() const
00231 {
00232 return meritFunctionName;
00233 }