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