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   static RCP<const Teuchos::ParameterList> validPL;
00347   if(is_null(validPL)) {
00348     RCP<Teuchos::ParameterList>
00349       pl = Teuchos::rcp(new Teuchos::ParameterList());
00350     stateFunctionScalingValidator = rcp(
00351       new StringToIntegralParameterEntryValidator<EStateFunctionScaling>(
00352         tuple<std::string>(
00353           "None",
00354           "Row Sum"
00355           ),
00356         tuple<std::string>(
00357           "Do not scale the state function f(...) in this class.",
00358 
00359           "Scale the state function f(...) and all its derivatives\n"
00360           "using the row sum scaling from the initial Jacobian\n"
00361           "W=d(f)/d(x).  Note, this only works with Epetra_CrsMatrix\n"
00362           "currently."
00363           ),
00364         tuple<EStateFunctionScaling>(
00365           STATE_FUNC_SCALING_NONE,
00366           STATE_FUNC_SCALING_ROW_SUM
00367           ),
00368         StateFunctionScaling_name
00369         )
00370       );
00371     pl->set(StateFunctionScaling_name,StateFunctionScaling_default,
00372       "Determines if and how the state function f(...) and all of its\n"
00373       "derivatives are scaled.  The scaling is done explicitly so there should\n"
00374       "be no impact on the meaning of inner products or tolerances for\n"
00375       "linear solves.",
00376       stateFunctionScalingValidator
00377       );
00378     Teuchos::setupVerboseObjectSublist(&*pl);
00379     validPL = pl;
00380   }
00381   return validPL;
00382 }
00383 
00384 
00385 // Overridden from ModelEvaulator.
00386 
00387 
00388 int EpetraModelEvaluator::Np() const
00389 {
00390   return p_space_.size();
00391 }
00392 
00393 
00394 int EpetraModelEvaluator::Ng() const
00395 {
00396   return g_space_.size();
00397 }
00398 
00399 
00400 RCP<const VectorSpaceBase<double> >
00401 EpetraModelEvaluator::get_x_space() const
00402 {
00403   return x_space_;
00404 }
00405 
00406 
00407 RCP<const VectorSpaceBase<double> >
00408 EpetraModelEvaluator::get_f_space() const
00409 {
00410   return f_space_;
00411 }
00412 
00413 
00414 RCP<const VectorSpaceBase<double> >
00415 EpetraModelEvaluator::get_p_space(int l) const
00416 {
00417 #ifdef TEUCHOS_DEBUG
00418   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, this->Np() );
00419 #endif
00420   return p_space_[l];
00421 }
00422 
00423 
00424 RCP<const Teuchos::Array<std::string> >
00425 EpetraModelEvaluator::get_p_names(int l) const
00426 {
00427 #ifdef TEUCHOS_DEBUG
00428   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, this->Np() );
00429 #endif
00430   return epetraModel_->get_p_names(l);
00431 }
00432 
00433 
00434 RCP<const VectorSpaceBase<double> >
00435 EpetraModelEvaluator::get_g_space(int j) const
00436 {
00437   TEST_FOR_EXCEPT( ! ( 0 <= j && j < this->Ng() ) );
00438   return g_space_[j];
00439 }
00440 
00441 
00442 ModelEvaluatorBase::InArgs<double>
00443 EpetraModelEvaluator::getNominalValues() const
00444 {
00445   updateNominalValuesAndBounds();
00446   return nominalValues_;
00447 }
00448 
00449 
00450 ModelEvaluatorBase::InArgs<double>
00451 EpetraModelEvaluator::getLowerBounds() const
00452 {
00453   updateNominalValuesAndBounds();
00454   return lowerBounds_;
00455 }
00456 
00457 
00458 ModelEvaluatorBase::InArgs<double>
00459 EpetraModelEvaluator::getUpperBounds() const
00460 {
00461   updateNominalValuesAndBounds();
00462   return upperBounds_;
00463 }
00464 
00465 
00466 RCP<LinearOpWithSolveBase<double> >
00467 EpetraModelEvaluator::create_W() const
00468 {
00469   TEST_FOR_EXCEPTION(
00470     W_factory_.get()==NULL, std::logic_error
00471     ,"Thyra::EpetraModelEvaluator::create_W(): "
00472     "Error, the client did not set a LinearOpWithSolveFactoryBase"
00473     " object for W!"
00474     );
00475   W_factory_->setOStream(this->getOStream());
00476   W_factory_->setVerbLevel(this->getVerbLevel());
00477   return W_factory_->createOp();
00478 }
00479 
00480 
00481 RCP<LinearOpBase<double> >
00482 EpetraModelEvaluator::create_W_op() const
00483 {
00484   return this->create_epetra_W_op();
00485 }
00486 
00487 
00488 RCP<const LinearOpWithSolveFactoryBase<double> >
00489 EpetraModelEvaluator::get_W_factory() const
00490 {
00491   return W_factory_;
00492 }
00493 
00494 
00495 ModelEvaluatorBase::InArgs<double> EpetraModelEvaluator::createInArgs() const
00496 {
00497   if (!currentInArgsOutArgs_)
00498     updateInArgsOutArgs();
00499   return prototypeInArgs_;
00500 }
00501 
00502 
00503 void EpetraModelEvaluator::reportFinalPoint(
00504   const ModelEvaluatorBase::InArgs<double> &finalPoint,
00505   const bool wasSolved
00506   )
00507 {
00508   finalPoint_ = this->createInArgs();
00509   finalPoint_.setArgs(finalPoint);
00510   finalPointWasSolved_ = wasSolved;
00511 }
00512 
00513 
00514 // Private functions overridden from ModelEvaulatorDefaultBase
00515 
00516 
00517 RCP<LinearOpBase<double> >
00518 EpetraModelEvaluator::create_DfDp_op_impl(int l) const
00519 {
00520   TEST_FOR_EXCEPT(true);
00521   return Teuchos::null;
00522 }
00523 
00524 
00525 RCP<LinearOpBase<double> >
00526 EpetraModelEvaluator::create_DgDx_dot_op_impl(int j) const
00527 {
00528   TEST_FOR_EXCEPT(true);
00529   return Teuchos::null;
00530 }
00531 
00532 
00533 RCP<LinearOpBase<double> >
00534 EpetraModelEvaluator::create_DgDx_op_impl(int j) const
00535 {
00536   TEST_FOR_EXCEPT(true);
00537   return Teuchos::null;
00538 }
00539 
00540 
00541 RCP<LinearOpBase<double> >
00542 EpetraModelEvaluator::create_DgDp_op_impl( int j, int l ) const
00543 {
00544   TEST_FOR_EXCEPT(true);
00545   return Teuchos::null;
00546 }
00547 
00548 
00549 ModelEvaluatorBase::OutArgs<double>
00550 EpetraModelEvaluator::createOutArgsImpl() const
00551 {
00552   if (!currentInArgsOutArgs_)
00553     updateInArgsOutArgs();
00554   return prototypeOutArgs_;
00555 }
00556 
00557 
00558 void EpetraModelEvaluator::evalModelImpl(
00559   const ModelEvaluatorBase::InArgs<double>& inArgs_in,
00560   const ModelEvaluatorBase::OutArgs<double>& outArgs
00561   ) const
00562 {
00563 
00564   using Teuchos::rcp;
00565   using Teuchos::rcp_const_cast;
00566   using Teuchos::rcp_dynamic_cast;
00567   using Teuchos::OSTab;
00568   using Teuchos::includesVerbLevel;
00569   typedef EpetraExt::ModelEvaluator EME;
00570 
00571   //
00572   // A) Initial setup
00573   //
00574 
00575   // Make sure that we are fully initialized!
00576   this->updateNominalValuesAndBounds();
00577 
00578   // Make sure we grab the initial guess first!
00579   InArgs<double> inArgs = this->getNominalValues();
00580   // Now, copy the parameters from the input inArgs_in object to the inArgs
00581   // object.  Any input objects that are set in inArgs_in will overwrite those
00582   // in inArgs that will already contain the nominal values.  This will insure
00583   // that all input parameters are set and those that are not set by the
00584   // client will be at their nominal values (as defined by the underlying
00585   // EpetraExt::ModelEvaluator object).  The full set of Thyra input arguments
00586   // must be set before these can be translated into Epetra input arguments.
00587   inArgs.setArgs(inArgs_in);
00588 
00589   // Print the header and the values of the inArgs and outArgs objects!
00590   typedef double Scalar; // Needed for below macro!
00591   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00592     "Thyra::EpetraModelEvaluator",inArgs,outArgs,Teuchos::null
00593     );
00594 
00595   // State function Scaling
00596   const bool firstTimeStateFuncScaling
00597     = (
00598       stateFunctionScaling_ != STATE_FUNC_SCALING_NONE
00599       && is_null(stateFunctionScalingVec_)
00600       );
00601   
00602   typedef Teuchos::VerboseObjectTempState<LinearOpWithSolveFactoryBase<double> > VOTSLOWSF;
00603   VOTSLOWSF W_factory_outputTempState(W_factory_,out,verbLevel);
00604 
00605   Teuchos::Time timer("");
00606 
00607   //
00608   // B) Prepressess the InArgs and OutArgs in preparation to call
00609   // the underlying EpetraExt::ModelEvaluator
00610   //
00611 
00612   //
00613   // B.1) Translate InArgs from Thyra to Epetra objects
00614   //
00615   
00616   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00617     *out << "\nSetting-up/creating input arguments ...\n";
00618   timer.start(true);
00619 
00620   // Unwrap input Thyra objects to get Epetra objects
00621   EME::InArgs epetraScaledInArgs = epetraModel_->createInArgs();
00622   convertInArgsFromThyraToEpetra( inArgs, &epetraScaledInArgs );
00623 
00624   // Unscale the input Epetra objects which will be passed to the underlying
00625   // EpetraExt::ModelEvaluator object.
00626   EME::InArgs epetraInArgs = epetraModel_->createInArgs();
00627   EpetraExt::unscaleModelVars(
00628     epetraScaledInArgs, epetraInArgsScaling_, &epetraInArgs,
00629     out.get(), verbLevel
00630     );
00631 
00632   timer.stop();
00633   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00634     OSTab(out).o() << "\nTime to setup InArgs = "<<timer.totalElapsedTime()<<" sec\n";
00635 
00636   //
00637   // B.2) Convert from Thyra to Epetra OutArgs
00638   //
00639   
00640   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00641     *out << "\nSetting-up/creating output arguments ...\n";
00642   timer.start(true);
00643   
00644   // The unscaled Epetra OutArgs that will be passed to the
00645   // underlying EpetraExt::ModelEvaluator object
00646   EME::OutArgs epetraUnscaledOutArgs = epetraModel_->createOutArgs();
00647 
00648   // Various objects that are needed later (see documentation in
00649   // the function convertOutArgsFromThyraToEpetra(...)
00650   RCP<LinearOpWithSolveBase<double> > W;
00651   RCP<LinearOpBase<double> > W_op;
00652   RCP<const LinearOpBase<double> > fwdW;
00653   RCP<EpetraLinearOp> efwdW;
00654   RCP<Epetra_Operator> eW;
00655   
00656   // Convert from Thyra to Epetra OutArgs and grap some of the intermediate
00657   // objects accessed along the way that are needed later.
00658   convertOutArgsFromThyraToEpetra(
00659     outArgs,
00660     &epetraUnscaledOutArgs,
00661     &W, &W_op, &fwdW, &efwdW, &eW
00662     );
00663 
00664   //
00665   // B.3) Setup OutArgs to computing scaling if needed
00666   //
00667 
00668   if (firstTimeStateFuncScaling) {
00669     preEvalScalingSetup(&epetraInArgs,&epetraUnscaledOutArgs,out,verbLevel);
00670   }
00671 
00672   timer.stop();
00673   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00674     OSTab(out).o()
00675       << "\nTime to setup OutArgs = "
00676       << timer.totalElapsedTime() <<" sec\n";
00677 
00678   //
00679   // C) Evaluate the underlying EpetraExt model to compute the Epetra outputs
00680   //
00681   
00682   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00683     *out << "\nEvaluating the Epetra output functions ...\n";
00684   timer.start(true);
00685 
00686   epetraModel_->evalModel(epetraInArgs,epetraUnscaledOutArgs);
00687 
00688   timer.stop();
00689   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00690     OSTab(out).o()
00691       << "\nTime to evaluate Epetra output functions = "
00692       << timer.totalElapsedTime() <<" sec\n";
00693 
00694   //
00695   // D) Postprocess the output objects
00696   //
00697 
00698   //
00699   // D.1) Compute the scaling factors if needed
00700   //
00701   
00702   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00703     *out << "\nCompute scale factors if needed ...\n";
00704   timer.start(true);
00705 
00706   if (firstTimeStateFuncScaling) {
00707     postEvalScalingSetup(epetraUnscaledOutArgs,out,verbLevel);
00708   }
00709 
00710   timer.stop();
00711   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00712     OSTab(out).o()
00713       << "\nTime to compute scale factors = "
00714       << timer.totalElapsedTime() <<" sec\n";
00715 
00716   //
00717   // D.2) Scale the output Epetra objects
00718   //
00719 
00720   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00721     *out << "\nScale the output objects ...\n";
00722   timer.start(true);
00723 
00724   EME::OutArgs epetraOutArgs = epetraModel_->createOutArgs();
00725   bool allFuncsWhereScaled = false;
00726   EpetraExt::scaleModelFuncs(
00727     epetraUnscaledOutArgs, epetraInArgsScaling_, epetraOutArgsScaling_,
00728     &epetraOutArgs, &allFuncsWhereScaled,
00729     out.get(), verbLevel
00730     );
00731   TEST_FOR_EXCEPTION(
00732     !allFuncsWhereScaled, std::logic_error,
00733     "Error, we can not currently handle epetra output objects that could not be"
00734     " scaled.  Special code will have to be added to handle this (i.e. using"
00735     " implicit diagonal and multiplied linear operators to implicitly do"
00736     " the scaling."
00737     );
00738 
00739   timer.stop();
00740   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00741     OSTab(out).o()
00742       << "\nTime to scale the output objects = "
00743       << timer.totalElapsedTime() << " sec\n";
00744 
00745   //
00746   // D.3) Convert any Epetra objects to Thyra OutArgs objects that still need to
00747   // be converted
00748   //
00749 
00750   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00751     *out << "\nFinish processing and wrapping the output objects ...\n";
00752   timer.start(true);
00753 
00754   finishConvertingOutArgsFromEpetraToThyra(
00755     epetraOutArgs, W, W_op, fwdW, efwdW, eW,
00756     outArgs
00757     );
00758 
00759   timer.stop();
00760   if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00761     OSTab(out).o()
00762       << "\nTime to finish processing and wrapping the output objects = "
00763       << timer.totalElapsedTime() <<" sec\n";
00764 
00765   //
00766   // E) Print footer to end the function
00767   //
00768 
00769   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00770   
00771 }
00772 
00773 
00774 // private
00775 
00776 
00777 void EpetraModelEvaluator::convertInArgsFromEpetraToThyra(
00778   const EpetraExt::ModelEvaluator::InArgs &epetraInArgs,
00779   ModelEvaluatorBase::InArgs<double> *inArgs
00780   ) const
00781 {
00782   
00783   using Teuchos::implicit_cast;
00784   typedef ModelEvaluatorBase MEB;
00785 
00786   TEST_FOR_EXCEPT(!inArgs);
00787 
00788   if(inArgs->supports(MEB::IN_ARG_x)) {
00789     inArgs->set_x( create_Vector( epetraInArgs.get_x(), x_space_ ) );
00790   }
00791   
00792   if(inArgs->supports(MEB::IN_ARG_x_dot)) {
00793     inArgs->set_x_dot( create_Vector( epetraInArgs.get_x_dot(), x_space_ ) );
00794   }
00795 
00796   const int Np = inArgs->Np();
00797   for( int l = 0; l < Np; ++l ) {
00798     inArgs->set_p( l, create_Vector( epetraInArgs.get_p(l), p_space_[l] ) );
00799   }
00800   
00801   if(inArgs->supports(MEB::IN_ARG_t)) {
00802     inArgs->set_t(epetraInArgs.get_t());
00803   }
00804   
00805 }
00806 
00807 
00808 void EpetraModelEvaluator::convertInArgsFromThyraToEpetra(
00809   const ModelEvaluatorBase::InArgs<double> &inArgs,
00810   EpetraExt::ModelEvaluator::InArgs *epetraInArgs
00811   ) const
00812 {
00813 
00814   using Teuchos::rcp;
00815   using Teuchos::rcp_const_cast;
00816   using Teuchos::Polynomial;
00817 
00818   TEST_FOR_EXCEPT(0==epetraInArgs);
00819 
00820   RCP<const VectorBase<double> > x_dot;
00821   if( inArgs.supports(IN_ARG_x_dot) && (x_dot = inArgs.get_x_dot()).get() ) {
00822     RCP<const Epetra_Vector> e_x_dot = get_Epetra_Vector(*x_map_,x_dot);
00823     epetraInArgs->set_x_dot(e_x_dot);
00824   }
00825 
00826   RCP<const VectorBase<double> > x;
00827   if( inArgs.supports(IN_ARG_x) && (x = inArgs.get_x()).get() ) {
00828     RCP<const Epetra_Vector> e_x = get_Epetra_Vector(*x_map_,x);
00829     epetraInArgs->set_x(e_x);
00830   }
00831 
00832   RCP<const VectorBase<double> > p_l;
00833   for(int l = 0;  l < inArgs.Np(); ++l ) {
00834     p_l = inArgs.get_p(l);
00835     if(p_l.get()) epetraInArgs->set_p(l,get_Epetra_Vector(*p_map_[l],p_l));
00836   }
00837 
00838   RCP<const Polynomial< VectorBase<double> > > x_dot_poly;
00839   RCP<Epetra_Vector> epetra_ptr;
00840   if(
00841     inArgs.supports(IN_ARG_x_dot_poly)
00842     && (x_dot_poly = inArgs.get_x_dot_poly()).get()
00843     )
00844   {
00845     RCP<Polynomial<Epetra_Vector> > epetra_x_dot_poly = 
00846       rcp(new Polynomial<Epetra_Vector>(x_dot_poly->degree()));
00847     for (unsigned int i=0; i<=x_dot_poly->degree(); i++) {
00848       epetra_ptr = rcp_const_cast<Epetra_Vector>(
00849         get_Epetra_Vector(*x_map_, x_dot_poly->getCoefficient(i)) );
00850       epetra_x_dot_poly->setCoefficientPtr(i,epetra_ptr);
00851     }
00852     epetraInArgs->set_x_dot_poly(epetra_x_dot_poly);
00853   }
00854   
00855   RCP<const Polynomial< VectorBase<double> > > x_poly;
00856   if(
00857     inArgs.supports(IN_ARG_x_poly)
00858     && (x_poly = inArgs.get_x_poly()).get()
00859     )
00860   {
00861     RCP<Polynomial<Epetra_Vector> > epetra_x_poly = 
00862       rcp(new Polynomial<Epetra_Vector>(x_poly->degree()));
00863     for (unsigned int i=0; i<=x_poly->degree(); i++) {
00864       epetra_ptr = rcp_const_cast<Epetra_Vector>(
00865         get_Epetra_Vector(*x_map_, x_poly->getCoefficient(i)) );
00866       epetra_x_poly->setCoefficientPtr(i,epetra_ptr);
00867     }
00868     epetraInArgs->set_x_poly(epetra_x_poly);
00869   }
00870 
00871   if( inArgs.supports(IN_ARG_t) )
00872     epetraInArgs->set_t(inArgs.get_t());
00873   
00874   if( inArgs.supports(IN_ARG_alpha) )
00875     epetraInArgs->set_alpha(inArgs.get_alpha());
00876   
00877   if( inArgs.supports(IN_ARG_beta) )
00878     epetraInArgs->set_beta(inArgs.get_beta());
00879 
00880 }
00881 
00882 
00883 void EpetraModelEvaluator::convertOutArgsFromThyraToEpetra(
00884   const ModelEvaluatorBase::OutArgs<double> &outArgs,
00885   EpetraExt::ModelEvaluator::OutArgs *epetraUnscaledOutArgs_inout,
00886   RCP<LinearOpWithSolveBase<double> > *W_out,
00887   RCP<LinearOpBase<double> > *W_op_out,
00888   RCP<const LinearOpBase<double> > *fwdW_out,
00889   RCP<EpetraLinearOp> *efwdW_out,
00890   RCP<Epetra_Operator> *eW_out
00891   ) const
00892 {
00893 
00894   using Teuchos::rcp;
00895   using Teuchos::rcp_const_cast;
00896   using Teuchos::rcp_dynamic_cast;
00897   using Teuchos::OSTab;
00898   using Teuchos::implicit_cast;
00899   using Thyra::get_Epetra_Vector;
00900   typedef EpetraExt::ModelEvaluator EME;
00901 
00902   // Assert input
00903 #ifdef TEUCHOS_DEBUG
00904   TEUCHOS_ASSERT(epetraUnscaledOutArgs_inout);
00905   TEUCHOS_ASSERT(W_out);
00906   TEUCHOS_ASSERT(W_op_out);
00907   TEUCHOS_ASSERT(fwdW_out);
00908   TEUCHOS_ASSERT(efwdW_out);
00909   TEUCHOS_ASSERT(eW_out);
00910 #endif
00911 
00912   // Create easy to use references
00913   EpetraExt::ModelEvaluator::OutArgs &epetraUnscaledOutArgs = *epetraUnscaledOutArgs_inout;
00914   RCP<LinearOpWithSolveBase<double> > &W = *W_out;
00915   RCP<LinearOpBase<double> > &W_op = *W_op_out;
00916   RCP<const LinearOpBase<double> > &fwdW = *fwdW_out;
00917   RCP<EpetraLinearOp> &efwdW = *efwdW_out;
00918   RCP<Epetra_Operator> &eW = *eW_out;
00919 
00920   // f
00921   { 
00922     RCP<VectorBase<double> > f;
00923     if( outArgs.supports(OUT_ARG_f) && (f = outArgs.get_f()).get() )
00924       epetraUnscaledOutArgs.set_f(get_Epetra_Vector(*f_map_,f));
00925   }
00926     
00927   // g
00928   {
00929     RCP<VectorBase<double> > g_j;
00930     for(int j = 0;  j < outArgs.Ng(); ++j ) {
00931       g_j = outArgs.get_g(j);
00932       if(g_j.get()) epetraUnscaledOutArgs.set_g(j,get_Epetra_Vector(*g_map_[j],g_j));
00933     }
00934   }
00935   
00936   // W and W_op
00937   {
00938 
00939     if( outArgs.supports(OUT_ARG_W) && (W = outArgs.get_W()).get() ) {
00940       Thyra::uninitializeOp<double>(*W_factory_,&*W,&fwdW);
00941       if(fwdW.get()) {
00942         efwdW = rcp_const_cast<EpetraLinearOp>(
00943           rcp_dynamic_cast<const EpetraLinearOp>(fwdW,true));
00944       }
00945       else {
00946         efwdW = this->create_epetra_W_op();
00947         fwdW = efwdW;
00948       }
00949     }
00950     if( outArgs.supports(OUT_ARG_W_op) && (W_op = outArgs.get_W_op()).get() ) {
00951       if( W_op.get() && !efwdW.get() )
00952         efwdW = rcp_const_cast<EpetraLinearOp>(
00953           rcp_dynamic_cast<const EpetraLinearOp>(W_op,true));
00954     }
00955     if(efwdW.get()) {
00956       // By the time we get here, if we have an object in efwdW, then it
00957       // should already be embeadded with an underlying Epetra_Operator object
00958       // that was allocated by the EpetraExt::ModelEvaluator object.
00959       // Therefore, we should just have to grab this object and be on our way.
00960       eW = efwdW->epetra_op();
00961       epetraUnscaledOutArgs.set_W(eW);
00962     }
00963     // NOTE: Above, if both W and W_op are set and have been through at least
00964     // one prior evaluation (and therefore have Epetra_Operator objects embedded
00965     // in them), then we will use the Epetra_Operator embedded in W to pass to
00966     // the EpetraExt::ModelEvaluator object and ignore the Epetra_Operator
00967     // object in W_op.  In the standard use case, these will be the same
00968     // Epetra_Operator objects.  However, it is possible that the client could
00969     // use this interface in such a way that these would have different
00970     // Epetra_Operator objects embedded in them.  In this (very unlikely) case,
00971     // the Epetra_Operator embedded in W_op will be discarded!  This might be
00972     // surprising to a client but it is very unlikely that this will ever be a
00973     // problem, but the issue is duly noted here!  Only dangerous programming
00974     // use of this interface would cause any problem.
00975     
00976     // Note: The following derivative objects update in place!
00977 
00978   }
00979 
00980   // DfDp
00981   {
00982     Derivative<double> DfDp_l;
00983     for(int l = 0;  l < outArgs.Np(); ++l ) {
00984       if( !outArgs.supports(OUT_ARG_DfDp,l).none()
00985         && !(DfDp_l = outArgs.get_DfDp(l)).isEmpty() )
00986       {
00987         epetraUnscaledOutArgs.set_DfDp(l,convert(DfDp_l,f_map_,p_map_[l]));
00988       }
00989     }
00990   }
00991 
00992   // DgDx_dot
00993   {
00994     Derivative<double> DgDx_dot_j;
00995     for(int j = 0;  j < outArgs.Ng(); ++j ) {
00996       if( !outArgs.supports(OUT_ARG_DgDx_dot,j).none()
00997         && !(DgDx_dot_j = outArgs.get_DgDx_dot(j)).isEmpty() )
00998       {
00999         epetraUnscaledOutArgs.set_DgDx_dot(j,convert(DgDx_dot_j,g_map_[j],x_map_));
01000       }
01001     }
01002   }
01003 
01004   // DgDx
01005   {
01006     Derivative<double> DgDx_j;
01007     for(int j = 0;  j < outArgs.Ng(); ++j ) {
01008       if( !outArgs.supports(OUT_ARG_DgDx,j).none()
01009         && !(DgDx_j = outArgs.get_DgDx(j)).isEmpty() )
01010       {
01011         epetraUnscaledOutArgs.set_DgDx(j,convert(DgDx_j,g_map_[j],x_map_));
01012       }
01013     }
01014   }
01015 
01016   // DgDp
01017   {
01018     DerivativeSupport DgDp_j_l_support;
01019     Derivative<double> DgDp_j_l;
01020     for (int j = 0;  j < outArgs.Ng(); ++j ) {
01021       for (int l = 0;  l < outArgs.Np(); ++l ) {
01022         if (!(DgDp_j_l_support = outArgs.supports(OUT_ARG_DgDp,j,l)).none()
01023           && !(DgDp_j_l = outArgs.get_DgDp(j,l)).isEmpty() )
01024         {
01025           epetraUnscaledOutArgs.set_DgDp(j,l,convert(DgDp_j_l,g_map_[j],p_map_[l]));
01026         }
01027       }
01028     }
01029   }
01030 
01031   // f_poly
01032   RCP<const Teuchos::Polynomial< VectorBase<double> > > f_poly;
01033   if( outArgs.supports(OUT_ARG_f_poly) && (f_poly = outArgs.get_f_poly()).get() )
01034   {
01035     RCP<Teuchos::Polynomial<Epetra_Vector> > epetra_f_poly = 
01036       Teuchos::rcp(new Teuchos::Polynomial<Epetra_Vector>(f_poly->degree()));
01037     for (unsigned int i=0; i<=f_poly->degree(); i++) {
01038       RCP<Epetra_Vector> epetra_ptr
01039         = Teuchos::rcp_const_cast<Epetra_Vector>(get_Epetra_Vector(*f_map_,
01040             f_poly->getCoefficient(i)));
01041       epetra_f_poly->setCoefficientPtr(i,epetra_ptr);
01042     }
01043     epetraUnscaledOutArgs.set_f_poly(epetra_f_poly);
01044   }
01045 
01046 }
01047 
01048 
01049 void EpetraModelEvaluator::preEvalScalingSetup(
01050   EpetraExt::ModelEvaluator::InArgs *epetraInArgs_inout,
01051   EpetraExt::ModelEvaluator::OutArgs *epetraUnscaledOutArgs_inout,
01052   const RCP<Teuchos::FancyOStream> &out,
01053   const Teuchos::EVerbosityLevel verbLevel
01054   ) const
01055 {
01056   
01057   typedef EpetraExt::ModelEvaluator EME;
01058   
01059 #ifdef TEUCHOS_DEBUG
01060   TEUCHOS_ASSERT(epetraInArgs_inout);
01061   TEUCHOS_ASSERT(epetraUnscaledOutArgs_inout);
01062 #endif
01063 
01064   EpetraExt::ModelEvaluator::InArgs
01065     &epetraInArgs = *epetraInArgs_inout;
01066   EpetraExt::ModelEvaluator::OutArgs
01067     &epetraUnscaledOutArgs = *epetraUnscaledOutArgs_inout;
01068 
01069   if (
01070     ( stateFunctionScaling_ == STATE_FUNC_SCALING_ROW_SUM )
01071     &&
01072     (
01073       epetraUnscaledOutArgs.supports(EME::OUT_ARG_f) 
01074       &&
01075       epetraUnscaledOutArgs.funcOrDerivesAreSet(EME::OUT_ARG_f)
01076       )
01077     &&
01078     (
01079       epetraUnscaledOutArgs.supports(EME::OUT_ARG_W)
01080       &&
01081       is_null(epetraUnscaledOutArgs.get_W())
01082       )
01083     )
01084   {
01085     // This is the first pass through with scaling turned on and the client
01086     // turned on automatic scaling but did not ask for W.  We must compute W
01087     // in order to compute the scale factors so we must allocate a temporary W
01088     // just to compute the scale factors and then throw it away.  If the
01089     // client wants to evaluate W at the same point, then it should have
01090     // passed W in but that is not our problem here.  The ModelEvaluator
01091     // relies on the client to set up the calls to allow for efficient
01092     // evaluation.
01093 
01094     if(out.get() && verbLevel >= Teuchos::VERB_LOW)
01095       *out
01096         << "\nCreating a temporary Epetra W to compute scale factors"
01097         << " for f(...) ...\n";
01098     epetraUnscaledOutArgs.set_W(create_and_assert_W(*epetraModel_));
01099     if( epetraInArgs.supports(EME::IN_ARG_beta) )
01100       epetraInArgs.set_beta(1.0);
01101     if( epetraInArgs.supports(EME::IN_ARG_alpha) )
01102       epetraInArgs.set_alpha(0.0);
01103   }
01104   
01105 }
01106 
01107 
01108 void EpetraModelEvaluator::postEvalScalingSetup(
01109   const EpetraExt::ModelEvaluator::OutArgs &epetraUnscaledOutArgs,
01110   const RCP<Teuchos::FancyOStream> &out,
01111   const Teuchos::EVerbosityLevel verbLevel
01112   ) const
01113 {
01114 
01115   using Teuchos::OSTab;
01116   using Teuchos::rcp;
01117   using Teuchos::rcp_const_cast;
01118   using Teuchos::includesVerbLevel;
01119 
01120   // Compute the scale factors for the state function f(...)
01121   switch(stateFunctionScaling_) {
01122 
01123     case STATE_FUNC_SCALING_ROW_SUM: {
01124 
01125       // Compute the inverse row-sum scaling from W
01126 
01127       const RCP<Epetra_RowMatrix>
01128         ermW = get_Epetra_RowMatrix(epetraUnscaledOutArgs);
01129       // Note: Above, we get the Epetra W object directly from the Epetra
01130       // OutArgs object since this might be a temporary matrix just to
01131       // compute scaling factors.  In this case, the stack funtion variable
01132       // eW might be empty!
01133 
01134       RCP<Epetra_Vector>
01135         invRowSums = rcp(new Epetra_Vector(ermW->OperatorRangeMap()));
01136       // Above: From the documentation is seems that the RangeMap should be
01137       // okay but who knows for sure!
01138 
01139       ermW->InvRowSums(*invRowSums);
01140 
01141       if (out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW)) {
01142         *out
01143           << "\nComputed inverse row sum scaling from W that"
01144           " will be used to scale f(...) and its derivatives:\n";
01145         double minVal = 0, maxVal = 0, avgVal = 0;
01146         invRowSums->MinValue(&minVal);
01147         invRowSums->MaxValue(&maxVal);
01148         invRowSums->MeanValue(&avgVal);
01149         OSTab tab(out);
01150         *out
01151           << "min(invRowSums) = " << minVal << "\n"
01152           << "max(invRowSums) = " << maxVal << "\n"
01153           << "avg(invRowSums) = " << avgVal << "\n";
01154       }
01155 
01156       stateFunctionScalingVec_ = invRowSums;
01157 
01158       break;
01159 
01160     }
01161 
01162     default:
01163       TEST_FOR_EXCEPT("Should never get here!");
01164 
01165   }
01166 
01167   epetraOutArgsScaling_ = epetraModel_->createOutArgs();
01168 
01169   epetraOutArgsScaling_.set_f(
01170     rcp_const_cast<Epetra_Vector>(stateFunctionScalingVec_) );
01171 
01172 }
01173 
01174 
01175 void EpetraModelEvaluator::finishConvertingOutArgsFromEpetraToThyra(
01176   const EpetraExt::ModelEvaluator::OutArgs &epetraOutArgs,
01177   RCP<LinearOpWithSolveBase<double> > &W,
01178   RCP<LinearOpBase<double> > &W_op,
01179   RCP<const LinearOpBase<double> > &fwdW,
01180   RCP<EpetraLinearOp> &efwdW,
01181   RCP<Epetra_Operator> &eW,
01182   const ModelEvaluatorBase::OutArgs<double> &outArgs
01183   ) const
01184 {
01185 
01186   using Teuchos::rcp_dynamic_cast;
01187   typedef EpetraExt::ModelEvaluator EME;
01188 
01189   if(efwdW.get()) {
01190     efwdW->setFullyInitialized(true); 
01191     // NOTE: Above will directly update W_op also if W.get()==NULL!
01192   }
01193   
01194   if( W.get() ) {
01195     Thyra::initializeOp<double>(*W_factory_,fwdW,&*W);
01196     W->setOStream(this->getOStream());
01197   }
01198 
01199   if( W_op.get() ) {
01200     if( W_op.shares_resource(efwdW) ) {
01201       // W_op was already updated above since *efwdW is the same object as *W_op
01202     }
01203     else {
01204       rcp_dynamic_cast<EpetraLinearOp>(W_op,true)->setFullyInitialized(true);
01205     }
01206   }
01207 
01208 }
01209 
01210 
01211 void EpetraModelEvaluator::updateNominalValuesAndBounds() const
01212 {
01213 
01214   using Teuchos::rcp;
01215   using Teuchos::implicit_cast;
01216   typedef ModelEvaluatorBase MEB;
01217   typedef EpetraExt::ModelEvaluator EME;
01218 
01219   if( !nominalValuesAndBoundsAreUpdated_ ) {
01220 
01221     // Gather the nominal values and bounds into Epetra InArgs objects
01222 
01223     EME::InArgs epetraOrigNominalValues;
01224     EpetraExt::gatherModelNominalValues(
01225       *epetraModel_, &epetraOrigNominalValues );
01226 
01227     EME::InArgs epetraOrigLowerBounds;
01228     EME::InArgs epetraOrigUpperBounds;
01229     EpetraExt::gatherModelBounds(
01230       *epetraModel_, &epetraOrigLowerBounds, &epetraOrigUpperBounds );
01231 
01232     // Set up Epetra InArgs scaling object
01233 
01234     epetraInArgsScaling_ = epetraModel_->createInArgs();
01235 
01236     if( !is_null(stateVariableScalingVec_) ) {
01237       invStateVariableScalingVec_
01238         = EpetraExt::createInverseModelScalingVector(stateVariableScalingVec_);
01239       if( epetraOrigNominalValues.supports(EME::IN_ARG_x_dot) ) {
01240         epetraInArgsScaling_.set_x_dot(invStateVariableScalingVec_);
01241       }
01242       if( epetraOrigNominalValues.supports(EME::IN_ARG_x) ) {
01243         epetraInArgsScaling_.set_x(invStateVariableScalingVec_);
01244       }
01245     }
01246     
01247     // Scale the original variables and bounds
01248 
01249     EME::InArgs epetraScaledNominalValues = epetraModel_->createInArgs();
01250     EpetraExt::scaleModelVars(
01251       epetraOrigNominalValues, epetraInArgsScaling_, &epetraScaledNominalValues
01252       );
01253 
01254     EME::InArgs epetraScaledLowerBounds = epetraModel_->createInArgs();
01255     EME::InArgs epetraScaledUpperBounds = epetraModel_->createInArgs();
01256     EpetraExt::scaleModelBounds(
01257       epetraOrigLowerBounds, epetraOrigUpperBounds, epetraModel_->getInfBound(),
01258       epetraInArgsScaling_,
01259       &epetraScaledLowerBounds, &epetraScaledUpperBounds
01260       );
01261 
01262     // Wrap the scaled epetra InArgs objects as Thyra InArgs objects!
01263 
01264     nominalValues_ = this->createInArgs();
01265     lowerBounds_ = this->createInArgs();
01266     upperBounds_ = this->createInArgs();
01267     convertInArgsFromEpetraToThyra(epetraScaledNominalValues, &nominalValues_);
01268     convertInArgsFromEpetraToThyra(epetraScaledLowerBounds, &lowerBounds_);
01269     convertInArgsFromEpetraToThyra(epetraScaledUpperBounds, &upperBounds_);
01270 
01271     nominalValuesAndBoundsAreUpdated_ = true;
01272 
01273   }
01274   else {
01275 
01276     // The nominal values and bounds should already be updated an should have
01277     // the currect scaling!
01278 
01279   }
01280 
01281 }
01282 
01283 
01284 void EpetraModelEvaluator::updateInArgsOutArgs() const
01285 {
01286 
01287   typedef EpetraExt::ModelEvaluator EME;
01288 
01289   const EpetraExt::ModelEvaluator &epetraModel = *epetraModel_;
01290   EME::InArgs  epetraInArgs  = epetraModel.createInArgs();
01291   EME::OutArgs epetraOutArgs = epetraModel.createOutArgs();
01292   const int Np = epetraOutArgs.Np();
01293   const int Ng = epetraOutArgs.Ng();
01294 
01295   //
01296   // InArgs
01297   //
01298 
01299   InArgsSetup<double> inArgs;
01300   inArgs.setModelEvalDescription(this->description());
01301   inArgs.set_Np(epetraInArgs.Np());
01302   inArgs.setSupports(IN_ARG_x_dot, epetraInArgs.supports(EME::IN_ARG_x_dot));
01303   inArgs.setSupports(IN_ARG_x, epetraInArgs.supports(EME::IN_ARG_x));
01304   inArgs.setSupports(IN_ARG_x_dot_poly,
01305     epetraInArgs.supports(EME::IN_ARG_x_dot_poly));
01306   inArgs.setSupports(IN_ARG_x_poly, epetraInArgs.supports(EME::IN_ARG_x_poly));
01307   inArgs.setSupports(IN_ARG_t, epetraInArgs.supports(EME::IN_ARG_t));
01308   inArgs.setSupports(IN_ARG_alpha, epetraInArgs.supports(EME::IN_ARG_alpha));
01309   inArgs.setSupports(IN_ARG_beta, epetraInArgs.supports(EME::IN_ARG_beta));
01310   prototypeInArgs_ = inArgs;
01311 
01312   //
01313   // OutArgs
01314   //
01315 
01316   OutArgsSetup<double> outArgs;
01317   outArgs.setModelEvalDescription(this->description());
01318   outArgs.set_Np_Ng(Np, Ng);
01319   // f
01320   outArgs.setSupports(OUT_ARG_f, epetraOutArgs.supports(EME::OUT_ARG_f));
01321   if (outArgs.supports(OUT_ARG_f)) {
01322     // W
01323     outArgs.setSupports(
01324       OUT_ARG_W, epetraOutArgs.supports(EME::OUT_ARG_W)&&!is_null(W_factory_));
01325     outArgs.setSupports(OUT_ARG_W_op,  epetraOutArgs.supports(EME::OUT_ARG_W));
01326     outArgs.set_W_properties(convert(epetraOutArgs.get_W_properties()));
01327     // DfDp
01328     for(int l=0; l<Np; ++l) {
01329       outArgs.setSupports(OUT_ARG_DfDp, l,
01330         convert(epetraOutArgs.supports(EME::OUT_ARG_DfDp, l)));
01331       if(!outArgs.supports(OUT_ARG_DfDp, l).none())
01332         outArgs.set_DfDp_properties(l,
01333           convert(epetraOutArgs.get_DfDp_properties(l)));
01334     }
01335   }
01336   // DgDx_dot and DgDx
01337   for(int j=0; j<Ng; ++j) {
01338     if (inArgs.supports(IN_ARG_x_dot))
01339       outArgs.setSupports(OUT_ARG_DgDx_dot, j,
01340         convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx_dot, j)));
01341     if(!outArgs.supports(OUT_ARG_DgDx_dot, j).none())
01342       outArgs.set_DgDx_dot_properties(j,
01343         convert(epetraOutArgs.get_DgDx_dot_properties(j)));
01344     if (inArgs.supports(IN_ARG_x))
01345       outArgs.setSupports(OUT_ARG_DgDx, j,
01346         convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx, j)));
01347     if(!outArgs.supports(OUT_ARG_DgDx, j).none())
01348       outArgs.set_DgDx_properties(j,
01349         convert(epetraOutArgs.get_DgDx_properties(j)));
01350   }
01351   // DgDp
01352   for(int j=0; j<Ng; ++j) for(int l=0; l<Np; ++l) {
01353     const EME::DerivativeSupport epetra_DgDp_j_l_support =
01354       epetraOutArgs.supports(EME::OUT_ARG_DgDp, j, l);
01355     outArgs.setSupports(OUT_ARG_DgDp, j, l,
01356       convert(epetra_DgDp_j_l_support));
01357     if(!outArgs.supports(OUT_ARG_DgDp, j, l).none())
01358       outArgs.set_DgDp_properties(j, l,
01359         convert(epetraOutArgs.get_DgDp_properties(j, l)));
01360   }
01361   outArgs.setSupports(OUT_ARG_f_poly,
01362     epetraOutArgs.supports(EME::OUT_ARG_f_poly));
01363   prototypeOutArgs_ = outArgs;
01364 
01365   // We are current!
01366   currentInArgsOutArgs_ = true;
01367 
01368 }
01369 
01370 
01371 RCP<EpetraLinearOp>
01372 EpetraModelEvaluator::create_epetra_W_op() const
01373 {
01374   return Thyra::partialNonconstEpetraLinearOp(
01375     this->get_f_space(), this->get_x_space(),
01376     create_and_assert_W(*epetraModel_)
01377     );
01378 }
01379 
01380 
01381 } // namespace Thyra
01382 
01383 
01384 //
01385 // Non-member utility functions
01386 //
01387 
01388 
01389 Teuchos::RCP<Thyra::EpetraModelEvaluator>
01390 Thyra::epetraModelEvaluator(
01391   const RCP<const EpetraExt::ModelEvaluator> &epetraModel,
01392   const RCP<LinearOpWithSolveFactoryBase<double> > &W_factory
01393   )
01394 {
01395   return Teuchos::rcp(new EpetraModelEvaluator(epetraModel,W_factory));
01396 }
01397 
01398 
01399 Thyra::ModelEvaluatorBase::EDerivativeMultiVectorOrientation
01400 Thyra::convert(
01401   const EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation &mvOrientation
01402   )
01403 {
01404   switch(mvOrientation) {
01405     case EpetraExt::ModelEvaluator::DERIV_MV_BY_COL :
01406       return ModelEvaluatorBase::DERIV_MV_BY_COL;
01407     case EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW :
01408       return ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW;
01409     default:
01410       TEST_FOR_EXCEPT(true);
01411   }
01412   return ModelEvaluatorBase::DERIV_MV_BY_COL; // Should never be called!
01413 }
01414 
01415 
01416 EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation
01417 Thyra::convert(
01418   const ModelEvaluatorBase::EDerivativeMultiVectorOrientation &mvOrientation
01419   )
01420 {
01421   switch(mvOrientation) {
01422     case ModelEvaluatorBase::DERIV_MV_BY_COL :
01423       return EpetraExt::ModelEvaluator::DERIV_MV_BY_COL;
01424     case ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW :
01425       return EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW;
01426     default:
01427       TEST_FOR_EXCEPT(true);
01428   }
01429   return EpetraExt::ModelEvaluator::DERIV_MV_BY_COL; // Should never be called!
01430 }
01431 
01432 
01433 Thyra::ModelEvaluatorBase::DerivativeProperties
01434 Thyra::convert(
01435   const EpetraExt::ModelEvaluator::DerivativeProperties &derivativeProperties
01436   )
01437 {
01438   ModelEvaluatorBase::EDerivativeLinearity linearity;
01439   switch(derivativeProperties.linearity) {
01440     case EpetraExt::ModelEvaluator::DERIV_LINEARITY_UNKNOWN:
01441       linearity = ModelEvaluatorBase::DERIV_LINEARITY_UNKNOWN;
01442       break;
01443     case EpetraExt::ModelEvaluator::DERIV_LINEARITY_CONST:
01444       linearity = ModelEvaluatorBase::DERIV_LINEARITY_CONST;
01445       break;
01446     case  EpetraExt::ModelEvaluator::DERIV_LINEARITY_NONCONST:
01447       linearity = ModelEvaluatorBase::DERIV_LINEARITY_NONCONST;
01448       break;
01449     default:
01450       TEST_FOR_EXCEPT(true);
01451   }
01452   ModelEvaluatorBase::ERankStatus rank;
01453   switch(derivativeProperties.rank) {
01454     case EpetraExt::ModelEvaluator::DERIV_RANK_UNKNOWN:
01455       rank = ModelEvaluatorBase::DERIV_RANK_UNKNOWN;
01456       break;
01457     case EpetraExt::ModelEvaluator::DERIV_RANK_FULL:
01458       rank = ModelEvaluatorBase::DERIV_RANK_FULL;
01459       break;
01460     case EpetraExt::ModelEvaluator::DERIV_RANK_DEFICIENT:
01461       rank = ModelEvaluatorBase::DERIV_RANK_DEFICIENT;
01462       break;
01463     default:
01464       TEST_FOR_EXCEPT(true);
01465   }
01466   return ModelEvaluatorBase::DerivativeProperties(
01467     linearity,rank,derivativeProperties.supportsAdjoint);
01468 }
01469 
01470 
01471 Thyra::ModelEvaluatorBase::DerivativeSupport
01472 Thyra::convert(
01473   const EpetraExt::ModelEvaluator::DerivativeSupport &derivativeSupport
01474   )
01475 {
01476   ModelEvaluatorBase::DerivativeSupport ds;
01477   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_LINEAR_OP))
01478     ds.plus(ModelEvaluatorBase::DERIV_LINEAR_OP);
01479   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_MV_BY_COL))
01480     ds.plus(ModelEvaluatorBase::DERIV_MV_BY_COL);
01481   if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW))
01482     ds.plus(ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW);
01483   return ds;
01484 }
01485 
01486 
01487 EpetraExt::ModelEvaluator::Derivative
01488 Thyra::convert(
01489   const ModelEvaluatorBase::Derivative<double> &derivative,
01490   const RCP<const Epetra_Map> &fnc_map,
01491   const RCP<const Epetra_Map> &var_map
01492   )
01493 {
01494   typedef ModelEvaluatorBase MEB;
01495   if(derivative.getLinearOp().get()) {
01496     return EpetraExt::ModelEvaluator::Derivative(
01497       Teuchos::rcp_const_cast<Epetra_Operator>(
01498         Teuchos::dyn_cast<const EpetraLinearOp>(*derivative.getLinearOp()).epetra_op()
01499         )
01500       );
01501   }
01502   else if(derivative.getDerivativeMultiVector().getMultiVector().get()) {
01503     return EpetraExt::ModelEvaluator::Derivative(
01504       EpetraExt::ModelEvaluator::DerivativeMultiVector(
01505         get_Epetra_MultiVector(
01506           ( derivative.getDerivativeMultiVector().getOrientation() == MEB::DERIV_MV_BY_COL
01507             ? *fnc_map
01508             : *var_map
01509             )
01510           ,derivative.getDerivativeMultiVector().getMultiVector()
01511           )
01512         ,convert(derivative.getDerivativeMultiVector().getOrientation())
01513         )
01514       );
01515   }
01516   return EpetraExt::ModelEvaluator::Derivative();
01517 }

Generated on Wed May 12 21:42:59 2010 for EpetraExt/Thyra Adapters by  doxygen 1.4.7