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_LAPACK_Group.H"
00044 #include "NOX_Abstract_MultiVector.H"
00045 #include "Teuchos_BLAS_wrappers.hpp"
00046 #include "Teuchos_LAPACK_wrappers.hpp"
00047
00048 NOX::LAPACK::Group::Group(NOX::LAPACK::Interface& interface):
00049 xVector(interface.getInitialGuess()),
00050 fVector(xVector, ShapeCopy),
00051 newtonVector(xVector, ShapeCopy),
00052 gradientVector(xVector, ShapeCopy),
00053 jacSolver(xVector.length()),
00054 problemInterface(interface)
00055 {
00056 resetIsValid();
00057 }
00058
00059 NOX::LAPACK::Group::Group(const NOX::LAPACK::Group& source, NOX::CopyType type) :
00060 xVector(source.xVector, type),
00061 fVector(source.fVector, type),
00062 newtonVector(source.newtonVector, type),
00063 gradientVector(source.gradientVector, type),
00064 jacSolver(source.jacSolver),
00065 problemInterface(source.problemInterface)
00066 {
00067
00068 switch (type) {
00069
00070 case NOX::DeepCopy:
00071
00072 isValidF = source.isValidF;
00073 isValidGradient = source.isValidGradient;
00074 isValidNewton = source.isValidNewton;
00075 isValidJacobian = source.isValidJacobian;
00076 break;
00077
00078 case NOX::ShapeCopy:
00079 resetIsValid();
00080 break;
00081
00082 default:
00083 cerr << "NOX:LAPACK::Group - invalid CopyType for copy constructor." << endl;
00084 throw "NOX LAPACK Error";
00085 }
00086
00087 }
00088
00089 NOX::LAPACK::Group::~Group()
00090 {
00091 }
00092
00093 void NOX::LAPACK::Group::resetIsValid()
00094 {
00095 isValidF = false;
00096 isValidJacobian = false;
00097 isValidGradient = false;
00098 isValidNewton = false;
00099 jacSolver.reset();
00100 }
00101
00102 Teuchos::RCP<NOX::Abstract::Group> NOX::LAPACK::Group::
00103 clone(NOX::CopyType type) const
00104 {
00105 Teuchos::RCP<NOX::Abstract::Group> newgrp =
00106 Teuchos::rcp(new NOX::LAPACK::Group(*this, type));
00107 return newgrp;
00108 }
00109
00110 NOX::Abstract::Group& NOX::LAPACK::Group::operator=(const NOX::Abstract::Group& source)
00111 {
00112 return operator=(dynamic_cast<const Group&> (source));
00113 }
00114
00115 NOX::Abstract::Group& NOX::LAPACK::Group::operator=(const Group& source)
00116 {
00117 if (this != &source) {
00118
00119
00120 xVector = source.xVector;
00121
00122
00123 isValidF = source.isValidF;
00124 isValidGradient = source.isValidGradient;
00125 isValidNewton = source.isValidNewton;
00126 isValidJacobian = source.isValidJacobian;
00127
00128
00129 if (isValidF) {
00130 fVector = source.fVector;
00131 }
00132
00133 if (isValidGradient)
00134 gradientVector = source.gradientVector;
00135
00136 if (isValidNewton)
00137 newtonVector = source.newtonVector;
00138
00139 if (isValidJacobian)
00140 jacSolver = source.jacSolver;
00141
00142 }
00143
00144 return *this;
00145 }
00146
00147 void NOX::LAPACK::Group::setX(const NOX::Abstract::Vector& y)
00148 {
00149 setX(dynamic_cast<const Vector&> (y));
00150 }
00151
00152 void NOX::LAPACK::Group::setX(const Vector& y)
00153 {
00154 resetIsValid();
00155 xVector = y;
00156 }
00157
00158 void NOX::LAPACK::Group::computeX(const NOX::Abstract::Group& grp,
00159 const NOX::Abstract::Vector& d,
00160 double step)
00161 {
00162
00163 const Group& lapackgrp = dynamic_cast<const Group&> (grp);
00164 const Vector& lapackd = dynamic_cast<const Vector&> (d);
00165 computeX(lapackgrp, lapackd, step);
00166 }
00167
00168 void NOX::LAPACK::Group::computeX(const Group& grp, const Vector& d, double step)
00169 {
00170 resetIsValid();
00171 xVector.update(1.0, grp.xVector, step, d);
00172 }
00173
00174 NOX::Abstract::Group::ReturnType NOX::LAPACK::Group::computeF()
00175 {
00176 if (isValidF)
00177 return NOX::Abstract::Group::Ok;
00178
00179 isValidF = problemInterface.computeF(fVector, xVector);
00180
00181 return (isValidF) ? (NOX::Abstract::Group::Ok) : (NOX::Abstract::Group::Failed);
00182 }
00183
00184 NOX::Abstract::Group::ReturnType NOX::LAPACK::Group::computeJacobian()
00185 {
00186
00187 if (isValidJacobian)
00188 return NOX::Abstract::Group::Ok;
00189
00190 isValidJacobian =
00191 problemInterface.computeJacobian(jacSolver.getMatrix(), xVector);
00192
00193 return (isValidJacobian) ? (NOX::Abstract::Group::Ok) : (NOX::Abstract::Group::Failed);
00194 }
00195
00196 NOX::Abstract::Group::ReturnType NOX::LAPACK::Group::computeGradient()
00197 {
00198 if (isValidGradient)
00199 return NOX::Abstract::Group::Ok;
00200
00201 if (!isF()) {
00202 cerr << "ERROR: NOX::LAPACK::Group::computeGrad() - F is out of date wrt X!" << endl;
00203 return NOX::Abstract::Group::BadDependency;
00204 }
00205
00206 if (!isJacobian()) {
00207 cerr << "ERROR: NOX::LAPACK::Group::computeGrad() - Jacobian is out of date wrt X!" << endl;
00208 return NOX::Abstract::Group::BadDependency;
00209 }
00210
00211
00212 jacSolver.apply(true, 1, &fVector(0), &gradientVector(0));
00213
00214
00215 isValidGradient = true;
00216
00217
00218 return NOX::Abstract::Group::Ok;
00219 }
00220
00221 NOX::Abstract::Group::ReturnType NOX::LAPACK::Group::computeNewton(Teuchos::ParameterList& p)
00222 {
00223 if (isNewton())
00224 return NOX::Abstract::Group::Ok;
00225
00226 if (!isF()) {
00227 cerr << "ERROR: NOX::Example::Group::computeNewton() - invalid F" << endl;
00228 throw "NOX Error";
00229 }
00230
00231 if (!isJacobian()) {
00232 cerr << "ERROR: NOX::Example::Group::computeNewton() - invalid Jacobian" << endl;
00233 throw "NOX Error";
00234 }
00235
00236 NOX::Abstract::Group::ReturnType status = applyJacobianInverse(p, fVector, newtonVector);
00237 isValidNewton = (status == NOX::Abstract::Group::Ok);
00238
00239
00240 newtonVector.scale(-1.0);
00241
00242
00243 return status;
00244 }
00245
00246 NOX::Abstract::Group::ReturnType
00247 NOX::LAPACK::Group::applyJacobian(const Abstract::Vector& input,
00248 NOX::Abstract::Vector& result) const
00249 {
00250 const Vector& lapackinput = dynamic_cast<const Vector&> (input);
00251 Vector& lapackresult = dynamic_cast<Vector&> (result);
00252 return applyJacobian(lapackinput, lapackresult);
00253 }
00254
00255 NOX::Abstract::Group::ReturnType
00256 NOX::LAPACK::Group::applyJacobian(const Vector& input, Vector& result) const
00257 {
00258
00259 if (!isJacobian())
00260 return NOX::Abstract::Group::BadDependency;
00261
00262
00263 jacSolver.apply(false, 1, &input(0), &result(0));
00264
00265 return NOX::Abstract::Group::Ok;
00266 }
00267
00268 NOX::Abstract::Group::ReturnType
00269 NOX::LAPACK::Group::applyJacobianTranspose(const NOX::Abstract::Vector& input,
00270 NOX::Abstract::Vector& result) const
00271 {
00272 const Vector& lapackinput = dynamic_cast<const Vector&> (input);
00273 Vector& lapackresult = dynamic_cast<Vector&> (result);
00274 return applyJacobianTranspose(lapackinput, lapackresult);
00275 }
00276
00277 NOX::Abstract::Group::ReturnType
00278 NOX::LAPACK::Group::applyJacobianTranspose(const Vector& input, Vector& result) const
00279 {
00280
00281 if (!isJacobian())
00282 return NOX::Abstract::Group::BadDependency;
00283
00284
00285 jacSolver.apply(true, 1, &input(0), &result(0));
00286
00287 return NOX::Abstract::Group::Ok;
00288 }
00289
00290 NOX::Abstract::Group::ReturnType
00291 NOX::LAPACK::Group::applyJacobianInverse(Teuchos::ParameterList& p,
00292 const Abstract::Vector& input,
00293 NOX::Abstract::Vector& result) const
00294 {
00295 const Vector& lapackinput = dynamic_cast<const Vector&> (input);
00296 Vector& lapackresult = dynamic_cast<Vector&> (result);
00297 return applyJacobianInverse(p, lapackinput, lapackresult);
00298 }
00299
00300 NOX::Abstract::Group::ReturnType
00301 NOX::LAPACK::Group::applyJacobianInverse(Teuchos::ParameterList& p,
00302 const Vector& input,
00303 Vector& result) const
00304 {
00305
00306 if (!isJacobian()) {
00307 cerr << "ERROR: NOX::LAPACK::Group::applyJacobianInverse() - invalid Jacobian" << endl;
00308 throw "NOX Error";
00309 }
00310
00311
00312 result = input;
00313 bool res = jacSolver.solve(false, 1, &result(0));
00314
00315 return res ? (NOX::Abstract::Group::Ok) : (NOX::Abstract::Group::Failed);
00316 }
00317
00318 NOX::Abstract::Group::ReturnType
00319 NOX::LAPACK::Group::applyJacobianInverseMultiVector(
00320 Teuchos::ParameterList& p,
00321 const NOX::Abstract::MultiVector& input,
00322 NOX::Abstract::MultiVector& result) const
00323 {
00324
00325 if (!isJacobian()) {
00326 cerr << "ERROR: NOX::LAPACK::Group::applyJacobianInverseMultiVector() "
00327 << "- invalid Jacobian" << endl;
00328 throw "NOX Error";
00329 }
00330
00331
00332 int nVecs = input.numVectors();
00333
00334 int m = jacSolver.getMatrix().numRows();
00335
00336
00337 NOX::LAPACK::Matrix<double> B(m,nVecs);
00338 const NOX::LAPACK::Vector* constVecPtr;
00339 for (int j=0; j<nVecs; j++) {
00340 constVecPtr = dynamic_cast<const NOX::LAPACK::Vector*>(&(input[j]));
00341 for (int i=0; i<m; i++)
00342 B(i,j) = (*constVecPtr)(i);
00343 }
00344
00345
00346 bool res = jacSolver.solve(false, nVecs, &B(0,0));
00347
00348 if (!res)
00349 return NOX::Abstract::Group::Failed;
00350
00351
00352 NOX::LAPACK::Vector* vecPtr;
00353 for (int j=0; j<nVecs; j++) {
00354 vecPtr = dynamic_cast<NOX::LAPACK::Vector*>(&(result[j]));
00355 for (int i=0; i<m; i++)
00356 (*vecPtr)(i) = B(i,j);
00357 }
00358
00359 return NOX::Abstract::Group::Ok;
00360 }
00361
00362 bool NOX::LAPACK::Group::isF() const
00363 {
00364 return isValidF;
00365 }
00366
00367 bool NOX::LAPACK::Group::isJacobian() const
00368 {
00369 return isValidJacobian;
00370 }
00371
00372 bool NOX::LAPACK::Group::isGradient() const
00373 {
00374 return isValidGradient;
00375 }
00376
00377 bool NOX::LAPACK::Group::isNewton() const
00378 {
00379 return isValidNewton;
00380 }
00381
00382 const NOX::Abstract::Vector& NOX::LAPACK::Group::getX() const
00383 {
00384 return xVector;
00385 }
00386
00387 const NOX::Abstract::Vector& NOX::LAPACK::Group::getF() const
00388 {
00389 return fVector;
00390 }
00391
00392 double NOX::LAPACK::Group::getNormF() const
00393 {
00394 if (isValidF)
00395 return fVector.norm();
00396
00397 cerr << "ERROR: NOX::LAPACK::Group::getNormF() "
00398 << "- invalid F, please call computeF() first." << endl;
00399 throw "NOX Error";
00400
00401 return 0.0;
00402 }
00403
00404 const NOX::Abstract::Vector& NOX::LAPACK::Group::getGradient() const
00405 {
00406 return gradientVector;
00407 }
00408
00409 const NOX::Abstract::Vector& NOX::LAPACK::Group::getNewton() const
00410 {
00411 return newtonVector;
00412 }
00413
00414
00415 void NOX::LAPACK::Group::print() const
00416 {
00417 cout << "x = " << xVector << "\n";
00418
00419 if (isValidF) {
00420 cout << "F(x) = " << fVector << "\n";
00421 }
00422 else
00423 cout << "F(x) has not been computed" << "\n";
00424
00425 cout << endl;
00426 }