Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_EpetraModelEvaluator.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "Thyra_EpetraModelEvaluator.hpp"
00030 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00031 #include "Thyra_EpetraThyraWrappers.hpp"
00032 #include "Thyra_EpetraLinearOp.hpp"
00033 #include "Thyra_DetachedMultiVectorView.hpp"
00034 #include "Thyra_ModelEvaluatorDelegatorBase.hpp" // Gives verbose macros!
00035 #include "EpetraExt_ModelEvaluatorScalingTools.h"
00036 #include "Epetra_RowMatrix.h"
00037 #include "Teuchos_Time.hpp"
00038 #include "Teuchos_implicit_cast.hpp"
00039 #include "Teuchos_Assert.hpp"
00040 #include "Teuchos_StandardParameterEntryValidators.hpp"
00041 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00042 
00043 
00044 namespace {
00045 
00046 
00047 const std::string StateFunctionScaling_name = "State Function Scaling";
00048 Teuchos::RCP<
00049   Teuchos::StringToIntegralParameterEntryValidator<
00050     Thyra::EpetraModelEvaluator::EStateFunctionScaling
00051     >
00052   >
00053 stateFunctionScalingValidator;
00054 const std::string StateFunctionScaling_default = "None";
00055 
00056 
00057 // Extract out the Epetra_RowMatrix from the set W in an Epetra OutArgs object
00058 Teuchos::RCP<Epetra_RowMatrix>
00059 get_Epetra_RowMatrix(
00060   const EpetraExt::ModelEvaluator::OutArgs &epetraOutArgs
00061   )
00062 {
00063   using Teuchos::RCP;
00064   const RCP<Epetra_Operator>
00065     eW = epetraOutArgs.get_W();
00066   const RCP<Epetra_RowMatrix>
00067     ermW = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(eW,false);
00068   TEST_FOR_EXCEPTION(
00069     is_null(ermW), std::logic_error,
00070     "Thyra::EpetraModelEvaluator::evalModel(...): Error, if\n"
00071     "scaling is turned on, the underlying Epetra_Operator created\n"
00072     "an initialized by the underlying epetra model evaluator\n"
00073     "\"" << epetraOutArgs.modelEvalDescription() << "\"\n"
00074     "must support the Epetra_RowMatrix interface through a dynamic cast.\n"
00075     "The concrete type " << Teuchos::typeName(*eW) << " does not support\n"
00076     "Epetra_RowMatrix!"
00077     );
00078   return ermW;
00079 }
00080 
00081 
00082 Teuchos::RCP<Epetra_Operator>
00083 create_and_assert_W( 
00084   const EpetraExt::ModelEvaluator &epetraModel
00085   )
00086 {
00087   using Teuchos::RCP;
00088   RCP<Epetra_Operator>
00089     eW = epetraModel.create_W();
00090   TEST_FOR_EXCEPTION(
00091     is_null(eW), std::logic_error,
00092     "Error, the call to create_W() returned null on the "
00093     "EpetraExt::ModelEvaluator object "
00094     "\"" << epetraModel.description() << "\".  This may mean that "
00095     "the underlying model does not support more than one copy of "
00096     "W at one time!" );
00097   return eW;
00098 }
00099 
00100 
00101 } // namespace
00102 
00103 
00104 namespace Thyra {
00105 
00106 
00107 // Constructors/initializers/accessors.
00108 
00109 
00110 EpetraModelEvaluator::EpetraModelEvaluator()
00111   :nominalValuesAndBoundsAreUpdated_(false), stateFunctionScaling_(STATE_FUNC_SCALING_NONE),
00112    currentInArgsOutArgs_(false), finalPointWasSolved_(false)
00113 {}
00114 
00115 
00116 EpetraModelEvaluator::EpetraModelEvaluator(
00117   const RCP<const EpetraExt::ModelEvaluator> &epetraModel,
00118   const RCP<LinearOpWithSolveFactoryBase<double> > &W_factory
00119   )
00120   :nominalValuesAndBoundsAreUpdated_(false), stateFunctionScaling_(STATE_FUNC_SCALING_NONE),
00121    currentInArgsOutArgs_(false), finalPointWasSolved_(false)
00122 {
00123   initialize(epetraModel,W_factory);
00124 }
00125 
00126 
00127 void EpetraModelEvaluator::initialize(
00128   const RCP<const EpetraExt::ModelEvaluator> &epetraModel,
00129   const RCP<LinearOpWithSolveFactoryBase<double> > &W_factory
00130   )
00131 {
00132   using Teuchos::implicit_cast;
00133   typedef ModelEvaluatorBase MEB;
00134   //
00135   epetraModel_ = epetraModel;
00136   //
00137   W_factory_ = W_factory;
00138   //
00139   x_map_ = epetraModel_->get_x_map();
00140   f_map_ = epetraModel_->get_f_map();
00141   if (!is_null(x_map_)) {
00142     x_space_ = create_VectorSpace(x_map_);
00143     f_space_ = create_VectorSpace(f_map_);
00144   }
00145   //
00146   EpetraExt::ModelEvaluator::InArgs inArgs = epetraModel_->createInArgs();
00147   p_map_.resize(inArgs.Np()); p_space_.resize(inArgs.Np());
00148   p_map_is_local_.resize(inArgs.Np(),false);
00149   for( int l = 0; l < implicit_cast<int>(p_space_.size()); ++l ) {
00150     RCP<const Epetra_Map>
00151       p_map_l = ( p_map_[l] = epetraModel_->get_p_map(l) );
00152 #ifdef TEUCHOS_DEBUG
00153     TEST_FOR_EXCEPTION(
00154       is_null(p_map_l), std::logic_error,
00155       "Error, the the map p["<<l<<"] for the model \""
00156       <<epetraModel->description()<<"\" can not be null!");
00157 #endif
00158 
00159     p_map_is_local_[l] = !p_map_l->DistributedGlobal();
00160     p_space_[l] = create_VectorSpace(p_map_l);
00161   }
00162   //
00163   EpetraExt::ModelEvaluator::OutArgs outArgs = epetraModel_->createOutArgs();
00164   g_map_.resize(outArgs.Ng()); g_space_.resize(outArgs.Ng());
00165   g_map_is_local_.resize(outArgs.Ng(),false);
00166   for( int j = 0; j < implicit_cast<int>(g_space_.size()); ++j ) {
00167     RCP<const Epetra_Map>
00168       g_map_j = ( g_map_[j] = epetraModel_->get_g_map(j) );
00169     g_map_is_local_[j] = !g_map_j->DistributedGlobal();
00170     g_space_[j] = create_VectorSpace( g_map_j );
00171   }
00172   //
00173   epetraInArgsScaling_ = epetraModel_->createInArgs();
00174   epetraOutArgsScaling_ = epetraModel_->createOutArgs();
00175   nominalValuesAndBoundsAreUpdated_ = false;
00176   finalPointWasSolved_ = false;
00177   stateFunctionScalingVec_ = Teuchos::null; // Must set new scaling!
00178   //
00179   currentInArgsOutArgs_ = false;
00180 }
00181 
00182 
00183 RCP<const EpetraExt::ModelEvaluator>
00184 EpetraModelEvaluator::getEpetraModel() const
00185 {
00186   return epetraModel_;
00187 }
00188 
00189 
00190 void EpetraModelEvaluator::setNominalValues(
00191   const ModelEvaluatorBase::InArgs<double>& nominalValues
00192  )
00193 {
00194   nominalValues_.setArgs(nominalValues);
00195   // Note: These must be the scaled values so we don't need to scale!
00196 }
00197 
00198 
00199 void EpetraModelEvaluator::setStateVariableScalingVec(
00200   const RCP<const Epetra_Vector> &stateVariableScalingVec
00201   )
00202 {
00203   typedef ModelEvaluatorBase MEB;
00204 #ifdef TEUCHOS_DEBUG
00205   TEST_FOR_EXCEPT( !this->createInArgs().supports(MEB::IN_ARG_x) );
00206 #endif  
00207   stateVariableScalingVec_ = stateVariableScalingVec.assert_not_null();
00208   invStateVariableScalingVec_ = Teuchos::null;
00209   nominalValuesAndBoundsAreUpdated_ = false;
00210 }
00211 
00212 
00213 RCP<const Epetra_Vector>
00214 EpetraModelEvaluator::getStateVariableScalingVec() const
00215 {
00216   return stateVariableScalingVec_;
00217 }
00218 
00219 
00220 RCP<const Epetra_Vector>
00221 EpetraModelEvaluator::getStateVariableInvScalingVec() const
00222 {
00223   updateNominalValuesAndBounds();
00224   return invStateVariableScalingVec_;
00225 }
00226 
00227 
00228 void EpetraModelEvaluator::setStateFunctionScalingVec(
00229   const RCP<const Epetra_Vector> &stateFunctionScalingVec
00230   )
00231 {
00232   stateFunctionScalingVec_ = stateFunctionScalingVec;
00233 }
00234 
00235 
00236 RCP<const Epetra_Vector>
00237 EpetraModelEvaluator::getStateFunctionScalingVec() const
00238 {
00239   return stateFunctionScalingVec_;
00240 }
00241 
00242 
00243 void EpetraModelEvaluator::uninitialize(
00244   RCP<const EpetraExt::ModelEvaluator> *epetraModel,
00245   RCP<LinearOpWithSolveFactoryBase<double> > *W_factory
00246   )
00247 {
00248   if(epetraModel) *epetraModel = epetraModel_;
00249   if(W_factory) *W_factory = W_factory_;
00250   epetraModel_ = Teuchos::null;
00251   W_factory_ = Teuchos::null;
00252   stateFunctionScalingVec_ = Teuchos::null;
00253   stateVariableScalingVec_ = Teuchos::null;
00254   invStateVariableScalingVec_ = Teuchos::null;
00255   currentInArgsOutArgs_ = false;
00256 }
00257 
00258 
00259 const ModelEvaluatorBase::InArgs<double>&
00260 EpetraModelEvaluator::getFinalPoint() const
00261 {
00262   return finalPoint_;
00263 }
00264 
00265 
00266 bool EpetraModelEvaluator::finalPointWasSolved() const
00267 {
00268   return finalPointWasSolved_;
00269 }
00270 
00271 
00272 // Public functions overridden from Teuchos::Describable
00273 
00274 
00275 std::string EpetraModelEvaluator::description() const
00276 {
00277   std::ostringstream oss;
00278   oss << "Thyra::EpetraModelEvaluator{";
00279   oss << "epetraModel=";
00280   if(epetraModel_.get())
00281     oss << "\'"<<epetraModel_->description()<<"\'";
00282   else
00283     oss << "NULL";
00284   oss << ",W_factory=";
00285   if(W_factory_.get())
00286     oss << "\'"<<W_factory_->description()<<"\'";
00287   else
00288     oss << "NULL";
00289   oss << "}";
00290   return oss.str();
00291 }
00292 
00293 
00294 // Overridden from Teuchos::ParameterListAcceptor
00295 
00296 
00297 void EpetraModelEvaluator::setParameterList(
00298   RCP<Teuchos::ParameterList> const& paramList
00299   )
00300 {
00301   TEST_FOR_EXCEPT(is_null(paramList));
00302   paramList->validateParameters(*getValidParameters(),0); // Just validate my params
00303   paramList_ = paramList;
00304   const EStateFunctionScaling stateFunctionScaling_old = stateFunctionScaling_; 
00305   stateFunctionScaling_ = stateFunctionScalingValidator->getIntegralValue(
00306     *paramList_, StateFunctionScaling_name, StateFunctionScaling_default
00307     );
00308   if( stateFunctionScaling_ != stateFunctionScaling_old )
00309     stateFunctionScalingVec_ = Teuchos::null;
00310   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00311 #ifdef TEUCHOS_DEBUG
00312   paramList_->validateParameters(*getValidParameters(),0);
00313 #endif // TEUCHOS_DEBUG
00314 }
00315 
00316 
00317 RCP<Teuchos::ParameterList>
00318 EpetraModelEvaluator::getNonconstParameterList()
00319 {
00320   return paramList_;
00321 }
00322 
00323 
00324 RCP<Teuchos::ParameterList>
00325 EpetraModelEvaluator::unsetParameterList()
00326 {
00327   RCP<Teuchos::ParameterList> _paramList = paramList_;
00328   paramList_ = Teuchos::null;
00329   return _paramList;
00330 }
00331 
00332 
00333 RCP<const Teuchos::ParameterList>
00334 EpetraModelEvaluator::getParameterList() const
00335 {
00336   return paramList_;
00337 }
00338 
00339 
00340 RCP<const Teuchos::ParameterList>
00341 EpetraModelEvaluator::getValidParameters() const
00342 {
00343   using Teuchos::rcp;
00344   using Teuchos::StringToIntegralParameterEntryValidator;
00345   using Teuchos::tuple;
00346   using Teuchos::rcp_implicit_cast;
00347   typedef Teuchos::ParameterEntryValidator PEV;
00348   static RCP<const Teuchos::ParameterList> validPL;
00349   if(is_null(validPL)) {
00350     RCP<Teuchos::ParameterList>
00351       pl = Teuchos::rcp(new Teuchos::ParameterList());
00352     stateFunctionScalingValidator = rcp(
00353       new StringToIntegralParameterEntryValidator<EStateFunctionScaling>(
00354         tuple<std::string>(
00355           "None",
00356           "Row Sum"
00357           ),
00358         tuple<std::string>(
00359           "Do not scale the state function f(...) in this class.",
00360 
00361           "Scale the state function f(...) and all its derivatives\n"
00362           "using the row sum scaling from the initial Jacobian\n"
00363           "W=d(f)/d(x).  Note, this only works with Epetra_CrsMatrix\n"
00364           "currently."
00365           ),
00366         tuple<EStateFunctionScaling>(
00367           STATE_FUNC_SCALING_NONE,
00368           STATE_FUNC_SCALING_ROW_SUM
00369           ),
00370         StateFunctionScaling_name
00371         )
00372       );
00373     pl->set(StateFunctionScaling_name,StateFunctionScaling_default,
00374       "Determines if and how the state function f(...) and all of its\n"
00375       "derivatives are scaled.  The scaling is done explicitly so there should\n"
00376       "be no impact on the meaning of inner products or tolerances for\n"
00377       "linear solves.",
00378       rcp_implicit_cast<const PEV>(stateFunctionScalingValidator)
00379       );
00380     Teuchos::setupVerboseObjectSublist(&*pl);
00381     validPL = pl;
00382   }
00383   return validPL;
00384 }
00385 
00386 
00387 // Overridden from ModelEvaulator.
00388 
00389 
00390 int EpetraModelEvaluator::Np() const
00391 {
00392   return p_space_.size();
00393 }
00394 
00395 
00396 int EpetraModelEvaluator::Ng() const
00397 {
00398   return g_space_.size();
00399 }
00400 
00401 
00402 RCP<const VectorSpaceBase<double> >
00403 EpetraModelEvaluator::get_x_space() const
00404 {
00405   return x_space_;
00406 }
00407 
00408 
00409 RCP<const VectorSpaceBase<double> >
00410 EpetraModelEvaluator::get_f_space() const
00411 {
00412   return f_space_;
00413 }
00414 
00415 
00416 RCP<const VectorSpaceBase<double> >
00417 EpetraModelEvaluator::get_p_space(int l) const
00418 {
00419 #ifdef TEUCHOS_DEBUG
00420   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, this->Np() );
00421 #endif
00422   return p_space_[l];
00423 }
00424 
00425 
00426 RCP<const Teuchos::Array<std::string> >
00427 EpetraModelEvaluator::get_p_names(int l) const
00428 {
00429 #ifdef TEUCHOS_DEBUG
00430   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, this->Np() );
00431 #endif
00432   return epetraModel_->get_p_names(l);
00433 }
00434 
00435 
00436 RCP<const VectorSpaceBase<double> >
00437 EpetraModelEvaluator::get_g_space(int j) const
00438 {
00439   TEST_FOR_EXCEPT( ! ( 0 <= j && j < this->Ng() ) );
00440   return g_space_[j];
00441 }
00442 
00443 
00444 ModelEvaluatorBase::InArgs<double>
00445 EpetraModelEvaluator::getNominalValues() const
00446 {
00447   updateNominalValuesAndBounds();
00448   return nominalValues_;
00449 }
00450 
00451 
00452 ModelEvaluatorBase::InArgs<double>
00453 EpetraModelEvaluator::getLowerBounds() const
00454 {
00455   updateNominalValuesAndBounds();
00456   return lowerBounds_;
00457 }
00458 
00459 
00460 ModelEvaluatorBase::InArgs<double>
00461 EpetraModelEvaluator::getUpperBounds() const
00462 {
00463   updateNominalValuesAndBounds();
00464   return upperBounds_;
00465 }
00466 
00467 
00468 RCP<LinearOpWithSolveBase<double> >
00469 EpetraModelEvaluator::create_W() const
00470 {
00471   TEST_FOR_EXCEPTION(
00472     W_factory_.get()==NULL, std::logic_error
00473     ,"Thyra::EpetraModelEvaluator::create_W(): "
00474     "Error, the client did not set a LinearOpWithSolveFactoryBase"
00475     " object for W!"
00476     );
00477   W_factory_->setOStream(this->getOStream());
00478   W_factory_->setVerbLevel(this->getVerbLevel());
00479   return W_factory_->createOp();
00480 }
00481 
00482 
00483 RCP<LinearOpBase<double> >
00484 EpetraModelEvaluator::create_W_op() const
00485 {
00486   return this->create_epetra_W_op();
00487 }
00488 
00489 
00490 RCP<const LinearOpWithSolveFactoryBase<double> >
00491 EpetraModelEvaluator::get_W_factory() const
00492 {
00493   return W_factory_;
00494 }
00495 
00496 
00497 ModelEvaluatorBase::InArgs<double> EpetraModelEvaluator::createInArgs() const
00498 {
00499   if (!currentInArgsOutArgs_)
00500     updateInArgsOutArgs();
00501   return prototypeInArgs_;
00502 }
00503 
00504 
00505 void EpetraModelEvaluator::reportFinalPoint(
00506   const ModelEvaluatorBase::InArgs<double> &finalPoint,
00507   const bool wasSolved
00508   )
00509 {
00510   finalPoint_ = this->createInArgs();
00511   finalPoint_.setArgs(finalPoint);
00512   finalPointWasSolved_ = wasSolved;
00513 }
00514 
00515 
00516 // Private functions overridden from ModelEvaulatorDefaultBase
00517 
00518 
00519 RCP<LinearOpBase<double> >
00520 EpetraModelEvaluator::create_DfDp_op_impl(int l) const
00521 {
00522   TEST_FOR_EXCEPT(true);
00523   return Teuchos::null;
00524 }
00525 
00526 
00527 RCP<LinearOpBase<double> >
00528 EpetraModelEvaluator::create_DgDx_dot_op_impl(int j) const
00529 {
00530   TEST_FOR_EXCEPT(true);
00531   return Teuchos::null;
00532 }
00533 
00534 
00535 RCP<LinearOpBase<double> >
00536 EpetraModelEvaluator::create_DgDx_op_impl(int j) const
00537 {
00538   TEST_FOR_EXCEPT(true);
00539   return Teuchos::null;
00540 }
00541 
00542 
00543 RCP<LinearOpBase<double> >
00544 EpetraModelEvaluator::create_DgDp_op_impl( int j, int l ) const
00545 {
00546   TEST_FOR_EXCEPT(true);
00547   return Teuchos::null;
00548 }
00549 
00550 
00551 ModelEvaluatorBase::OutArgs<double>
00552 EpetraModelEvaluator::createOutArgsImpl() const
00553 {
00554   if (!currentInArgsOutArgs_)
00555     updateInArgsOutArgs();
00556   return prototypeOutArgs_;
00557 }
00558 
00559 
00560 void EpetraModelEvaluator::evalModelImpl(
00561   const ModelEvaluatorBase::InArgs<double>& inArgs_in,
00562   const ModelEvaluatorBase::OutArgs<double>& outArgs
00563   ) const
00564 {
00565 
00566   using Teuchos::rcp;
00567   using Teuchos::rcp_const_cast;
00568   using Teuchos::rcp_dynamic_cast;
00569   using Teuchos::OSTab;
00570   using Teuchos::includesVerbLevel;
00571   typedef EpetraExt::ModelEvaluator EME;
00572 
00573   //
00574   // A) Initial setup
00575   //
00576 
00577   // Make sure that we are fully initialized!
00578   this->updateNominalValuesAndBounds();
00579 
00580   // Make sure we grab the initial guess first!
00581   InArgs<double> inArgs = this->getNominalValues();
00582   // Now, copy the parameters from the input inArgs_in object to the inArgs
00583   // object.  Any input objects that are set in inArgs_in will overwrite those
00584   // in inArgs that will already contain the nominal values.  This will insure
00585   // that all input parameters are set and those that are not set by the
00586   // client will be at their nominal values (as defined by the underlying
00587   // EpetraExt::ModelEvaluator object).  The full set of Thyra input arguments
00588   // must be set before these can be translated into Epetra input arguments.
00589   inArgs.setArgs(inArgs_in);
00590 
00591   // Print the header and the values of the inArgs and outArgs objects!
00592   typedef double Scalar; // Needed for below macro!
00593   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00594     "Thyra::EpetraModelEvaluator",inArgs,outArgs,Teuchos::null
00595     );
00596 
00597   // State function Scaling
00598   const bool firstTimeStateFuncScaling
00599     = (
00600       stateFunctionScaling_ != STATE_FUNC_SCALING_NONE
00601       && is_null(stateFunctionScalingVec_)
00602       );
00603   
00604   typedef Teuchos::VerboseObjectTempState<LinearOpWithSolveFactoryBase<double> > VOTSLOWSF;
00605   VOTSLOWSF W_factory_outputTempState(W_factory_,out,verbLevel);
00606 
00607   Teuchos::Time timer("");
00608 
00609   //
00610   // B) Prepressess the InArgs and OutArgs in preparation to call
00611   // the underlying EpetraExt::ModelEvaluator
00612   //
00613 
00614   //
00615   // B.1) Translate InArgs from Thyra to Epetra objects
00616   //
00617   
00618   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00619     *out << "\nSetting-up/creating input arguments ...\n";
00620   timer.start(true);
00621 
00622   // Unwrap input Thyra objects to get Epetra objects
00623   EME::InArgs epetraScaledInArgs = epetraModel_->createInArgs();
00624   convertInArgsFromThyraToEpetra( inArgs, &epetraScaledInArgs );
00625 
00626   // Unscale the input Epetra objects which will be passed to the underlying
00627   // EpetraExt::ModelEvaluator object.
00628   EME::InArgs epetraInArgs = epetraModel_->createInArgs();
00629   EpetraExt::unscaleModelVars(
00630     epetraScaledInArgs, epetraInArgsScaling_, &epetraInArgs,
00631     out.get(), verbLevel
00632     );
00633 
00634   timer.stop();
00635   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00636     OSTab(out).o() << "\nTime to setup InArgs = "<<timer.totalElapsedTime()<<" sec\n";
00637 
00638   //
00639   // B.2) Convert from Thyra to Epetra OutArgs
00640   //
00641   
00642   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00643     *out << "\nSetting-up/creating output arguments ...\n";
00644   timer.start(true);
00645   
00646   // The unscaled Epetra OutArgs that will be passed to the
00647   // underlying EpetraExt::ModelEvaluator object
00648   EME::OutArgs epetraUnscaledOutArgs = epetraModel_->createOutArgs();
00649 
00650   // Various objects that are needed later (see documentation in
00651   // the function convertOutArgsFromThyraToEpetra(...)
00652   RCP<LinearOpWithSolveBase<double> > W;
00653   RCP<LinearOpBase<double> > W_op;
00654   RCP<const LinearOpBase<double> > fwdW;
00655   RCP<EpetraLinearOp> efwdW;
00656   RCP<Epetra_Operator> eW;
00657   
00658   // Convert from Thyra to Epetra OutArgs and grap some of the intermediate
00659   // objects accessed along the way that are needed later.
00660   convertOutArgsFromThyraToEpetra(
00661     outArgs,
00662     &epetraUnscaledOutArgs,
00663     &W, &W_op, &fwdW, &efwdW, &eW
00664     );
00665 
00666   //
00667   // B.3) Setup OutArgs to computing scaling if needed
00668   //
00669 
00670   if (firstTimeStateFuncScaling) {
00671     preEvalScalingSetup(&epetraInArgs,&epetraUnscaledOutArgs,out,verbLevel);
00672   }
00673 
00674   timer.stop();
00675   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00676     OSTab(out).o()
00677       << "\nTime to setup OutArgs = "
00678       << timer.totalElapsedTime() <<" sec\n";
00679 
00680   //
00681   // C) Evaluate the underlying EpetraExt model to compute the Epetra outputs
00682   //
00683   
00684   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00685     *out << "\nEvaluating the Epetra output functions ...\n";
00686   timer.start(true);
00687 
00688   epetraModel_->evalModel(epetraInArgs,epetraUnscaledOutArgs);
00689 
00690   timer.stop();
00691   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00692     OSTab(out).o()
00693       << "\nTime to evaluate Epetra output functions = "
00694       << timer.totalElapsedTime() <<" sec\n";
00695 
00696   //
00697   // D) Postprocess the output objects
00698   //
00699 
00700   //
00701   // D.1) Compute the scaling factors if needed
00702   //
00703   
00704   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00705     *out << "\nCompute scale factors if needed ...\n";
00706   timer.start(true);
00707 
00708   if (firstTimeStateFuncScaling) {
00709     postEvalScalingSetup(epetraUnscaledOutArgs,out,verbLevel);
00710   }
00711 
00712   timer.stop();
00713   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00714     OSTab(out).o()
00715       << "\nTime to compute scale factors = "
00716       << timer.totalElapsedTime() <<" sec\n";
00717 
00718   //
00719   // D.2) Scale the output Epetra objects
00720   //
00721 
00722   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00723     *out << "\nScale the output objects ...\n";
00724   timer.start(true);
00725 
00726   EME::OutArgs epetraOutArgs = epetraModel_->createOutArgs();
00727   bool allFuncsWhereScaled = false;
00728   EpetraExt::scaleModelFuncs(
00729     epetraUnscaledOutArgs, epetraInArgsScaling_, epetraOutArgsScaling_,
00730     &epetraOutArgs, &allFuncsWhereScaled,
00731     out.get(), verbLevel
00732     );
00733   TEST_FOR_EXCEPTION(
00734     !allFuncsWhereScaled, std::logic_error,
00735     "Error, we can not currently handle epetra output objects that could not be"
00736     " scaled.  Special code will have to be added to handle this (i.e. using"
00737     " implicit diagonal and multiplied linear operators to implicitly do"
00738     " the scaling."
00739     );
00740 
00741   timer.stop();
00742   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00743     OSTab(out).o()
00744       << "\nTime to scale the output objects = "
00745       << timer.totalElapsedTime() << " sec\n";
00746 
00747   //
00748   // D.3) Convert any Epetra objects to Thyra OutArgs objects that still need to
00749   // be converted
00750   //
00751 
00752   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00753     *out << "\nFinish processing and wrapping the output objects ...\n";
00754   timer.start(true);
00755 
00756   finishConvertingOutArgsFromEpetraToThyra(
00757     epetraOutArgs, W, W_op, fwdW, efwdW, eW,
00758     outArgs
00759     );
00760 
00761   timer.stop();
00762   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00763     OSTab(out).o()
00764       << "\nTime to finish processing and wrapping the output objects = "
00765       << timer.totalElapsedTime() <<" sec\n";
00766 
00767   //
00768   // E) Print footer to end the function
00769   //
00770 
00771   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00772   
00773 }
00774 
00775 
00776 // private
00777 
00778 
00779 void EpetraModelEvaluator::convertInArgsFromEpetraToThyra(
00780   const EpetraExt::ModelEvaluator::InArgs &epetraInArgs,
00781   ModelEvaluatorBase::InArgs<double> *inArgs
00782   ) const
00783 {
00784   
00785   using Teuchos::implicit_cast;
00786   typedef ModelEvaluatorBase MEB;
00787 
00788   TEST_FOR_EXCEPT(!inArgs);
00789 
00790   if(inArgs->supports(MEB::IN_ARG_x)) {
00791     inArgs->set_x( create_Vector( epetraInArgs.get_x(), x_space_ ) );
00792   }
00793   
00794   if(inArgs->supports(MEB::IN_ARG_x_dot)) {
00795     inArgs->set_x_dot( create_Vector( epetraInArgs.get_x_dot(), x_space_ ) );
00796   }
00797 
00798   const int l_Np = inArgs->Np();
00799   for( int l = 0; l < l_Np; ++l ) {
00800     inArgs->set_p( l, create_Vector( epetraInArgs.get_p(l), p_space_[l] ) );
00801   }
00802   
00803   if(inArgs->supports(MEB::IN_ARG_t)) {
00804     inArgs->set_t(epetraInArgs.get_t());
00805   }
00806   
00807 }
00808 
00809 
00810 void EpetraModelEvaluator::convertInArgsFromThyraToEpetra(
00811   const ModelEvaluatorBase::InArgs<double> &inArgs,
00812   EpetraExt::ModelEvaluator::InArgs *epetraInArgs
00813   ) const
00814 {
00815 
00816   using Teuchos::rcp;
00817   using Teuchos::rcp_const_cast;
00818 #ifdef HAVE_THYRA_ME_POLYNOMIAL
00819   using Teuchos::Polynomial;
00820 #endif // HAVE_THYRA_ME_POLYNOMIAL
00821 
00822 
00823   TEST_FOR_EXCEPT(0==epetraInArgs);
00824 
00825   RCP<const VectorBase<double> > x_dot;
00826   if( inArgs.supports(IN_ARG_x_dot) && (x_dot = inArgs.get_x_dot()).get() ) {
00827     RCP<const Epetra_Vector> e_x_dot = get_Epetra_Vector(*x_map_,x_dot);
00828     epetraInArgs->set_x_dot(e_x_dot);
00829   }
00830 
00831   RCP<const VectorBase<double> > x;
00832   if( inArgs.supports(IN_ARG_x) && (x = inArgs.get_x()).get() ) {
00833     RCP<const Epetra_Vector> e_x = get_Epetra_Vector(*x_map_,x);
00834     epetraInArgs->set_x(e_x);
00835   }
00836 
00837   RCP<const VectorBase<double> > p_l;
00838   for(int l = 0;  l < inArgs.Np(); ++l ) {
00839     p_l = inArgs.get_p(l);
00840     if(p_l.get()) epetraInArgs->set_p(l,get_Epetra_Vector(*p_map_[l],p_l));
00841   }
00842 
00843 #ifdef HAVE_THYRA_ME_POLYNOMIAL
00844 
00845   RCP<const Polynomial< VectorBase<double> > > x_dot_poly;
00846   RCP<Epetra_Vector> epetra_ptr;
00847   if(
00848     inArgs.supports(IN_ARG_x_dot_poly)
00849     && (x_dot_poly = inArgs.get_x_dot_poly()).get()
00850     )
00851   {
00852     RCP<Polynomial<Epetra_Vector> > epetra_x_dot_poly = 
00853       rcp(new Polynomial<Epetra_Vector>(x_dot_poly->degree()));
00854     for (unsigned int i=0; i<=x_dot_poly->degree(); i++) {
00855       epetra_ptr = rcp_const_cast<Epetra_Vector>(
00856         get_Epetra_Vector(*x_map_, x_dot_poly->getCoefficient(i)) );
00857       epetra_x_dot_poly->setCoefficientPtr(i,epetra_ptr);
00858     }
00859     epetraInArgs->set_x_dot_poly(epetra_x_dot_poly);
00860   }
00861   
00862   RCP<const Polynomial< VectorBase<double> > > x_poly;
00863   if(
00864     inArgs.supports(IN_ARG_x_poly)
00865     && (x_poly = inArgs.get_x_poly()).get()
00866     )
00867   {
00868     RCP<Polynomial<Epetra_Vector> > epetra_x_poly = 
00869       rcp(new Polynomial<Epetra_Vector>(x_poly->degree()));
00870     for (unsigned int i=0; i<=x_poly->degree(); i++) {
00871       epetra_ptr = rcp_const_cast<Epetra_Vector>(
00872         get_Epetra_Vector(*x_map_, x_poly->getCoefficient(i)) );
00873       epetra_x_poly->setCoefficientPtr(i,epetra_ptr);
00874     }
00875     epetraInArgs->set_x_poly(epetra_x_poly);
00876   }
00877 
00878 #endif // HAVE_THYRA_ME_POLYNOMIAL
00879 
00880   if( inArgs.supports(IN_ARG_t) )
00881     epetraInArgs->set_t(inArgs.get_t());
00882   
00883   if( inArgs.supports(IN_ARG_alpha) )
00884     epetraInArgs->set_alpha(inArgs.get_alpha());
00885   
00886   if( inArgs.supports(IN_ARG_beta) )
00887     epetraInArgs->set_beta(inArgs.get_beta());
00888 
00889 }
00890 
00891 
00892 void EpetraModelEvaluator::convertOutArgsFromThyraToEpetra(
00893   const ModelEvaluatorBase::OutArgs<double> &outArgs,
00894   EpetraExt::ModelEvaluator::OutArgs *epetraUnscaledOutArgs_inout,
00895   RCP<LinearOpWithSolveBase<double> > *W_out,
00896   RCP<LinearOpBase<double> > *W_op_out,
00897   RCP<const LinearOpBase<double> > *fwdW_out,
00898   RCP<EpetraLinearOp> *efwdW_out,
00899   RCP<Epetra_Operator> *eW_out
00900   ) const
00901 {
00902 
00903   using Teuchos::rcp;
00904   using Teuchos::rcp_const_cast;
00905   using Teuchos::rcp_dynamic_cast;
00906   using Teuchos::OSTab;
00907   using Teuchos::implicit_cast;
00908   using Thyra::get_Epetra_Vector;
00909   typedef EpetraExt::ModelEvaluator EME;
00910 
00911   // Assert input
00912 #ifdef TEUCHOS_DEBUG
00913   TEUCHOS_ASSERT(epetraUnscaledOutArgs_inout);
00914   TEUCHOS_ASSERT(W_out);
00915   TEUCHOS_ASSERT(W_op_out);
00916   TEUCHOS_ASSERT(fwdW_out);
00917   TEUCHOS_ASSERT(efwdW_out);
00918   TEUCHOS_ASSERT(eW_out);
00919 #endif
00920 
00921   // Create easy to use references
00922   EpetraExt::ModelEvaluator::OutArgs &epetraUnscaledOutArgs = *epetraUnscaledOutArgs_inout;
00923   RCP<LinearOpWithSolveBase<double> > &W = *W_out;
00924   RCP<LinearOpBase<double> > &W_op = *W_op_out;
00925   RCP<const LinearOpBase<double> > &fwdW = *fwdW_out;
00926   RCP<EpetraLinearOp> &efwdW = *efwdW_out;
00927   RCP<Epetra_Operator> &eW = *eW_out;
00928 
00929   // f
00930   { 
00931     RCP<VectorBase<double> > f;
00932     if( outArgs.supports(OUT_ARG_f) && (f = outArgs.get_f()).get() )
00933       epetraUnscaledOutArgs.set_f(get_Epetra_Vector(*f_map_,f));
00934   }
00935     
00936   // g
00937   {
00938     RCP<VectorBase<double> > g_j;
00939     for(int j = 0;  j < outArgs.Ng(); ++j ) {
00940       g_j = outArgs.get_g(j);
00941       if(g_j.get()) epetraUnscaledOutArgs.set_g(j,get_Epetra_Vector(*g_map_[j],g_j));
00942     }
00943   }
00944   
00945   // W and W_op
00946   {
00947 
00948     if( outArgs.supports(OUT_ARG_W) && (W = outArgs.get_W()).get() ) {
00949       Thyra::uninitializeOp<double>(*W_factory_, W.ptr(), Teuchos::outArg(fwdW));
00950       if(fwdW.get()) {
00951         efwdW = rcp_const_cast<EpetraLinearOp>(
00952           rcp_dynamic_cast<const EpetraLinearOp>(fwdW,true));
00953       }
00954       else {
00955         efwdW = this->create_epetra_W_op();
00956         fwdW = efwdW;
00957       }
00958     }
00959     if( outArgs.supports(OUT_ARG_W_op) && (W_op = outArgs.get_W_op()).get() ) {
00960       if( W_op.get() && !efwdW.get() )
00961         efwdW = rcp_const_cast<EpetraLinearOp>(
00962           rcp_dynamic_cast<const EpetraLinearOp>(W_op,true));
00963     }
00964     if(efwdW.get()) {
00965       // By the time we get here, if we have an object in efwdW, then it
00966       // should already be embeadded with an underlying Epetra_Operator object
00967       // that was allocated by the EpetraExt::ModelEvaluator object.
00968       // Therefore, we should just have to grab this object and be on our way.
00969       eW = efwdW->epetra_op();
00970       epetraUnscaledOutArgs.set_W(eW);
00971     }
00972     // NOTE: Above, if both W and W_op are set and have been through at least
00973     // one prior evaluation (and therefore have Epetra_Operator objects embedded
00974     // in them), then we will use the Epetra_Operator embedded in W to pass to
00975     // the EpetraExt::ModelEvaluator object and ignore the Epetra_Operator
00976     // object in W_op.  In the standard use case, these will be the same
00977     // Epetra_Operator objects.  However, it is possible that the client could
00978     // use this interface in such a way that these would have different
00979     // Epetra_Operator objects embedded in them.  In this (very unlikely) case,
00980     // the Epetra_Operator embedded in W_op will be discarded!  This might be
00981     // surprising to a client but it is very unlikely that this will ever be a
00982     // problem, but the issue is duly noted here!  Only dangerous programming
00983     // use of this interface would cause any problem.
00984     
00985     // Note: The following derivative objects update in place!
00986 
00987   }
00988 
00989   // DfDp
00990   {
00991     Derivative<double> DfDp_l;
00992     for(int l = 0;  l < outArgs.Np(); ++l ) {
00993       if( !outArgs.supports(OUT_ARG_DfDp,l).none()
00994         && !(DfDp_l = outArgs.get_DfDp(l)).isEmpty() )
00995       {
00996         epetraUnscaledOutArgs.set_DfDp(l,convert(DfDp_l,f_map_,p_map_[l]));
00997       }
00998     }
00999   }
01000 
01001   // DgDx_dot
01002   {
01003     Derivative<double> DgDx_dot_j;
01004     for(int j = 0;  j < outArgs.Ng(); ++j ) {
01005       if( !outArgs.supports(OUT_ARG_DgDx_dot,j).none()
01006         && !(DgDx_dot_j = outArgs.get_DgDx_dot(j)).isEmpty() )
01007       {
01008         epetraUnscaledOutArgs.set_DgDx_dot(j,convert(DgDx_dot_j,g_map_[j],x_map_));
01009       }
01010     }
01011   }
01012 
01013   // DgDx
01014   {
01015     Derivative<double> DgDx_j;
01016     for(int j = 0;  j < outArgs.Ng(); ++j ) {
01017       if( !outArgs.supports(OUT_ARG_DgDx,j).none()
01018         && !(DgDx_j = outArgs.get_DgDx(j)).isEmpty() )
01019       {
01020         epetraUnscaledOutArgs.set_DgDx(j,convert(DgDx_j,g_map_[j],x_map_));
01021       }
01022     }
01023   }
01024 
01025   // DgDp
01026   {
01027     DerivativeSupport DgDp_j_l_support;
01028     Derivative<double> DgDp_j_l;
01029     for (int j = 0;  j < outArgs.Ng(); ++j ) {
01030       for (int l = 0;  l < outArgs.Np(); ++l ) {
01031         if (!(DgDp_j_l_support = outArgs.supports(OUT_ARG_DgDp,j,l)).none()
01032           && !(DgDp_j_l = outArgs.get_DgDp(j,l)).isEmpty() )
01033         {
01034           epetraUnscaledOutArgs.set_DgDp(j,l,convert(DgDp_j_l,g_map_[j],p_map_[l]));
01035         }
01036       }
01037     }
01038   }
01039 
01040 #ifdef HAVE_THYRA_ME_POLYNOMIAL
01041 
01042   // f_poly
01043   RCP<const Teuchos::Polynomial< VectorBase<double> > > f_poly;
01044   if( outArgs.supports(OUT_ARG_f_poly) && (f_poly = outArgs.get_f_poly()).get() )
01045   {
01046     RCP<Teuchos::Polynomial<Epetra_Vector> > epetra_f_poly = 
01047       Teuchos::rcp(new Teuchos::Polynomial<Epetra_Vector>(f_poly->degree()));
01048     for (unsigned int i=0; i<=f_poly->degree(); i++) {
01049       RCP<Epetra_Vector> epetra_ptr
01050         = Teuchos::rcp_const_cast<Epetra_Vector>(get_Epetra_Vector(*f_map_,
01051             f_poly->getCoefficient(i)));
01052       epetra_f_poly->setCoefficientPtr(i,epetra_ptr);
01053     }
01054     epetraUnscaledOutArgs.set_f_poly(epetra_f_poly);
01055   }
01056 
01057 #endif // HAVE_THYRA_ME_POLYNOMIAL
01058 
01059 }
01060 
01061 
01062 void EpetraModelEvaluator::preEvalScalingSetup(
01063   EpetraExt::ModelEvaluator::InArgs *epetraInArgs_inout,
01064   EpetraExt::ModelEvaluator::OutArgs *epetraUnscaledOutArgs_inout,
01065   const RCP<Teuchos::FancyOStream> &out,
01066   const Teuchos::EVerbosityLevel verbLevel
01067   ) const
01068 {
01069   
01070   typedef EpetraExt::ModelEvaluator EME;
01071   
01072 #ifdef TEUCHOS_DEBUG
01073   TEUCHOS_ASSERT(epetraInArgs_inout);
01074   TEUCHOS_ASSERT(epetraUnscaledOutArgs_inout);
01075 #endif
01076 
01077   EpetraExt::ModelEvaluator::InArgs
01078     &epetraInArgs = *epetraInArgs_inout;
01079   EpetraExt::ModelEvaluator::OutArgs
01080     &epetraUnscaledOutArgs = *epetraUnscaledOutArgs_inout;
01081 
01082   if (
01083     ( stateFunctionScaling_ == STATE_FUNC_SCALING_ROW_SUM )
01084     &&
01085     (
01086       epetraUnscaledOutArgs.supports(EME::OUT_ARG_f) 
01087       &&
01088       epetraUnscaledOutArgs.funcOrDerivesAreSet(EME::OUT_ARG_f)
01089       )
01090     &&
01091     (
01092       epetraUnscaledOutArgs.supports(EME::OUT_ARG_W)
01093       &&
01094       is_null(epetraUnscaledOutArgs.get_W())
01095       )
01096     )
01097   {
01098     // This is the first pass through with scaling turned on and the client
01099     // turned on automatic scaling but did not ask for W.  We must compute W
01100     // in order to compute the scale factors so we must allocate a temporary W
01101     // just to compute the scale factors and then throw it away.  If the
01102     // client wants to evaluate W at the same point, then it should have
01103     // passed W in but that is not our problem here.  The ModelEvaluator
01104     // relies on the client to set up the calls to allow for efficient
01105     // evaluation.
01106 
01107     if(out.get() && verbLevel >= Teuchos::VERB_LOW)
01108       *out
01109         << "\nCreating a temporary Epetra W to compute scale factors"
01110         << " for f(...) ...\n";
01111     epetraUnscaledOutArgs.set_W(create_and_assert_W(*epetraModel_));
01112     if( epetraInArgs.supports(EME::IN_ARG_beta) )
01113       epetraInArgs.set_beta(1.0);
01114     if( epetraInArgs.supports(EME::IN_ARG_alpha) )
01115       epetraInArgs.set_alpha(0.0);
01116   }
01117   
01118 }
01119 
01120 
01121 void EpetraModelEvaluator::postEvalScalingSetup(
01122   const EpetraExt::ModelEvaluator::OutArgs &epetraUnscaledOutArgs,
01123   const RCP<Teuchos::FancyOStream> &out,
01124   const Teuchos::EVerbosityLevel verbLevel
01125   ) const
01126 {
01127 
01128   using Teuchos::OSTab;
01129   using Teuchos::rcp;
01130   using Teuchos::rcp_const_cast;
01131   using Teuchos::includesVerbLevel;
01132 
01133   // Compute the scale factors for the state function f(...)
01134   switch(stateFunctionScaling_) {
01135 
01136     case STATE_FUNC_SCALING_ROW_SUM: {
01137 
01138       // Compute the inverse row-sum scaling from W
01139 
01140       const RCP<Epetra_RowMatrix>
01141         ermW = get_Epetra_RowMatrix(epetraUnscaledOutArgs);
01142       // Note: Above, we get the Epetra W object directly from the Epetra
01143       // OutArgs object since this might be a temporary matrix just to
01144       // compute scaling factors.  In this case, the stack funtion variable
01145       // eW might be empty!
01146 
01147       RCP<Epetra_Vector>
01148         invRowSums = rcp(new Epetra_Vector(ermW->OperatorRangeMap()));
01149       // Above: From the documentation is seems that the RangeMap should be
01150       // okay but who knows for sure!
01151 
01152       ermW->InvRowSums(*invRowSums);
01153 
01154       if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW)) {
01155         *out
01156           << "\nComputed inverse row sum scaling from W that"
01157           " will be used to scale f(...) and its derivatives:\n";
01158         double minVal = 0, maxVal = 0, avgVal = 0;
01159         invRowSums->MinValue(&minVal);
01160         invRowSums->MaxValue(&maxVal);
01161         invRowSums->MeanValue(&avgVal);
01162         OSTab tab(out);
01163         *out
01164           << "min(invRowSums) = " << minVal << "\n"
01165           << "max(invRowSums) = " << maxVal << "\n"
01166           << "avg(invRowSums) = " << avgVal << "\n";
01167       }
01168 
01169       stateFunctionScalingVec_ = invRowSums;
01170 
01171       break;
01172 
01173     }
01174 
01175     default:
01176       TEST_FOR_EXCEPT("Should never get here!");
01177 
01178   }
01179 
01180   epetraOutArgsScaling_ = epetraModel_->createOutArgs();
01181 
01182   epetraOutArgsScaling_.set_f(
01183     rcp_const_cast<Epetra_Vector>(stateFunctionScalingVec_) );
01184 
01185 }
01186 
01187 
01188 void EpetraModelEvaluator::finishConvertingOutArgsFromEpetraToThyra(
01189   const EpetraExt::ModelEvaluator::OutArgs &epetraOutArgs,
01190   RCP<LinearOpWithSolveBase<double> > &W,
01191   RCP<LinearOpBase<double> > &W_op,
01192   RCP<const LinearOpBase<double> > &fwdW,
01193   RCP<EpetraLinearOp> &efwdW,
01194   RCP<Epetra_Operator> &eW,
01195   const ModelEvaluatorBase::OutArgs<double> &outArgs
01196   ) const
01197 {
01198 
01199   using Teuchos::rcp_dynamic_cast;
01200   typedef EpetraExt::ModelEvaluator EME;
01201 
01202   if(efwdW.get()) {
01203     efwdW->setFullyInitialized(true); 
01204     // NOTE: Above will directly update W_op also if W.get()==NULL!
01205   }
01206   
01207   if( W.get() ) {
01208     Thyra::initializeOp<double>(*W_factory_, fwdW, W.ptr());
01209     W->setOStream(this->getOStream());
01210   }
01211 
01212   if( W_op.get() ) {
01213     if( W_op.shares_resource(efwdW) ) {
01214       // W_op was already updated above since *efwdW is the same object as *W_op
01215     }
01216     else {
01217       rcp_dynamic_cast<EpetraLinearOp>(W_op,true)->setFullyInitialized(true);
01218     }
01219   }
01220 
01221 }
01222 
01223 
01224 void EpetraModelEvaluator::updateNominalValuesAndBounds() const
01225 {
01226 
01227   using Teuchos::rcp;
01228   using Teuchos::implicit_cast;
01229   typedef ModelEvaluatorBase MEB;
01230   typedef EpetraExt::ModelEvaluator EME;
01231 
01232   if( !nominalValuesAndBoundsAreUpdated_ ) {
01233 
01234     // Gather the nominal values and bounds into Epetra InArgs objects
01235 
01236     EME::InArgs epetraOrigNominalValues;
01237     EpetraExt::gatherModelNominalValues(
01238       *epetraModel_, &epetraOrigNominalValues );
01239 
01240     EME::InArgs epetraOrigLowerBounds;
01241     EME::InArgs epetraOrigUpperBounds;
01242     EpetraExt::gatherModelBounds(
01243       *epetraModel_, &epetraOrigLowerBounds, &epetraOrigUpperBounds );
01244 
01245     // Set up Epetra InArgs scaling object
01246 
01247     epetraInArgsScaling_ = epetraModel_->createInArgs();
01248 
01249     if( !is_null(stateVariableScalingVec_) ) {
01250       invStateVariableScalingVec_
01251         = EpetraExt::createInverseModelScalingVector(stateVariableScalingVec_);
01252       if( epetraOrigNominalValues.supports(EME::IN_ARG_x_dot) ) {
01253         epetraInArgsScaling_.set_x_dot(invStateVariableScalingVec_);
01254       }
01255       if( epetraOrigNominalValues.supports(EME::IN_ARG_x) ) {
01256         epetraInArgsScaling_.set_x(invStateVariableScalingVec_);
01257       }
01258     }
01259     
01260     // Scale the original variables and bounds
01261 
01262     EME::InArgs epetraScaledNominalValues = epetraModel_->createInArgs();
01263     EpetraExt::scaleModelVars(
01264       epetraOrigNominalValues, epetraInArgsScaling_, &epetraScaledNominalValues
01265       );
01266 
01267     EME::InArgs epetraScaledLowerBounds = epetraModel_->createInArgs();
01268     EME::InArgs epetraScaledUpperBounds = epetraModel_->createInArgs();
01269     EpetraExt::scaleModelBounds(
01270       epetraOrigLowerBounds, epetraOrigUpperBounds, epetraModel_->getInfBound(),
01271       epetraInArgsScaling_,
01272       &epetraScaledLowerBounds, &epetraScaledUpperBounds
01273       );
01274 
01275     // Wrap the scaled epetra InArgs objects as Thyra InArgs objects!
01276 
01277     nominalValues_ = this->createInArgs();
01278     lowerBounds_ = this->createInArgs();
01279     upperBounds_ = this->createInArgs();
01280     convertInArgsFromEpetraToThyra(epetraScaledNominalValues, &nominalValues_);
01281     convertInArgsFromEpetraToThyra(epetraScaledLowerBounds, &lowerBounds_);
01282     convertInArgsFromEpetraToThyra(epetraScaledUpperBounds, &upperBounds_);
01283 
01284     nominalValuesAndBoundsAreUpdated_ = true;
01285 
01286   }
01287   else {
01288 
01289     // The nominal values and bounds should already be updated an should have
01290     // the currect scaling!
01291 
01292   }
01293 
01294 }
01295 
01296 
01297 void EpetraModelEvaluator::updateInArgsOutArgs() const
01298 {
01299 
01300   typedef EpetraExt::ModelEvaluator EME;
01301 
01302   const EpetraExt::ModelEvaluator &epetraModel = *epetraModel_;
01303   EME::InArgs  epetraInArgs  = epetraModel.createInArgs();
01304   EME::OutArgs epetraOutArgs = epetraModel.createOutArgs();
01305   const int l_Np = epetraOutArgs.Np();
01306   const int l_Ng = epetraOutArgs.Ng();
01307 
01308   //
01309   // InArgs
01310   //
01311 
01312   InArgsSetup<double> inArgs;
01313   inArgs.setModelEvalDescription(this->description());
01314   inArgs.set_Np(epetraInArgs.Np());
01315   inArgs.setSupports(IN_ARG_x_dot, epetraInArgs.supports(EME::IN_ARG_x_dot));
01316   inArgs.setSupports(IN_ARG_x, epetraInArgs.supports(EME::IN_ARG_x));
01317 #ifdef HAVE_THYRA_ME_POLYNOMIAL
01318   inArgs.setSupports(IN_ARG_x_dot_poly,
01319     epetraInArgs.supports(EME::IN_ARG_x_dot_poly));
01320   inArgs.setSupports(IN_ARG_x_poly, epetraInArgs.supports(EME::IN_ARG_x_poly));
01321 #endif // HAVE_THYRA_ME_POLYNOMIAL
01322   inArgs.setSupports(IN_ARG_t, epetraInArgs.supports(EME::IN_ARG_t));
01323   inArgs.setSupports(IN_ARG_alpha, epetraInArgs.supports(EME::IN_ARG_alpha));
01324   inArgs.setSupports(IN_ARG_beta, epetraInArgs.supports(EME::IN_ARG_beta));
01325   prototypeInArgs_ = inArgs;
01326 
01327   //
01328   // OutArgs
01329   //
01330 
01331   OutArgsSetup<double> outArgs;
01332   outArgs.setModelEvalDescription(this->description());
01333   outArgs.set_Np_Ng(l_Np, l_Ng);
01334   // f
01335   outArgs.setSupports(OUT_ARG_f, epetraOutArgs.supports(EME::OUT_ARG_f));
01336   if (outArgs.supports(OUT_ARG_f)) {
01337     // W
01338     outArgs.setSupports(
01339       OUT_ARG_W, epetraOutArgs.supports(EME::OUT_ARG_W)&&!is_null(W_factory_));
01340     outArgs.setSupports(OUT_ARG_W_op,  epetraOutArgs.supports(EME::OUT_ARG_W));
01341     outArgs.set_W_properties(convert(epetraOutArgs.get_W_properties()));
01342     // DfDp
01343     for(int l=0; l<l_Np; ++l) {
01344       outArgs.setSupports(OUT_ARG_DfDp, l,
01345         convert(epetraOutArgs.supports(EME::OUT_ARG_DfDp, l)));
01346       if(!outArgs.supports(OUT_ARG_DfDp, l).none())
01347         outArgs.set_DfDp_properties(l,
01348           convert(epetraOutArgs.get_DfDp_properties(l)));
01349     }
01350   }
01351   // DgDx_dot and DgDx
01352   for(int j=0; j<l_Ng; ++j) {
01353     if (inArgs.supports(IN_ARG_x_dot))
01354       outArgs.setSupports(OUT_ARG_DgDx_dot, j,
01355         convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx_dot, j)));
01356     if(!outArgs.supports(OUT_ARG_DgDx_dot, j).none())
01357       outArgs.set_DgDx_dot_properties(j,
01358         convert(epetraOutArgs.get_DgDx_dot_properties(j)));
01359     if (inArgs.supports(IN_ARG_x))
01360       outArgs.setSupports(OUT_ARG_DgDx, j,
01361         convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx, j)));
01362     if(!outArgs.supports(OUT_ARG_DgDx, j).none())
01363       outArgs.set_DgDx_properties(j,
01364         convert(epetraOutArgs.get_DgDx_properties(j)));
01365   }
01366   // DgDp
01367   for(int j=0; j < l_Ng; ++j) for(int l=0; l < l_Np; ++l) {
01368     const EME::DerivativeSupport epetra_DgDp_j_l_support =
01369       epetraOutArgs.supports(EME::OUT_ARG_DgDp, j, l);
01370     outArgs.setSupports(OUT_ARG_DgDp, j, l,
01371       convert(epetra_DgDp_j_l_support));
01372     if(!outArgs.supports(OUT_ARG_DgDp, j, l).none())
01373       outArgs.set_DgDp_properties(j, l,
01374         convert(epetraOutArgs.get_DgDp_properties(j, l)));
01375   }
01376 #ifdef HAVE_THYRA_ME_POLYNOMIAL
01377   outArgs.setSupports(OUT_ARG_f_poly,
01378     epetraOutArgs.supports(EME::OUT_ARG_f_poly));
01379 #endif // HAVE_THYRA_ME_POLYNOMIAL
01380   prototypeOutArgs_ = outArgs;
01381 
01382   // We are current!
01383   currentInArgsOutArgs_ = true;
01384 
01385 }
01386 
01387 
01388 RCP<EpetraLinearOp>
01389 EpetraModelEvaluator::create_epetra_W_op() const
01390 {
01391   return Thyra::partialNonconstEpetraLinearOp(
01392     this->get_f_space(), this->get_x_space(),
01393     create_and_assert_W(*epetraModel_)
01394     );
01395 }
01396 
01397 
01398 } // namespace Thyra
01399 
01400 
01401 //
01402 // Non-member utility functions
01403 //
01404 
01405 
01406 Teuchos::RCP<Thyra::EpetraModelEvaluator>
01407 Thyra::epetraModelEvaluator(
01408   const RCP<const EpetraExt::ModelEvaluator> &epetraModel,
01409   const RCP<LinearOpWithSolveFactoryBase<double> > &W_factory
01410   )
01411 {
01412   return Teuchos::rcp(new EpetraModelEvaluator(epetraModel,W_factory));
01413 }
01414 
01415 
01416 Thyra::ModelEvaluatorBase::EDerivativeMultiVectorOrientation
01417 Thyra::convert(
01418   const EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation &mvOrientation
01419   )
01420 {
01421   switch(mvOrientation) {
01422     case EpetraExt::ModelEvaluator::DERIV_MV_BY_COL :
01423       return ModelEvaluatorBase::DERIV_MV_BY_COL;
01424     case EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW :
01425       return ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW;
01426     default:
01427       TEST_FOR_EXCEPT(true);
01428   }
01429   return ModelEvaluatorBase::DERIV_MV_BY_COL; // Should never be called!
01430 }
01431 
01432 
01433 EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation
01434 Thyra::convert(
01435   const ModelEvaluatorBase::EDerivativeMultiVectorOrientation &mvOrientation
01436   )
01437 {
01438   switch(mvOrientation) {
01439     case ModelEvaluatorBase::DERIV_MV_BY_COL :
01440       return EpetraExt::ModelEvaluator::DERIV_MV_BY_COL;
01441     case ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW :
01442       return EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW;
01443     default:
01444       TEST_FOR_EXCEPT(true);
01445   }
01446   return EpetraExt::ModelEvaluator::DERIV_MV_BY_COL; // Should never be called!
01447 }
01448 
01449 
01450 Thyra::ModelEvaluatorBase::DerivativeProperties
01451 Thyra::convert(
01452   const EpetraExt::ModelEvaluator::DerivativeProperties &derivativeProperties
01453   )
01454 {
01455   ModelEvaluatorBase::EDerivativeLinearity linearity;
01456   switch(derivativeProperties.linearity) {
01457     case EpetraExt::ModelEvaluator::DERIV_LINEARITY_UNKNOWN:
01458       linearity = ModelEvaluatorBase::DERIV_LINEARITY_UNKNOWN;
01459       break;
01460     case EpetraExt::ModelEvaluator::DERIV_LINEARITY_CONST:
01461       linearity = ModelEvaluatorBase::DERIV_LINEARITY_CONST;
01462       break;
01463     case  EpetraExt::ModelEvaluator::DERIV_LINEARITY_NONCONST:
01464       linearity = ModelEvaluatorBase::DERIV_LINEARITY_NONCONST;
01465       break;
01466     default:
01467       TEST_FOR_EXCEPT(true);
01468   }
01469   ModelEvaluatorBase::ERankStatus rank;
01470   switch(derivativeProperties.rank) {
01471     case EpetraExt::ModelEvaluator::DERIV_RANK_UNKNOWN:
01472       rank = ModelEvaluatorBase::DERIV_RANK_UNKNOWN;
01473       break;
01474     case EpetraExt::ModelEvaluator::DERIV_RANK_FULL:
01475       rank = ModelEvaluatorBase::DERIV_RANK_FULL;
01476       break;
01477     case EpetraExt::ModelEvaluator::DERIV_RANK_DEFICIENT:
01478       rank = ModelEvaluatorBase::DERIV_RANK_DEFICIENT;
01479       break;
01480     default:
01481       TEST_FOR_EXCEPT(true);
01482   }
01483   return ModelEvaluatorBase::DerivativeProperties(
01484     linearity,rank,derivativeProperties.supportsAdjoint);
01485 }
01486 
01487 
01488 Thyra::ModelEvaluatorBase::DerivativeSupport
01489 Thyra::convert(
01490   const EpetraExt::ModelEvaluator::DerivativeSupport &derivativeSupport
01491   )
01492 {
01493   ModelEvaluatorBase::DerivativeSupport ds;
01494   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_LINEAR_OP))
01495     ds.plus(ModelEvaluatorBase::DERIV_LINEAR_OP);
01496   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_MV_BY_COL))
01497     ds.plus(ModelEvaluatorBase::DERIV_MV_BY_COL);
01498   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW))
01499     ds.plus(ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW);
01500   return ds;
01501 }
01502 
01503 
01504 EpetraExt::ModelEvaluator::Derivative
01505 Thyra::convert(
01506   const ModelEvaluatorBase::Derivative<double> &derivative,
01507   const RCP<const Epetra_Map> &fnc_map,
01508   const RCP<const Epetra_Map> &var_map
01509   )
01510 {
01511   typedef ModelEvaluatorBase MEB;
01512   if(derivative.getLinearOp().get()) {
01513     return EpetraExt::ModelEvaluator::Derivative(
01514       Teuchos::rcp_const_cast<Epetra_Operator>(
01515         Teuchos::dyn_cast<const EpetraLinearOp>(*derivative.getLinearOp()).epetra_op()
01516         )
01517       );
01518   }
01519   else if(derivative.getDerivativeMultiVector().getMultiVector().get()) {
01520     return EpetraExt::ModelEvaluator::Derivative(
01521       EpetraExt::ModelEvaluator::DerivativeMultiVector(
01522         get_Epetra_MultiVector(
01523           ( derivative.getDerivativeMultiVector().getOrientation() == MEB::DERIV_MV_BY_COL
01524             ? *fnc_map
01525             : *var_map
01526             )
01527           ,derivative.getDerivativeMultiVector().getMultiVector()
01528           )
01529         ,convert(derivative.getDerivativeMultiVector().getOrientation())
01530         )
01531       );
01532   }
01533   return EpetraExt::ModelEvaluator::Derivative();
01534 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines