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