Thyra_EpetraModelEvaluator.cpp

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 Np = inArgs->Np();
00799   for( int l = 0; 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   using Teuchos::Polynomial;
00819 
00820   TEST_FOR_EXCEPT(0==epetraInArgs);
00821 
00822   RCP<const VectorBase<double> > x_dot;
00823   if( inArgs.supports(IN_ARG_x_dot) && (x_dot = inArgs.get_x_dot()).get() ) {
00824     RCP<const Epetra_Vector> e_x_dot = get_Epetra_Vector(*x_map_,x_dot);
00825     epetraInArgs->set_x_dot(e_x_dot);
00826   }
00827 
00828   RCP<const VectorBase<double> > x;
00829   if( inArgs.supports(IN_ARG_x) && (x = inArgs.get_x()).get() ) {
00830     RCP<const Epetra_Vector> e_x = get_Epetra_Vector(*x_map_,x);
00831     epetraInArgs->set_x(e_x);
00832   }
00833 
00834   RCP<const VectorBase<double> > p_l;
00835   for(int l = 0;  l < inArgs.Np(); ++l ) {
00836     p_l = inArgs.get_p(l);
00837     if(p_l.get()) epetraInArgs->set_p(l,get_Epetra_Vector(*p_map_[l],p_l));
00838   }
00839 
00840   RCP<const Polynomial< VectorBase<double> > > x_dot_poly;
00841   RCP<Epetra_Vector> epetra_ptr;
00842   if(
00843     inArgs.supports(IN_ARG_x_dot_poly)
00844     && (x_dot_poly = inArgs.get_x_dot_poly()).get()
00845     )
00846   {
00847     RCP<Polynomial<Epetra_Vector> > epetra_x_dot_poly = 
00848       rcp(new Polynomial<Epetra_Vector>(x_dot_poly->degree()));
00849     for (unsigned int i=0; i<=x_dot_poly->degree(); i++) {
00850       epetra_ptr = rcp_const_cast<Epetra_Vector>(
00851         get_Epetra_Vector(*x_map_, x_dot_poly->getCoefficient(i)) );
00852       epetra_x_dot_poly->setCoefficientPtr(i,epetra_ptr);
00853     }
00854     epetraInArgs->set_x_dot_poly(epetra_x_dot_poly);
00855   }
00856   
00857   RCP<const Polynomial< VectorBase<double> > > x_poly;
00858   if(
00859     inArgs.supports(IN_ARG_x_poly)
00860     && (x_poly = inArgs.get_x_poly()).get()
00861     )
00862   {
00863     RCP<Polynomial<Epetra_Vector> > epetra_x_poly = 
00864       rcp(new Polynomial<Epetra_Vector>(x_poly->degree()));
00865     for (unsigned int i=0; i<=x_poly->degree(); i++) {
00866       epetra_ptr = rcp_const_cast<Epetra_Vector>(
00867         get_Epetra_Vector(*x_map_, x_poly->getCoefficient(i)) );
00868       epetra_x_poly->setCoefficientPtr(i,epetra_ptr);
00869     }
00870     epetraInArgs->set_x_poly(epetra_x_poly);
00871   }
00872 
00873   if( inArgs.supports(IN_ARG_t) )
00874     epetraInArgs->set_t(inArgs.get_t());
00875   
00876   if( inArgs.supports(IN_ARG_alpha) )
00877     epetraInArgs->set_alpha(inArgs.get_alpha());
00878   
00879   if( inArgs.supports(IN_ARG_beta) )
00880     epetraInArgs->set_beta(inArgs.get_beta());
00881 
00882 }
00883 
00884 
00885 void EpetraModelEvaluator::convertOutArgsFromThyraToEpetra(
00886   const ModelEvaluatorBase::OutArgs<double> &outArgs,
00887   EpetraExt::ModelEvaluator::OutArgs *epetraUnscaledOutArgs_inout,
00888   RCP<LinearOpWithSolveBase<double> > *W_out,
00889   RCP<LinearOpBase<double> > *W_op_out,
00890   RCP<const LinearOpBase<double> > *fwdW_out,
00891   RCP<EpetraLinearOp> *efwdW_out,
00892   RCP<Epetra_Operator> *eW_out
00893   ) const
00894 {
00895 
00896   using Teuchos::rcp;
00897   using Teuchos::rcp_const_cast;
00898   using Teuchos::rcp_dynamic_cast;
00899   using Teuchos::OSTab;
00900   using Teuchos::implicit_cast;
00901   using Thyra::get_Epetra_Vector;
00902   typedef EpetraExt::ModelEvaluator EME;
00903 
00904   // Assert input
00905 #ifdef TEUCHOS_DEBUG
00906   TEUCHOS_ASSERT(epetraUnscaledOutArgs_inout);
00907   TEUCHOS_ASSERT(W_out);
00908   TEUCHOS_ASSERT(W_op_out);
00909   TEUCHOS_ASSERT(fwdW_out);
00910   TEUCHOS_ASSERT(efwdW_out);
00911   TEUCHOS_ASSERT(eW_out);
00912 #endif
00913 
00914   // Create easy to use references
00915   EpetraExt::ModelEvaluator::OutArgs &epetraUnscaledOutArgs = *epetraUnscaledOutArgs_inout;
00916   RCP<LinearOpWithSolveBase<double> > &W = *W_out;
00917   RCP<LinearOpBase<double> > &W_op = *W_op_out;
00918   RCP<const LinearOpBase<double> > &fwdW = *fwdW_out;
00919   RCP<EpetraLinearOp> &efwdW = *efwdW_out;
00920   RCP<Epetra_Operator> &eW = *eW_out;
00921 
00922   // f
00923   { 
00924     RCP<VectorBase<double> > f;
00925     if( outArgs.supports(OUT_ARG_f) && (f = outArgs.get_f()).get() )
00926       epetraUnscaledOutArgs.set_f(get_Epetra_Vector(*f_map_,f));
00927   }
00928     
00929   // g
00930   {
00931     RCP<VectorBase<double> > g_j;
00932     for(int j = 0;  j < outArgs.Ng(); ++j ) {
00933       g_j = outArgs.get_g(j);
00934       if(g_j.get()) epetraUnscaledOutArgs.set_g(j,get_Epetra_Vector(*g_map_[j],g_j));
00935     }
00936   }
00937   
00938   // W and W_op
00939   {
00940 
00941     if( outArgs.supports(OUT_ARG_W) && (W = outArgs.get_W()).get() ) {
00942       Thyra::uninitializeOp<double>(*W_factory_,&*W,&fwdW);
00943       if(fwdW.get()) {
00944         efwdW = rcp_const_cast<EpetraLinearOp>(
00945           rcp_dynamic_cast<const EpetraLinearOp>(fwdW,true));
00946       }
00947       else {
00948         efwdW = this->create_epetra_W_op();
00949         fwdW = efwdW;
00950       }
00951     }
00952     if( outArgs.supports(OUT_ARG_W_op) && (W_op = outArgs.get_W_op()).get() ) {
00953       if( W_op.get() && !efwdW.get() )
00954         efwdW = rcp_const_cast<EpetraLinearOp>(
00955           rcp_dynamic_cast<const EpetraLinearOp>(W_op,true));
00956     }
00957     if(efwdW.get()) {
00958       // By the time we get here, if we have an object in efwdW, then it
00959       // should already be embeadded with an underlying Epetra_Operator object
00960       // that was allocated by the EpetraExt::ModelEvaluator object.
00961       // Therefore, we should just have to grab this object and be on our way.
00962       eW = efwdW->epetra_op();
00963       epetraUnscaledOutArgs.set_W(eW);
00964     }
00965     // NOTE: Above, if both W and W_op are set and have been through at least
00966     // one prior evaluation (and therefore have Epetra_Operator objects embedded
00967     // in them), then we will use the Epetra_Operator embedded in W to pass to
00968     // the EpetraExt::ModelEvaluator object and ignore the Epetra_Operator
00969     // object in W_op.  In the standard use case, these will be the same
00970     // Epetra_Operator objects.  However, it is possible that the client could
00971     // use this interface in such a way that these would have different
00972     // Epetra_Operator objects embedded in them.  In this (very unlikely) case,
00973     // the Epetra_Operator embedded in W_op will be discarded!  This might be
00974     // surprising to a client but it is very unlikely that this will ever be a
00975     // problem, but the issue is duly noted here!  Only dangerous programming
00976     // use of this interface would cause any problem.
00977     
00978     // Note: The following derivative objects update in place!
00979 
00980   }
00981 
00982   // DfDp
00983   {
00984     Derivative<double> DfDp_l;
00985     for(int l = 0;  l < outArgs.Np(); ++l ) {
00986       if( !outArgs.supports(OUT_ARG_DfDp,l).none()
00987         && !(DfDp_l = outArgs.get_DfDp(l)).isEmpty() )
00988       {
00989         epetraUnscaledOutArgs.set_DfDp(l,convert(DfDp_l,f_map_,p_map_[l]));
00990       }
00991     }
00992   }
00993 
00994   // DgDx_dot
00995   {
00996     Derivative<double> DgDx_dot_j;
00997     for(int j = 0;  j < outArgs.Ng(); ++j ) {
00998       if( !outArgs.supports(OUT_ARG_DgDx_dot,j).none()
00999         && !(DgDx_dot_j = outArgs.get_DgDx_dot(j)).isEmpty() )
01000       {
01001         epetraUnscaledOutArgs.set_DgDx_dot(j,convert(DgDx_dot_j,g_map_[j],x_map_));
01002       }
01003     }
01004   }
01005 
01006   // DgDx
01007   {
01008     Derivative<double> DgDx_j;
01009     for(int j = 0;  j < outArgs.Ng(); ++j ) {
01010       if( !outArgs.supports(OUT_ARG_DgDx,j).none()
01011         && !(DgDx_j = outArgs.get_DgDx(j)).isEmpty() )
01012       {
01013         epetraUnscaledOutArgs.set_DgDx(j,convert(DgDx_j,g_map_[j],x_map_));
01014       }
01015     }
01016   }
01017 
01018   // DgDp
01019   {
01020     DerivativeSupport DgDp_j_l_support;
01021     Derivative<double> DgDp_j_l;
01022     for (int j = 0;  j < outArgs.Ng(); ++j ) {
01023       for (int l = 0;  l < outArgs.Np(); ++l ) {
01024         if (!(DgDp_j_l_support = outArgs.supports(OUT_ARG_DgDp,j,l)).none()
01025           && !(DgDp_j_l = outArgs.get_DgDp(j,l)).isEmpty() )
01026         {
01027           epetraUnscaledOutArgs.set_DgDp(j,l,convert(DgDp_j_l,g_map_[j],p_map_[l]));
01028         }
01029       }
01030     }
01031   }
01032 
01033   // f_poly
01034   RCP<const Teuchos::Polynomial< VectorBase<double> > > f_poly;
01035   if( outArgs.supports(OUT_ARG_f_poly) && (f_poly = outArgs.get_f_poly()).get() )
01036   {
01037     RCP<Teuchos::Polynomial<Epetra_Vector> > epetra_f_poly = 
01038       Teuchos::rcp(new Teuchos::Polynomial<Epetra_Vector>(f_poly->degree()));
01039     for (unsigned int i=0; i<=f_poly->degree(); i++) {
01040       RCP<Epetra_Vector> epetra_ptr
01041         = Teuchos::rcp_const_cast<Epetra_Vector>(get_Epetra_Vector(*f_map_,
01042             f_poly->getCoefficient(i)));
01043       epetra_f_poly->setCoefficientPtr(i,epetra_ptr);
01044     }
01045     epetraUnscaledOutArgs.set_f_poly(epetra_f_poly);
01046   }
01047 
01048 }
01049 
01050 
01051 void EpetraModelEvaluator::preEvalScalingSetup(
01052   EpetraExt::ModelEvaluator::InArgs *epetraInArgs_inout,
01053   EpetraExt::ModelEvaluator::OutArgs *epetraUnscaledOutArgs_inout,
01054   const RCP<Teuchos::FancyOStream> &out,
01055   const Teuchos::EVerbosityLevel verbLevel
01056   ) const
01057 {
01058   
01059   typedef EpetraExt::ModelEvaluator EME;
01060   
01061 #ifdef TEUCHOS_DEBUG
01062   TEUCHOS_ASSERT(epetraInArgs_inout);
01063   TEUCHOS_ASSERT(epetraUnscaledOutArgs_inout);
01064 #endif
01065 
01066   EpetraExt::ModelEvaluator::InArgs
01067     &epetraInArgs = *epetraInArgs_inout;
01068   EpetraExt::ModelEvaluator::OutArgs
01069     &epetraUnscaledOutArgs = *epetraUnscaledOutArgs_inout;
01070 
01071   if (
01072     ( stateFunctionScaling_ == STATE_FUNC_SCALING_ROW_SUM )
01073     &&
01074     (
01075       epetraUnscaledOutArgs.supports(EME::OUT_ARG_f) 
01076       &&
01077       epetraUnscaledOutArgs.funcOrDerivesAreSet(EME::OUT_ARG_f)
01078       )
01079     &&
01080     (
01081       epetraUnscaledOutArgs.supports(EME::OUT_ARG_W)
01082       &&
01083       is_null(epetraUnscaledOutArgs.get_W())
01084       )
01085     )
01086   {
01087     // This is the first pass through with scaling turned on and the client
01088     // turned on automatic scaling but did not ask for W.  We must compute W
01089     // in order to compute the scale factors so we must allocate a temporary W
01090     // just to compute the scale factors and then throw it away.  If the
01091     // client wants to evaluate W at the same point, then it should have
01092     // passed W in but that is not our problem here.  The ModelEvaluator
01093     // relies on the client to set up the calls to allow for efficient
01094     // evaluation.
01095 
01096     if(out.get() && verbLevel >= Teuchos::VERB_LOW)
01097       *out
01098         << "\nCreating a temporary Epetra W to compute scale factors"
01099         << " for f(...) ...\n";
01100     epetraUnscaledOutArgs.set_W(create_and_assert_W(*epetraModel_));
01101     if( epetraInArgs.supports(EME::IN_ARG_beta) )
01102       epetraInArgs.set_beta(1.0);
01103     if( epetraInArgs.supports(EME::IN_ARG_alpha) )
01104       epetraInArgs.set_alpha(0.0);
01105   }
01106   
01107 }
01108 
01109 
01110 void EpetraModelEvaluator::postEvalScalingSetup(
01111   const EpetraExt::ModelEvaluator::OutArgs &epetraUnscaledOutArgs,
01112   const RCP<Teuchos::FancyOStream> &out,
01113   const Teuchos::EVerbosityLevel verbLevel
01114   ) const
01115 {
01116 
01117   using Teuchos::OSTab;
01118   using Teuchos::rcp;
01119   using Teuchos::rcp_const_cast;
01120   using Teuchos::includesVerbLevel;
01121 
01122   // Compute the scale factors for the state function f(...)
01123   switch(stateFunctionScaling_) {
01124 
01125     case STATE_FUNC_SCALING_ROW_SUM: {
01126 
01127       // Compute the inverse row-sum scaling from W
01128 
01129       const RCP<Epetra_RowMatrix>
01130         ermW = get_Epetra_RowMatrix(epetraUnscaledOutArgs);
01131       // Note: Above, we get the Epetra W object directly from the Epetra
01132       // OutArgs object since this might be a temporary matrix just to
01133       // compute scaling factors.  In this case, the stack funtion variable
01134       // eW might be empty!
01135 
01136       RCP<Epetra_Vector>
01137         invRowSums = rcp(new Epetra_Vector(ermW->OperatorRangeMap()));
01138       // Above: From the documentation is seems that the RangeMap should be
01139       // okay but who knows for sure!
01140 
01141       ermW->InvRowSums(*invRowSums);
01142 
01143       if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW)) {
01144         *out
01145           << "\nComputed inverse row sum scaling from W that"
01146           " will be used to scale f(...) and its derivatives:\n";
01147         double minVal = 0, maxVal = 0, avgVal = 0;
01148         invRowSums->MinValue(&minVal);
01149         invRowSums->MaxValue(&maxVal);
01150         invRowSums->MeanValue(&avgVal);
01151         OSTab tab(out);
01152         *out
01153           << "min(invRowSums) = " << minVal << "\n"
01154           << "max(invRowSums) = " << maxVal << "\n"
01155           << "avg(invRowSums) = " << avgVal << "\n";
01156       }
01157 
01158       stateFunctionScalingVec_ = invRowSums;
01159 
01160       break;
01161 
01162     }
01163 
01164     default:
01165       TEST_FOR_EXCEPT("Should never get here!");
01166 
01167   }
01168 
01169   epetraOutArgsScaling_ = epetraModel_->createOutArgs();
01170 
01171   epetraOutArgsScaling_.set_f(
01172     rcp_const_cast<Epetra_Vector>(stateFunctionScalingVec_) );
01173 
01174 }
01175 
01176 
01177 void EpetraModelEvaluator::finishConvertingOutArgsFromEpetraToThyra(
01178   const EpetraExt::ModelEvaluator::OutArgs &epetraOutArgs,
01179   RCP<LinearOpWithSolveBase<double> > &W,
01180   RCP<LinearOpBase<double> > &W_op,
01181   RCP<const LinearOpBase<double> > &fwdW,
01182   RCP<EpetraLinearOp> &efwdW,
01183   RCP<Epetra_Operator> &eW,
01184   const ModelEvaluatorBase::OutArgs<double> &outArgs
01185   ) const
01186 {
01187 
01188   using Teuchos::rcp_dynamic_cast;
01189   typedef EpetraExt::ModelEvaluator EME;
01190 
01191   if(efwdW.get()) {
01192     efwdW->setFullyInitialized(true); 
01193     // NOTE: Above will directly update W_op also if W.get()==NULL!
01194   }
01195   
01196   if( W.get() ) {
01197     Thyra::initializeOp<double>(*W_factory_,fwdW,&*W);
01198     W->setOStream(this->getOStream());
01199   }
01200 
01201   if( W_op.get() ) {
01202     if( W_op.shares_resource(efwdW) ) {
01203       // W_op was already updated above since *efwdW is the same object as *W_op
01204     }
01205     else {
01206       rcp_dynamic_cast<EpetraLinearOp>(W_op,true)->setFullyInitialized(true);
01207     }
01208   }
01209 
01210 }
01211 
01212 
01213 void EpetraModelEvaluator::updateNominalValuesAndBounds() const
01214 {
01215 
01216   using Teuchos::rcp;
01217   using Teuchos::implicit_cast;
01218   typedef ModelEvaluatorBase MEB;
01219   typedef EpetraExt::ModelEvaluator EME;
01220 
01221   if( !nominalValuesAndBoundsAreUpdated_ ) {
01222 
01223     // Gather the nominal values and bounds into Epetra InArgs objects
01224 
01225     EME::InArgs epetraOrigNominalValues;
01226     EpetraExt::gatherModelNominalValues(
01227       *epetraModel_, &epetraOrigNominalValues );
01228 
01229     EME::InArgs epetraOrigLowerBounds;
01230     EME::InArgs epetraOrigUpperBounds;
01231     EpetraExt::gatherModelBounds(
01232       *epetraModel_, &epetraOrigLowerBounds, &epetraOrigUpperBounds );
01233 
01234     // Set up Epetra InArgs scaling object
01235 
01236     epetraInArgsScaling_ = epetraModel_->createInArgs();
01237 
01238     if( !is_null(stateVariableScalingVec_) ) {
01239       invStateVariableScalingVec_
01240         = EpetraExt::createInverseModelScalingVector(stateVariableScalingVec_);
01241       if( epetraOrigNominalValues.supports(EME::IN_ARG_x_dot) ) {
01242         epetraInArgsScaling_.set_x_dot(invStateVariableScalingVec_);
01243       }
01244       if( epetraOrigNominalValues.supports(EME::IN_ARG_x) ) {
01245         epetraInArgsScaling_.set_x(invStateVariableScalingVec_);
01246       }
01247     }
01248     
01249     // Scale the original variables and bounds
01250 
01251     EME::InArgs epetraScaledNominalValues = epetraModel_->createInArgs();
01252     EpetraExt::scaleModelVars(
01253       epetraOrigNominalValues, epetraInArgsScaling_, &epetraScaledNominalValues
01254       );
01255 
01256     EME::InArgs epetraScaledLowerBounds = epetraModel_->createInArgs();
01257     EME::InArgs epetraScaledUpperBounds = epetraModel_->createInArgs();
01258     EpetraExt::scaleModelBounds(
01259       epetraOrigLowerBounds, epetraOrigUpperBounds, epetraModel_->getInfBound(),
01260       epetraInArgsScaling_,
01261       &epetraScaledLowerBounds, &epetraScaledUpperBounds
01262       );
01263 
01264     // Wrap the scaled epetra InArgs objects as Thyra InArgs objects!
01265 
01266     nominalValues_ = this->createInArgs();
01267     lowerBounds_ = this->createInArgs();
01268     upperBounds_ = this->createInArgs();
01269     convertInArgsFromEpetraToThyra(epetraScaledNominalValues, &nominalValues_);
01270     convertInArgsFromEpetraToThyra(epetraScaledLowerBounds, &lowerBounds_);
01271     convertInArgsFromEpetraToThyra(epetraScaledUpperBounds, &upperBounds_);
01272 
01273     nominalValuesAndBoundsAreUpdated_ = true;
01274 
01275   }
01276   else {
01277 
01278     // The nominal values and bounds should already be updated an should have
01279     // the currect scaling!
01280 
01281   }
01282 
01283 }
01284 
01285 
01286 void EpetraModelEvaluator::updateInArgsOutArgs() const
01287 {
01288 
01289   typedef EpetraExt::ModelEvaluator EME;
01290 
01291   const EpetraExt::ModelEvaluator &epetraModel = *epetraModel_;
01292   EME::InArgs  epetraInArgs  = epetraModel.createInArgs();
01293   EME::OutArgs epetraOutArgs = epetraModel.createOutArgs();
01294   const int Np = epetraOutArgs.Np();
01295   const int Ng = epetraOutArgs.Ng();
01296 
01297   //
01298   // InArgs
01299   //
01300 
01301   InArgsSetup<double> inArgs;
01302   inArgs.setModelEvalDescription(this->description());
01303   inArgs.set_Np(epetraInArgs.Np());
01304   inArgs.setSupports(IN_ARG_x_dot, epetraInArgs.supports(EME::IN_ARG_x_dot));
01305   inArgs.setSupports(IN_ARG_x, epetraInArgs.supports(EME::IN_ARG_x));
01306   inArgs.setSupports(IN_ARG_x_dot_poly,
01307     epetraInArgs.supports(EME::IN_ARG_x_dot_poly));
01308   inArgs.setSupports(IN_ARG_x_poly, epetraInArgs.supports(EME::IN_ARG_x_poly));
01309   inArgs.setSupports(IN_ARG_t, epetraInArgs.supports(EME::IN_ARG_t));
01310   inArgs.setSupports(IN_ARG_alpha, epetraInArgs.supports(EME::IN_ARG_alpha));
01311   inArgs.setSupports(IN_ARG_beta, epetraInArgs.supports(EME::IN_ARG_beta));
01312   prototypeInArgs_ = inArgs;
01313 
01314   //
01315   // OutArgs
01316   //
01317 
01318   OutArgsSetup<double> outArgs;
01319   outArgs.setModelEvalDescription(this->description());
01320   outArgs.set_Np_Ng(Np, Ng);
01321   // f
01322   outArgs.setSupports(OUT_ARG_f, epetraOutArgs.supports(EME::OUT_ARG_f));
01323   if (outArgs.supports(OUT_ARG_f)) {
01324     // W
01325     outArgs.setSupports(
01326       OUT_ARG_W, epetraOutArgs.supports(EME::OUT_ARG_W)&&!is_null(W_factory_));
01327     outArgs.setSupports(OUT_ARG_W_op,  epetraOutArgs.supports(EME::OUT_ARG_W));
01328     outArgs.set_W_properties(convert(epetraOutArgs.get_W_properties()));
01329     // DfDp
01330     for(int l=0; l<Np; ++l) {
01331       outArgs.setSupports(OUT_ARG_DfDp, l,
01332         convert(epetraOutArgs.supports(EME::OUT_ARG_DfDp, l)));
01333       if(!outArgs.supports(OUT_ARG_DfDp, l).none())
01334         outArgs.set_DfDp_properties(l,
01335           convert(epetraOutArgs.get_DfDp_properties(l)));
01336     }
01337   }
01338   // DgDx_dot and DgDx
01339   for(int j=0; j<Ng; ++j) {
01340     if (inArgs.supports(IN_ARG_x_dot))
01341       outArgs.setSupports(OUT_ARG_DgDx_dot, j,
01342         convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx_dot, j)));
01343     if(!outArgs.supports(OUT_ARG_DgDx_dot, j).none())
01344       outArgs.set_DgDx_dot_properties(j,
01345         convert(epetraOutArgs.get_DgDx_dot_properties(j)));
01346     if (inArgs.supports(IN_ARG_x))
01347       outArgs.setSupports(OUT_ARG_DgDx, j,
01348         convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx, j)));
01349     if(!outArgs.supports(OUT_ARG_DgDx, j).none())
01350       outArgs.set_DgDx_properties(j,
01351         convert(epetraOutArgs.get_DgDx_properties(j)));
01352   }
01353   // DgDp
01354   for(int j=0; j<Ng; ++j) for(int l=0; l<Np; ++l) {
01355     const EME::DerivativeSupport epetra_DgDp_j_l_support =
01356       epetraOutArgs.supports(EME::OUT_ARG_DgDp, j, l);
01357     outArgs.setSupports(OUT_ARG_DgDp, j, l,
01358       convert(epetra_DgDp_j_l_support));
01359     if(!outArgs.supports(OUT_ARG_DgDp, j, l).none())
01360       outArgs.set_DgDp_properties(j, l,
01361         convert(epetraOutArgs.get_DgDp_properties(j, l)));
01362   }
01363   outArgs.setSupports(OUT_ARG_f_poly,
01364     epetraOutArgs.supports(EME::OUT_ARG_f_poly));
01365   prototypeOutArgs_ = outArgs;
01366 
01367   // We are current!
01368   currentInArgsOutArgs_ = true;
01369 
01370 }
01371 
01372 
01373 RCP<EpetraLinearOp>
01374 EpetraModelEvaluator::create_epetra_W_op() const
01375 {
01376   return Thyra::partialNonconstEpetraLinearOp(
01377     this->get_f_space(), this->get_x_space(),
01378     create_and_assert_W(*epetraModel_)
01379     );
01380 }
01381 
01382 
01383 } // namespace Thyra
01384 
01385 
01386 //
01387 // Non-member utility functions
01388 //
01389 
01390 
01391 Teuchos::RCP<Thyra::EpetraModelEvaluator>
01392 Thyra::epetraModelEvaluator(
01393   const RCP<const EpetraExt::ModelEvaluator> &epetraModel,
01394   const RCP<LinearOpWithSolveFactoryBase<double> > &W_factory
01395   )
01396 {
01397   return Teuchos::rcp(new EpetraModelEvaluator(epetraModel,W_factory));
01398 }
01399 
01400 
01401 Thyra::ModelEvaluatorBase::EDerivativeMultiVectorOrientation
01402 Thyra::convert(
01403   const EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation &mvOrientation
01404   )
01405 {
01406   switch(mvOrientation) {
01407     case EpetraExt::ModelEvaluator::DERIV_MV_BY_COL :
01408       return ModelEvaluatorBase::DERIV_MV_BY_COL;
01409     case EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW :
01410       return ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW;
01411     default:
01412       TEST_FOR_EXCEPT(true);
01413   }
01414   return ModelEvaluatorBase::DERIV_MV_BY_COL; // Should never be called!
01415 }
01416 
01417 
01418 EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation
01419 Thyra::convert(
01420   const ModelEvaluatorBase::EDerivativeMultiVectorOrientation &mvOrientation
01421   )
01422 {
01423   switch(mvOrientation) {
01424     case ModelEvaluatorBase::DERIV_MV_BY_COL :
01425       return EpetraExt::ModelEvaluator::DERIV_MV_BY_COL;
01426     case ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW :
01427       return EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW;
01428     default:
01429       TEST_FOR_EXCEPT(true);
01430   }
01431   return EpetraExt::ModelEvaluator::DERIV_MV_BY_COL; // Should never be called!
01432 }
01433 
01434 
01435 Thyra::ModelEvaluatorBase::DerivativeProperties
01436 Thyra::convert(
01437   const EpetraExt::ModelEvaluator::DerivativeProperties &derivativeProperties
01438   )
01439 {
01440   ModelEvaluatorBase::EDerivativeLinearity linearity;
01441   switch(derivativeProperties.linearity) {
01442     case EpetraExt::ModelEvaluator::DERIV_LINEARITY_UNKNOWN:
01443       linearity = ModelEvaluatorBase::DERIV_LINEARITY_UNKNOWN;
01444       break;
01445     case EpetraExt::ModelEvaluator::DERIV_LINEARITY_CONST:
01446       linearity = ModelEvaluatorBase::DERIV_LINEARITY_CONST;
01447       break;
01448     case  EpetraExt::ModelEvaluator::DERIV_LINEARITY_NONCONST:
01449       linearity = ModelEvaluatorBase::DERIV_LINEARITY_NONCONST;
01450       break;
01451     default:
01452       TEST_FOR_EXCEPT(true);
01453   }
01454   ModelEvaluatorBase::ERankStatus rank;
01455   switch(derivativeProperties.rank) {
01456     case EpetraExt::ModelEvaluator::DERIV_RANK_UNKNOWN:
01457       rank = ModelEvaluatorBase::DERIV_RANK_UNKNOWN;
01458       break;
01459     case EpetraExt::ModelEvaluator::DERIV_RANK_FULL:
01460       rank = ModelEvaluatorBase::DERIV_RANK_FULL;
01461       break;
01462     case EpetraExt::ModelEvaluator::DERIV_RANK_DEFICIENT:
01463       rank = ModelEvaluatorBase::DERIV_RANK_DEFICIENT;
01464       break;
01465     default:
01466       TEST_FOR_EXCEPT(true);
01467   }
01468   return ModelEvaluatorBase::DerivativeProperties(
01469     linearity,rank,derivativeProperties.supportsAdjoint);
01470 }
01471 
01472 
01473 Thyra::ModelEvaluatorBase::DerivativeSupport
01474 Thyra::convert(
01475   const EpetraExt::ModelEvaluator::DerivativeSupport &derivativeSupport
01476   )
01477 {
01478   ModelEvaluatorBase::DerivativeSupport ds;
01479   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_LINEAR_OP))
01480     ds.plus(ModelEvaluatorBase::DERIV_LINEAR_OP);
01481   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_MV_BY_COL))
01482     ds.plus(ModelEvaluatorBase::DERIV_MV_BY_COL);
01483   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW))
01484     ds.plus(ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW);
01485   return ds;
01486 }
01487 
01488 
01489 EpetraExt::ModelEvaluator::Derivative
01490 Thyra::convert(
01491   const ModelEvaluatorBase::Derivative<double> &derivative,
01492   const RCP<const Epetra_Map> &fnc_map,
01493   const RCP<const Epetra_Map> &var_map
01494   )
01495 {
01496   typedef ModelEvaluatorBase MEB;
01497   if(derivative.getLinearOp().get()) {
01498     return EpetraExt::ModelEvaluator::Derivative(
01499       Teuchos::rcp_const_cast<Epetra_Operator>(
01500         Teuchos::dyn_cast<const EpetraLinearOp>(*derivative.getLinearOp()).epetra_op()
01501         )
01502       );
01503   }
01504   else if(derivative.getDerivativeMultiVector().getMultiVector().get()) {
01505     return EpetraExt::ModelEvaluator::Derivative(
01506       EpetraExt::ModelEvaluator::DerivativeMultiVector(
01507         get_Epetra_MultiVector(
01508           ( derivative.getDerivativeMultiVector().getOrientation() == MEB::DERIV_MV_BY_COL
01509             ? *fnc_map
01510             : *var_map
01511             )
01512           ,derivative.getDerivativeMultiVector().getMultiVector()
01513           )
01514         ,convert(derivative.getDerivativeMultiVector().getOrientation())
01515         )
01516       );
01517   }
01518   return EpetraExt::ModelEvaluator::Derivative();
01519 }

Generated on Thu Feb 11 16:28:43 2010 for EpetraExt/Thyra Adapters by  doxygen 1.4.7