Thyra_EpetraModelEvaluator.cpp

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

Generated on Tue Oct 20 12:47:19 2009 for EpetraExt/Thyra Adapters by doxygen 1.4.7