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 "Teuchos_TestForException.hpp"
00043 #include "Thyra_ModelEvaluator.hpp"
00044 #include "Thyra_SolveSupportTypes.hpp"
00045 #include "NOX_Common.H"
00046 #include "NOX_Thyra_Group.H"
00047 #include "NOX_Abstract_MultiVector.H"
00048 #include "NOX_Thyra_MultiVector.H"
00049
00050 NOX::Thyra::Group::
00051 Group(const NOX::Thyra::Vector& initial_guess,
00052 const Teuchos::RCP< const ::Thyra::ModelEvaluator<double> >& model):
00053 model_(model)
00054 {
00055 x_vec_ = Teuchos::rcp(new NOX::Thyra::Vector(initial_guess, DeepCopy));
00056 f_vec_ = Teuchos::rcp(new NOX::Thyra::Vector(initial_guess, ShapeCopy));
00057 newton_vec_ = Teuchos::rcp(new NOX::Thyra::Vector(initial_guess, ShapeCopy));
00058 gradient_vec_ =
00059 Teuchos::rcp(new NOX::Thyra::Vector(initial_guess, ShapeCopy));
00060
00061 shared_jacobian_ = Teuchos::rcp(new NOX::SharedObject< ::Thyra::LinearOpWithSolveBase<double, double>, NOX::Thyra::Group >(model->create_W()));
00062
00063 in_args_ = model_->createInArgs();
00064 out_args_ = model_->createOutArgs();
00065
00066 resetIsValidFlags();
00067 }
00068
00069 NOX::Thyra::Group::Group(const NOX::Thyra::Group& source, NOX::CopyType type) :
00070 model_(source.model_),
00071 shared_jacobian_(source.shared_jacobian_)
00072 {
00073
00074 x_vec_ = Teuchos::rcp(new NOX::Thyra::Vector(*source.x_vec_, type));
00075 f_vec_ = Teuchos::rcp(new NOX::Thyra::Vector(*source.f_vec_, type));
00076 newton_vec_ =
00077 Teuchos::rcp(new NOX::Thyra::Vector(*source.newton_vec_, type));
00078 gradient_vec_ =
00079 Teuchos::rcp(new NOX::Thyra::Vector(*source.gradient_vec_, type));
00080
00081 in_args_ = model_->createInArgs();
00082 out_args_ = model_->createOutArgs();
00083
00084 if (type == NOX::DeepCopy) {
00085 is_valid_f_ = source.is_valid_f_;
00086 is_valid_jacobian_ = source.is_valid_jacobian_;
00087 is_valid_newton_dir_ = source.is_valid_newton_dir_;
00088 is_valid_gradient_dir_ = source.is_valid_gradient_dir_;
00089
00090
00091 if (this->isJacobian())
00092 shared_jacobian_->getObject(this);
00093 }
00094 else if (type == NOX::ShapeCopy) {
00095 resetIsValidFlags();
00096 }
00097 else {
00098 TEST_FOR_EXCEPTION(true, std::logic_error,
00099 "NOX Error - Copy type is invalid!");
00100 }
00101
00102
00103 }
00104
00105 NOX::Thyra::Group::~Group()
00106 { }
00107
00108 void NOX::Thyra::Group::resetIsValidFlags()
00109 {
00110 is_valid_f_ = false;
00111 is_valid_jacobian_ = false;
00112 is_valid_newton_dir_ = false;
00113 is_valid_gradient_dir_ = false;
00114 }
00115
00116 Teuchos::RCP<NOX::Abstract::Group> NOX::Thyra::Group::
00117 clone(NOX::CopyType type) const
00118 {
00119 Teuchos::RCP<NOX::Abstract::Group> newgrp =
00120 Teuchos::rcp(new NOX::Thyra::Group(*this, type));
00121 return newgrp;
00122 }
00123
00124 NOX::Abstract::Group& NOX::Thyra::Group::operator=(const NOX::Abstract::Group& source)
00125 {
00126 return operator=(dynamic_cast<const NOX::Thyra::Group&> (source));
00127 }
00128
00129 NOX::Abstract::Group& NOX::Thyra::Group::operator=(const Group& source)
00130 {
00131
00132
00133 *x_vec_ = *source.x_vec_;
00134
00135 is_valid_f_ = source.is_valid_f_;
00136 is_valid_jacobian_ = source.is_valid_jacobian_;
00137 is_valid_newton_dir_ = source.is_valid_newton_dir_;
00138 is_valid_gradient_dir_ = source.is_valid_gradient_dir_;
00139
00140 if (this->isF())
00141 *f_vec_ = *(source.f_vec_);
00142
00143
00144 shared_jacobian_ = source.shared_jacobian_;
00145
00146 if (this->isNewton())
00147 *newton_vec_ = *(source.newton_vec_);
00148
00149 if (this->isGradient())
00150 *gradient_vec_ = *(source.gradient_vec_);
00151
00152
00153 if (this->isJacobian())
00154 shared_jacobian_->getObject(this);
00155
00156 return *this;
00157 }
00158
00159
00160 Teuchos::RCP<const ::Thyra::VectorBase<double> >
00161 NOX::Thyra::Group::get_current_x() const
00162 {
00163 if (is_null(x_vec_))
00164 return Teuchos::null;
00165 return x_vec_->getThyraRCPVector();
00166 }
00167
00168
00169 Teuchos::RCP< ::Thyra::LinearOpWithSolveBase<double> >
00170 NOX::Thyra::Group::getNonconstJacobian()
00171 {
00172 return shared_jacobian_->getObject(this);
00173 }
00174
00175
00176 Teuchos::RCP<const ::Thyra::LinearOpWithSolveBase<double> >
00177 NOX::Thyra::Group::getJacobian() const
00178 {
00179 return shared_jacobian_->getObject();
00180 }
00181
00182
00183 void NOX::Thyra::Group::setX(const NOX::Abstract::Vector& y)
00184 {
00185 setX(dynamic_cast<const NOX::Thyra::Vector&> (y));
00186 }
00187
00188
00189 void NOX::Thyra::Group::setX(const NOX::Thyra::Vector& y)
00190 {
00191 resetIsValidFlags();
00192 *x_vec_ = y;
00193 }
00194
00195
00196 void NOX::Thyra::Group::computeX(const NOX::Abstract::Group& grp,
00197 const NOX::Abstract::Vector& d,
00198 double step)
00199 {
00200 const NOX::Thyra::Group& thyra_grp =
00201 dynamic_cast<const NOX::Thyra::Group&> (grp);
00202
00203 const NOX::Thyra::Vector& thyra_d =
00204 dynamic_cast<const NOX::Thyra::Vector&> (d);
00205
00206 this->computeX(thyra_grp, thyra_d, step);
00207 }
00208
00209 void NOX::Thyra::Group::computeX(const NOX::Thyra::Group& grp,
00210 const NOX::Thyra::Vector& d,
00211 double step)
00212 {
00213 this->resetIsValidFlags();
00214 x_vec_->update(1.0, *(grp.x_vec_), step, d);
00215 }
00216
00217 NOX::Abstract::Group::ReturnType NOX::Thyra::Group::computeF()
00218 {
00219 if (this->isF())
00220 return NOX::Abstract::Group::Ok;
00221
00222 in_args_.set_x(x_vec_->getThyraRCPVector().assert_not_null());
00223 out_args_.set_f(f_vec_->getThyraRCPVector().assert_not_null());
00224
00225 model_->evalModel(in_args_, out_args_);
00226 in_args_.set_x(Teuchos::null);
00227 out_args_.set_f(Teuchos::null);
00228
00229 is_valid_f_ = true;
00230
00231 if (out_args_.isFailed())
00232 return NOX::Abstract::Group::Failed;
00233
00234 return NOX::Abstract::Group::Ok;
00235 }
00236
00237 NOX::Abstract::Group::ReturnType NOX::Thyra::Group::computeJacobian()
00238 {
00239 if (this->isJacobian())
00240 return NOX::Abstract::Group::Ok;
00241
00242 in_args_.set_x(x_vec_->getThyraRCPVector());
00243 out_args_.set_W(shared_jacobian_->getObject(this));
00244 model_->evalModel(in_args_, out_args_);
00245 in_args_.set_x(Teuchos::null);
00246 out_args_.set_W(Teuchos::null);
00247
00248 is_valid_jacobian_ = true;
00249
00250 if (out_args_.isFailed())
00251 return NOX::Abstract::Group::Failed;
00252
00253 return NOX::Abstract::Group::Ok;
00254 }
00255
00256 NOX::Abstract::Group::ReturnType NOX::Thyra::Group::computeGradient()
00257 {
00258 if ( ::Thyra::opSupported(*shared_jacobian_->getObject(), ::Thyra::TRANS) ) {
00259 TEST_FOR_EXCEPTION(true, std::logic_error,
00260 "NOX Error - compute gradient not implemented yet!");
00261 return NOX::Abstract::Group::Ok;
00262 }
00263 return NOX::Abstract::Group::Failed;
00264 }
00265
00266 NOX::Abstract::Group::ReturnType NOX::Thyra::Group::
00267 computeNewton(Teuchos::ParameterList& p)
00268 {
00269 NOX::Abstract::Group::ReturnType status =
00270 this->applyJacobianInverse(p, *f_vec_, *newton_vec_);
00271 newton_vec_->scale(-1.0);
00272
00273 return status;
00274 }
00275
00276 NOX::Abstract::Group::ReturnType
00277 NOX::Thyra::Group::applyJacobian(const Abstract::Vector& input,
00278 NOX::Abstract::Vector& result) const
00279 {
00280 using Teuchos::dyn_cast;
00281 return applyJacobian(dyn_cast<const NOX::Thyra::Vector>(input),
00282 dyn_cast<NOX::Thyra::Vector>(result));
00283
00284 }
00285
00286 NOX::Abstract::Group::ReturnType
00287 NOX::Thyra::Group::applyJacobian(const Vector& input, Vector& result) const
00288 {
00289 if ( !(this->isJacobian()) ) {
00290 TEST_FOR_EXCEPTION(true, std::logic_error,
00291 "NOX Error - Jacobian is not valid. " <<
00292 "Call computeJacobian before calling applyJacobian!");
00293 }
00294
00295 ::Thyra::apply(*shared_jacobian_->getObject(), ::Thyra::NOTRANS,
00296 input.getThyraVector(), &result.getThyraVector());
00297
00298 return NOX::Abstract::Group::Ok;
00299 }
00300
00301 NOX::Abstract::Group::ReturnType
00302 NOX::Thyra::Group::applyJacobianMultiVector(
00303 const NOX::Abstract::MultiVector& input,
00304 NOX::Abstract::MultiVector& result) const
00305 {
00306 if ( !(this->isJacobian()) ) {
00307 TEST_FOR_EXCEPTION(true, std::logic_error,
00308 "NOX Error - Jacobian is not valid. " <<
00309 "Call computeJacobian before calling applyJacobian!");
00310 }
00311
00312 const NOX::Thyra::MultiVector& nt_input =
00313 Teuchos::dyn_cast<const NOX::Thyra::MultiVector>(input);
00314 NOX::Thyra::MultiVector& nt_result =
00315 Teuchos::dyn_cast<NOX::Thyra::MultiVector>(result);
00316
00317 ::Thyra::apply(*shared_jacobian_->getObject(),
00318 ::Thyra::NOTRANS,
00319 *nt_input.getThyraMultiVector(),
00320 nt_result.getThyraMultiVector().get());
00321
00322 return NOX::Abstract::Group::Ok;
00323 }
00324
00325 NOX::Abstract::Group::ReturnType
00326 NOX::Thyra::Group::applyJacobianTranspose(const NOX::Abstract::Vector& input,
00327 NOX::Abstract::Vector& result) const
00328 {
00329 using Teuchos::dyn_cast;
00330 return applyJacobianTranspose(dyn_cast<const NOX::Thyra::Vector>(input),
00331 dyn_cast<NOX::Thyra::Vector>(result));
00332 }
00333
00334 NOX::Abstract::Group::ReturnType
00335 NOX::Thyra::Group::applyJacobianTranspose(const NOX::Thyra::Vector& input,
00336 NOX::Thyra::Vector& result) const
00337 {
00338 if ( !(this->isJacobian()) ) {
00339 TEST_FOR_EXCEPTION(true, std::logic_error,
00340 "NOX Error - Jacobian is not valid. " <<
00341 "Call computeJacobian before calling applyJacobian!");
00342 }
00343
00344 if ( ::Thyra::opSupported(*shared_jacobian_->getObject(), ::Thyra::TRANS) ) {
00345 ::Thyra::apply(*shared_jacobian_->getObject(), ::Thyra::TRANS,
00346 input.getThyraVector(), &result.getThyraVector());
00347 return NOX::Abstract::Group::Ok;
00348 }
00349 return NOX::Abstract::Group::Failed;
00350 }
00351
00352 NOX::Abstract::Group::ReturnType
00353 NOX::Thyra::Group::applyJacobianTransposeMultiVector(
00354 const NOX::Abstract::MultiVector& input,
00355 NOX::Abstract::MultiVector& result) const
00356 {
00357 if ( !(this->isJacobian()) ) {
00358 TEST_FOR_EXCEPTION(true, std::logic_error,
00359 "NOX Error - Jacobian is not valid. " <<
00360 "Call computeJacobian before calling applyJacobian!");
00361 }
00362
00363 if (! ::Thyra::opSupported(*shared_jacobian_->getObject(), ::Thyra::TRANS) )
00364 return NOX::Abstract::Group::Failed;
00365
00366 const NOX::Thyra::MultiVector& nt_input =
00367 Teuchos::dyn_cast<const NOX::Thyra::MultiVector>(input);
00368 NOX::Thyra::MultiVector& nt_result =
00369 Teuchos::dyn_cast<NOX::Thyra::MultiVector>(result);
00370
00371 ::Thyra::apply(*shared_jacobian_->getObject(),
00372 ::Thyra::TRANS,
00373 *nt_input.getThyraMultiVector(),
00374 nt_result.getThyraMultiVector().get());
00375
00376 return NOX::Abstract::Group::Ok;
00377 }
00378
00379 NOX::Abstract::Group::ReturnType
00380 NOX::Thyra::Group::applyJacobianInverse(Teuchos::ParameterList& p,
00381 const Abstract::Vector& input,
00382 NOX::Abstract::Vector& result) const
00383 {
00384 using Teuchos::dyn_cast;
00385 return applyJacobianInverse(p, dyn_cast<const NOX::Thyra::Vector>(input),
00386 dyn_cast<NOX::Thyra::Vector>(result));
00387 }
00388
00389 NOX::Abstract::Group::ReturnType
00390 NOX::Thyra::Group::applyJacobianInverse(Teuchos::ParameterList& p,
00391 const NOX::Thyra::Vector& input,
00392 NOX::Thyra::Vector& result) const
00393 {
00394 return applyJacobianInverseMultiVector( p, input.getThyraVector(),
00395 result.getThyraVector() );
00396 }
00397
00398 NOX::Abstract::Group::ReturnType NOX::Thyra::Group::
00399 applyJacobianInverseMultiVector(Teuchos::ParameterList& p,
00400 const NOX::Abstract::MultiVector& input,
00401 NOX::Abstract::MultiVector& result) const
00402 {
00403 const NOX::Thyra::MultiVector& nt_input =
00404 Teuchos::dyn_cast<const NOX::Thyra::MultiVector>(input);
00405 NOX::Thyra::MultiVector& nt_result =
00406 Teuchos::dyn_cast<NOX::Thyra::MultiVector>(result);
00407
00408 return applyJacobianInverseMultiVector(p,
00409 *nt_input.getThyraMultiVector(),
00410 *nt_result.getThyraMultiVector());
00411 }
00412
00413 bool NOX::Thyra::Group::isF() const
00414 {
00415 return is_valid_f_;
00416 }
00417
00418 bool NOX::Thyra::Group::isJacobian() const
00419 {
00420 return ((shared_jacobian_->isOwner(this)) && (is_valid_jacobian_));
00421 }
00422
00423 bool NOX::Thyra::Group::isNewton() const
00424 {
00425 return is_valid_newton_dir_;
00426 }
00427
00428 bool NOX::Thyra::Group::isGradient() const
00429 {
00430 return is_valid_gradient_dir_;
00431 }
00432
00433 const NOX::Abstract::Vector& NOX::Thyra::Group::getX() const
00434 {
00435 return *x_vec_;
00436 }
00437
00438 const NOX::Abstract::Vector& NOX::Thyra::Group::getF() const
00439 {
00440 return *f_vec_;
00441 }
00442
00443 double NOX::Thyra::Group::getNormF() const
00444 {
00445 if ( this->isF() )
00446 return f_vec_->norm();
00447
00448 cerr << "ERROR: NOX::Thyra::Group::getNormF() "
00449 << "- F is not up to date. Please call computeF()!" << endl;
00450 throw "NOX Error";
00451
00452 return 0.0;
00453 }
00454
00455 const NOX::Abstract::Vector& NOX::Thyra::Group::getNewton() const
00456 {
00457 return *newton_vec_;
00458 }
00459
00460 const NOX::Abstract::Vector& NOX::Thyra::Group::getGradient() const
00461 {
00462 return *gradient_vec_;
00463 }
00464
00465
00466 void NOX::Thyra::Group::print() const
00467 {
00468 TEST_FOR_EXCEPT(true);
00469 }
00470
00471
00472 NOX::Abstract::Group::ReturnType NOX::Thyra::Group::
00473 applyJacobianInverseMultiVector(Teuchos::ParameterList& p,
00474 const ::Thyra::MultiVectorBase<double>& input,
00475 ::Thyra::MultiVectorBase<double>& result) const
00476 {
00477 ::Thyra::SolveCriteria<double> solveCriteria;
00478 solveCriteria.requestedTol = p.get("Tolerance", 1.0e-6);
00479
00480 std::string numer_measure = p.get("Solve Measure Numerator",
00481 "Norm Residual");
00482 std::string denom_measure = p.get("Solve Measure Denominator",
00483 "Norm Initial Residual");
00484 solveCriteria.solveMeasureType =
00485 ::Thyra::SolveMeasureType(getThyraNormType(numer_measure),
00486 getThyraNormType(denom_measure));
00487
00488
00489 ::Thyra::assign(&result, 0.0);
00490
00491 const ::Thyra::SolveStatus<double> solve_status =
00492 ::Thyra::solve(*shared_jacobian_->getObject(),
00493 ::Thyra::NOTRANS, input, &result,
00494 &solveCriteria);
00495
00496
00497
00498
00499 if (solve_status.solveStatus == ::Thyra::SOLVE_STATUS_CONVERGED)
00500 return NOX::Abstract::Group::Ok;
00501 else if (solve_status.solveStatus == ::Thyra::SOLVE_STATUS_UNCONVERGED)
00502 return NOX::Abstract::Group::NotConverged;
00503
00504 return NOX::Abstract::Group::Failed;
00505 }
00506
00507 ::Thyra::ESolveMeasureNormType
00508 NOX::Thyra::Group::getThyraNormType(const string& name) const
00509 {
00510 if (name == "None")
00511 return ::Thyra::SOLVE_MEASURE_ONE;
00512 else if (name == "Norm Residual")
00513 return ::Thyra::SOLVE_MEASURE_NORM_RESIDUAL;
00514 else if (name == "Norm Solution")
00515 return ::Thyra::SOLVE_MEASURE_NORM_SOLUTION;
00516 else if (name == "Norm Initial Residual")
00517 return ::Thyra::SOLVE_MEASURE_NORM_INIT_RESIDUAL;
00518 else if (name == "Norm RHS")
00519 return ::Thyra::SOLVE_MEASURE_NORM_RHS;
00520 else {
00521 TEST_FOR_EXCEPTION(true, std::logic_error,
00522 "NOX Error - unknown solve measure " << name);
00523 return ::Thyra::SOLVE_MEASURE_ONE;
00524 }
00525 }