Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_DefaultMultiPeriodModelEvaluator.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // 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 #ifndef THYRA_DEFAULT_MULTI_PERIOD_MODEL_EVALUATOR_HPP
00030 #define THYRA_DEFAULT_MULTI_PERIOD_MODEL_EVALUATOR_HPP
00031 
00032 
00033 #include "Thyra_ModelEvaluatorDefaultBase.hpp"
00034 #include "Thyra_DefaultProductVectorSpace.hpp"
00035 #include "Thyra_DefaultBlockedTriangularLinearOpWithSolveFactory.hpp" // Default implementation
00036 #include "Thyra_DefaultBlockedLinearOp.hpp"
00037 #include "Thyra_ModelEvaluatorDelegatorBase.hpp"
00038 #include "Thyra_ProductVectorBase.hpp"
00039 #include "Teuchos_implicit_cast.hpp"
00040 #include "Teuchos_AbstractFactory.hpp" // Interface
00041 #include "Teuchos_AbstractFactoryStd.hpp" // Implementation
00042 
00043 
00044 namespace Thyra {
00045 
00046 
00117 template<class Scalar>
00118 class DefaultMultiPeriodModelEvaluator
00119   : virtual public ModelEvaluatorDefaultBase<Scalar>
00120 {
00121 public:
00122 
00125 
00127   DefaultMultiPeriodModelEvaluator();
00128 
00130   DefaultMultiPeriodModelEvaluator(
00131     const int N,
00132     const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
00133     const Array<int> &z_indexes,
00134     const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
00135     const int g_index,
00136     const Array<Scalar> g_weights,
00137     const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space = Teuchos::null,
00138     const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space = Teuchos::null,
00139     const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory = Teuchos::null
00140     );
00141 
00200   void initialize(
00201     const int N,
00202     const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
00203     const Array<int> &z_indexes,
00204     const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
00205     const int g_index,
00206     const Array<Scalar> g_weights,
00207     const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space = Teuchos::null,
00208     const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space = Teuchos::null,
00209     const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory = Teuchos::null
00210     );
00211 
00221   void reset_z(
00222     const Array<Array<RCP<const VectorBase<Scalar> > > > &z
00223     );
00224   
00226 
00229 
00231   int Np() const;
00233   int Ng() const;
00235   RCP<const VectorSpaceBase<Scalar> > get_x_space() const;
00237   RCP<const VectorSpaceBase<Scalar> > get_f_space() const;
00239   RCP<const VectorSpaceBase<Scalar> > get_p_space(int l) const;
00241   RCP<const Array<std::string> > get_p_names(int l) const;
00243   RCP<const VectorSpaceBase<Scalar> > get_g_space(int j) const;
00245   ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
00247   ModelEvaluatorBase::InArgs<Scalar> getLowerBounds() const;
00249   ModelEvaluatorBase::InArgs<Scalar> getUpperBounds() const;
00251   RCP<LinearOpBase<Scalar> > create_W_op() const;
00253   RCP<const LinearOpWithSolveFactoryBase<Scalar> > get_W_factory() const;
00255   ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;
00257   void reportFinalPoint(
00258     const ModelEvaluatorBase::InArgs<Scalar> &finalPoint,
00259     const bool wasSolved
00260     );
00261 
00263 
00264 private:
00265 
00266 
00269 
00271   RCP<LinearOpBase<Scalar> > create_DfDp_op_impl(int l) const;
00273   RCP<LinearOpBase<Scalar> > create_DgDx_dot_op_impl(int j) const;
00275   RCP<LinearOpBase<Scalar> > create_DgDx_op_impl(int j) const;
00277   RCP<LinearOpBase<Scalar> > create_DgDp_op_impl(int j, int l) const;
00279   ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
00281   void evalModelImpl(
00282     const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00283     const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00284     ) const;
00285 
00287 
00288 private:
00289 
00290   // //////////////////////////
00291   // Private types
00292 
00293   typedef Array<Scalar> g_weights_t;
00294   typedef Array<Array<RCP<const VectorBase<Scalar> > > > z_t;
00295 
00296   // /////////////////////////
00297   // Private data members
00298 
00299   RCP<ModelEvaluator<Scalar> > periodModel_;
00300   Array<RCP<ModelEvaluator<Scalar> > > periodModels_;
00301   Array<int> z_indexes_;
00302   Array<int> period_l_map_;
00303   z_t z_; // size == N
00304   int g_index_;
00305   g_weights_t g_weights_; // size == N
00306   RCP<const ProductVectorSpaceBase<Scalar> > x_bar_space_;
00307   RCP<const ProductVectorSpaceBase<Scalar> > f_bar_space_;
00308   RCP<LinearOpWithSolveFactoryBase<Scalar> > W_bar_factory_;
00309   int Np_;
00310   int Ng_;
00311   ModelEvaluatorBase::InArgs<Scalar> nominalValues_;
00312   ModelEvaluatorBase::InArgs<Scalar> lowerBounds_;
00313   ModelEvaluatorBase::InArgs<Scalar> upperBounds_;
00314 
00315   // /////////////////////////
00316   // Private member functions
00317 
00318   void set_z_indexes_and_create_period_l_map( const Array<int> &z_indexes );
00319 
00320   void wrapNominalValuesAndBounds();
00321 
00322   static
00323   RCP<ProductVectorBase<Scalar> >
00324   createProductVector(
00325     const RCP<const ProductVectorSpaceBase<Scalar> > &prodVecSpc
00326     );
00327 
00328   // Return the index of a "free" parameter in the period model given its
00329   // index in this mulit-period model.
00330   int period_l(const int l) const
00331     {
00332       TEST_FOR_EXCEPT( !( 0 <= l && l < Np_) );
00333       return period_l_map_[l];
00334     }
00335 
00336   int numPeriodZs() const { return z_indexes_.size(); }
00337 
00338   int N() const { return x_bar_space_->numBlocks(); }
00339   
00340 };
00341 
00342 
00343 // /////////////////////////////////
00344 // Implementations
00345 
00346 
00347 // Constructors/Intializers/Accessors
00348 
00349 
00350 template<class Scalar>
00351 DefaultMultiPeriodModelEvaluator<Scalar>::DefaultMultiPeriodModelEvaluator()
00352   :g_index_(-1), Np_(-1), Ng_(-1)
00353 {}
00354 
00355 
00356 template<class Scalar>
00357 DefaultMultiPeriodModelEvaluator<Scalar>::DefaultMultiPeriodModelEvaluator(
00358   const int N,
00359   const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
00360   const Array<int> &z_indexes,
00361   const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
00362   const int g_index,
00363   const Array<Scalar> g_weights,
00364   const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space,
00365   const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space,
00366   const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
00367   )
00368   :g_index_(-1), Np_(-1), Ng_(-1)
00369 {
00370   initialize(
00371     N, periodModels, z_indexes, z, g_index, g_weights,
00372     x_bar_space, f_bar_space, W_bar_factory
00373     );
00374 }
00375 
00376 
00377 template<class Scalar>
00378 void DefaultMultiPeriodModelEvaluator<Scalar>::initialize(
00379   const int N,
00380   const Array<RCP<ModelEvaluator<Scalar> > > &periodModels,
00381   const Array<int> &z_indexes,
00382   const Array<Array<RCP<const VectorBase<Scalar> > > > &z,
00383   const int g_index,
00384   const Array<Scalar> g_weights,
00385   const RCP<const ProductVectorSpaceBase<Scalar> > &x_bar_space,
00386   const RCP<const ProductVectorSpaceBase<Scalar> > &f_bar_space,
00387   const RCP<LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
00388   )
00389 {
00390 
00391   using Teuchos::implicit_cast;
00392   typedef Teuchos::ScalarTraits<Scalar> ST;
00393   typedef ModelEvaluatorBase MEB;
00394 
00395   TEST_FOR_EXCEPT( N <= 0 );
00396   TEST_FOR_EXCEPT( implicit_cast<int>(periodModels.size()) != N );
00397   MEB::InArgs<Scalar> periodInArgs = periodModels[0]->createInArgs();
00398   MEB::OutArgs<Scalar> periodOutArgs = periodModels[0]->createOutArgs();
00399   for ( int k = 0; k < implicit_cast<int>(z_indexes.size()); ++k ) {
00400     TEST_FOR_EXCEPT( !( 0 <= z_indexes[k] && z_indexes[k] < periodInArgs.Np() ) );
00401   }
00402   TEST_FOR_EXCEPT( implicit_cast<int>(z.size())!=N );
00403   TEST_FOR_EXCEPT( !( 0 <= g_index && g_index < periodOutArgs.Ng() ) );
00404   TEST_FOR_EXCEPT( implicit_cast<int>(g_weights.size())!=N );
00405   // ToDo: Assert that the different period models are compatible!
00406 
00407   Np_ = periodInArgs.Np() - z_indexes.size();
00408 
00409   Ng_ = 1;
00410 
00411   set_z_indexes_and_create_period_l_map(z_indexes);
00412   
00413   periodModel_ = periodModels[0]; // Assume basic structure!
00414   
00415   periodModels_ = periodModels;
00416   
00417   z_.resize(N);
00418   reset_z(z);
00419   
00420   g_index_ = g_index;
00421   g_weights_ = g_weights;
00422   
00423   if ( periodInArgs.supports(MEB::IN_ARG_x) ) {
00424     if( !is_null(x_bar_space) ) {
00425       TEST_FOR_EXCEPT(!(x_bar_space->numBlocks()==N));
00426       // ToDo: Check the constituent spaces more carefully against models[]->get_x_space().
00427       x_bar_space_ = x_bar_space;
00428     }
00429     else {
00430       x_bar_space_ = productVectorSpace(periodModel_->get_x_space(),N);
00431       // ToDo: Update to build from different models for different x spaces!
00432     }
00433   }
00434 
00435   if ( periodOutArgs.supports(MEB::OUT_ARG_f) ) {
00436     if( !is_null(f_bar_space) ) {
00437       TEST_FOR_EXCEPT(!(f_bar_space->numBlocks()==N));
00438       // ToDo: Check the constituent spaces more carefully against models[]->get_f_space().
00439       f_bar_space_ = f_bar_space;
00440     }
00441     else {
00442       f_bar_space_ = productVectorSpace(periodModel_->get_f_space(),N);
00443       // ToDo: Update to build from different models for different f spaces!
00444     }
00445   }
00446 
00447   if ( periodOutArgs.supports(MEB::OUT_ARG_W) ) {
00448     if ( !is_null(W_bar_factory) ) {
00449       // ToDo: Check the compatability of the W_bar objects created using this object!
00450       W_bar_factory_ = W_bar_factory;
00451     }
00452     else {
00453       W_bar_factory_ =
00454         defaultBlockedTriangularLinearOpWithSolveFactory<Scalar>(
00455           periodModel_->get_W_factory() );
00456     }
00457   }
00458 
00459   wrapNominalValuesAndBounds();
00460   
00461 }
00462 
00463 
00464 template<class Scalar>
00465 void DefaultMultiPeriodModelEvaluator<Scalar>::reset_z(
00466   const Array<Array<RCP<const VectorBase<Scalar> > > > &z
00467   )
00468 {
00469 
00470   using Teuchos::implicit_cast;
00471   
00472   const int N = z_.size();
00473 
00474 #ifdef TEUCHOS_DEBUG
00475   TEST_FOR_EXCEPT( N == 0 && "Error, must have called initialize() first!" );
00476   TEST_FOR_EXCEPT( implicit_cast<int>(z.size()) != N ); 
00477 #endif 
00478 
00479   for( int i = 0; i < N; ++i ) {
00480     const Array<RCP<const VectorBase<Scalar> > >  &z_i = z[i];
00481 #ifdef TEUCHOS_DEBUG
00482     TEST_FOR_EXCEPT( z_i.size() != z_indexes_.size() ); 
00483 #endif 
00484     z_[i] = z_i;
00485   }
00486   
00487 }
00488 
00489 
00490 // Public functions overridden from ModelEvaulator
00491 
00492 
00493 template<class Scalar>
00494 int DefaultMultiPeriodModelEvaluator<Scalar>::Np() const
00495 {
00496   return Np_;
00497 }
00498 
00499 
00500 template<class Scalar>
00501 int DefaultMultiPeriodModelEvaluator<Scalar>::Ng() const
00502 {
00503   return Ng_;
00504 }
00505 
00506 
00507 template<class Scalar>
00508 RCP<const VectorSpaceBase<Scalar> >
00509 DefaultMultiPeriodModelEvaluator<Scalar>::get_x_space() const
00510 {
00511   return x_bar_space_;
00512 }
00513 
00514 
00515 template<class Scalar>
00516 RCP<const VectorSpaceBase<Scalar> >
00517 DefaultMultiPeriodModelEvaluator<Scalar>::get_f_space() const
00518 {
00519   return f_bar_space_;
00520 }
00521 
00522 
00523 template<class Scalar>
00524 RCP<const VectorSpaceBase<Scalar> >
00525 DefaultMultiPeriodModelEvaluator<Scalar>::get_p_space(int l) const
00526 {
00527   return  periodModel_->get_p_space(period_l(l));
00528 }
00529 
00530 
00531 template<class Scalar>
00532 RCP<const Array<std::string> >
00533 DefaultMultiPeriodModelEvaluator<Scalar>::get_p_names(int l) const
00534 {
00535   return  periodModel_->get_p_names(period_l(l));
00536 }
00537 
00538 
00539 template<class Scalar>
00540 RCP<const VectorSpaceBase<Scalar> >
00541 DefaultMultiPeriodModelEvaluator<Scalar>::get_g_space(int j) const
00542 {
00543   TEST_FOR_EXCEPT(j!=0);
00544   return periodModel_->get_g_space(g_index_);
00545 }
00546 
00547 
00548 template<class Scalar>
00549 ModelEvaluatorBase::InArgs<Scalar>
00550 DefaultMultiPeriodModelEvaluator<Scalar>::getNominalValues() const
00551 {
00552   return nominalValues_;
00553 }
00554 
00555 
00556 template<class Scalar>
00557 ModelEvaluatorBase::InArgs<Scalar>
00558 DefaultMultiPeriodModelEvaluator<Scalar>::getLowerBounds() const
00559 {
00560   return lowerBounds_;
00561 }
00562 
00563 
00564 template<class Scalar>
00565 ModelEvaluatorBase::InArgs<Scalar>
00566 DefaultMultiPeriodModelEvaluator<Scalar>::getUpperBounds() const
00567 {
00568   return upperBounds_;
00569 }
00570 
00571 
00572 template<class Scalar>
00573 RCP<LinearOpBase<Scalar> >
00574 DefaultMultiPeriodModelEvaluator<Scalar>::create_W_op() const
00575 {
00576   // Set up the block structure ready to have the blocks filled!
00577   const RCP<PhysicallyBlockedLinearOpBase<Scalar> >
00578     W_op_bar = defaultBlockedLinearOp<Scalar>();
00579   W_op_bar->beginBlockFill(f_bar_space_,x_bar_space_);
00580   const int N = x_bar_space_->numBlocks();
00581   for ( int i = 0; i < N; ++i ) {
00582     W_op_bar->setNonconstBlock( i, i, periodModel_->create_W_op() );
00583   }
00584   W_op_bar->endBlockFill();
00585   return W_op_bar;
00586 }
00587 
00588 
00589 template<class Scalar>
00590 RCP<const LinearOpWithSolveFactoryBase<Scalar> >
00591 DefaultMultiPeriodModelEvaluator<Scalar>::get_W_factory() const
00592 {
00593   return W_bar_factory_;
00594 }
00595 
00596 
00597 template<class Scalar>
00598 ModelEvaluatorBase::InArgs<Scalar>
00599 DefaultMultiPeriodModelEvaluator<Scalar>::createInArgs() const
00600 {
00601   typedef ModelEvaluatorBase MEB;
00602   MEB::InArgs<Scalar> periodInArgs = periodModel_->createInArgs();
00603   MEB::InArgsSetup<Scalar> inArgs;
00604   inArgs.setModelEvalDescription(this->description());
00605   inArgs.set_Np(Np_);
00606   inArgs.setSupports( MEB::IN_ARG_x, periodInArgs.supports(MEB::IN_ARG_x) );
00607   return inArgs;
00608 }
00609 
00610 
00611 template<class Scalar>
00612 void DefaultMultiPeriodModelEvaluator<Scalar>::reportFinalPoint(
00613   const ModelEvaluatorBase::InArgs<Scalar>      &finalPoint
00614   ,const bool                                   wasSolved
00615   )
00616 {
00617   // We are just going to ignore the final point here.  It is not clear how to
00618   // report a "final" point back to the underlying *periodModel_ object since
00619   // we have so many different "points" that we could return (i.e. one for
00620   // each period).  I guess we could report back the final parameter values
00621   // (other than the z parameter) but there are multiple states x[i] and
00622   // period parameters z[i] that we can report back.
00623 }
00624 
00625 
00626 // Public functions overridden from ModelEvaulatorDefaultBase
00627 
00628 
00629 template<class Scalar>
00630 RCP<LinearOpBase<Scalar> >
00631 DefaultMultiPeriodModelEvaluator<Scalar>::create_DfDp_op_impl(int l) const
00632 {
00633   TEST_FOR_EXCEPT("This class does not support DfDp(l) as a linear operator yet.");
00634   return Teuchos::null;
00635 }
00636 
00637 
00638 template<class Scalar>
00639 RCP<LinearOpBase<Scalar> >
00640 DefaultMultiPeriodModelEvaluator<Scalar>::create_DgDx_dot_op_impl(int j) const
00641 {
00642   TEST_FOR_EXCEPT("This class does not support DgDx_dot(j) as a linear operator yet.");
00643   return Teuchos::null;
00644 }
00645 
00646 
00647 template<class Scalar>
00648 RCP<LinearOpBase<Scalar> >
00649 DefaultMultiPeriodModelEvaluator<Scalar>::create_DgDx_op_impl(int j) const
00650 {
00651   TEST_FOR_EXCEPT("This class does not support DgDx(j) as a linear operator yet.");
00652   return Teuchos::null;
00653 }
00654 
00655 
00656 template<class Scalar>
00657 RCP<LinearOpBase<Scalar> >
00658 DefaultMultiPeriodModelEvaluator<Scalar>::create_DgDp_op_impl(int j, int l) const
00659 {
00660   TEST_FOR_EXCEPT("This class does not support DgDp(j,l) as a linear operator yet.");
00661   return Teuchos::null;
00662 }
00663 
00664 
00665 template<class Scalar>
00666 ModelEvaluatorBase::OutArgs<Scalar>
00667 DefaultMultiPeriodModelEvaluator<Scalar>::createOutArgsImpl() const
00668 {
00669 
00670   typedef ModelEvaluatorBase MEB;
00671 
00672   MEB::OutArgs<Scalar> periodOutArgs = periodModel_->createOutArgs();
00673   MEB::OutArgsSetup<Scalar> outArgs;
00674 
00675   outArgs.setModelEvalDescription(this->description());
00676 
00677   outArgs.set_Np_Ng(Np_,Ng_);
00678 
00679   // f
00680   if (periodOutArgs.supports(MEB::OUT_ARG_f) ) {
00681     outArgs.setSupports(MEB::OUT_ARG_f);
00682   }
00683 
00684   // W_op
00685   if (periodOutArgs.supports(MEB::OUT_ARG_W_op) ) {
00686     outArgs.setSupports(MEB::OUT_ARG_W_op);
00687     outArgs.set_W_properties(periodOutArgs.get_W_properties());
00688   }
00689   // Note: We will not directly support the LOWSB form W as we will let the
00690   // default base class handle this given our W_factory!
00691 
00692   // DfDp(l)
00693   for ( int l = 0; l < Np_; ++l ) {
00694     const int period_l = this->period_l(l);
00695     const MEB::DerivativeSupport period_DfDp_l_support
00696       = periodOutArgs.supports(MEB::OUT_ARG_DfDp,period_l);
00697     if (!period_DfDp_l_support.none()) {
00698       outArgs.setSupports( MEB::OUT_ARG_DfDp, l, period_DfDp_l_support );
00699       outArgs.set_DfDp_properties(
00700         l, periodOutArgs.get_DfDp_properties(period_l) );
00701     }
00702   }
00703 
00704   // DgDx_dot
00705   const MEB::DerivativeSupport
00706     period_DgDx_dot_support = periodOutArgs.supports(MEB::OUT_ARG_DgDx_dot,g_index_);
00707   if (!period_DgDx_dot_support.none()) {
00708     outArgs.setSupports( MEB::OUT_ARG_DgDx_dot, 0, period_DgDx_dot_support );
00709     outArgs.set_DgDx_dot_properties(
00710       0, periodOutArgs.get_DgDx_dot_properties(g_index_) );
00711   }
00712 
00713   // DgDx
00714   const MEB::DerivativeSupport
00715     period_DgDx_support = periodOutArgs.supports(MEB::OUT_ARG_DgDx,g_index_);
00716   if (!period_DgDx_support.none()) {
00717     outArgs.setSupports( MEB::OUT_ARG_DgDx, 0, period_DgDx_support );
00718     outArgs.set_DgDx_properties(
00719       0, periodOutArgs.get_DgDx_properties(g_index_) );
00720   }
00721 
00722   // DgDp(l)
00723   for ( int l = 0; l < Np_; ++l ) {
00724     const int period_l = this->period_l(l);
00725     const MEB::DerivativeSupport period_DgDp_l_support
00726       = periodOutArgs.supports(MEB::OUT_ARG_DgDp, g_index_, period_l);
00727     if (!period_DgDp_l_support.none()) {
00728       outArgs.setSupports( MEB::OUT_ARG_DgDp, 0, l, period_DgDp_l_support );
00729       outArgs.set_DgDp_properties(
00730         0, l, periodOutArgs.get_DgDp_properties(g_index_,period_l) );
00731     }
00732   }
00733   
00734   return outArgs;
00735 
00736 }
00737 
00738 
00739 template<class Scalar>
00740 void DefaultMultiPeriodModelEvaluator<Scalar>::evalModelImpl(
00741   const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00742   const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00743   ) const
00744 {
00745 
00746   using Teuchos::rcp_dynamic_cast;
00747   typedef Teuchos::ScalarTraits<Scalar> ST;
00748   typedef ModelEvaluatorBase MEB;
00749   typedef Teuchos::VerboseObjectTempState<ModelEvaluatorBase> VOTSME;
00750 
00751   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
00752     "DefaultMultiPeriodModelEvaluator",inArgs,outArgs,periodModel_ );
00753   // ToDo: You will have to set the verbosity level for each of the
00754   // periodModels_[i] individually below!
00755   
00756   const int N = x_bar_space_->numBlocks();
00757   const int Np = this->Np_;
00758   //const int Ng = this->Ng_;
00759 
00760   //
00761   // A) Setup InArgs
00762   //
00763 
00764   RCP<const ProductVectorBase<Scalar> > x_bar;
00765   if (inArgs.supports(MEB::IN_ARG_x)) {
00766     x_bar = rcp_dynamic_cast<const ProductVectorBase<Scalar> >(
00767       inArgs.get_x(), true );
00768     TEST_FOR_EXCEPTION(
00769       is_null(x_bar), std::logic_error,
00770       "Error, if x is supported, it must be set!"
00771       );
00772   }
00773 
00774   //
00775   // B) Setup OutArgs
00776   //
00777   
00778   RCP<ProductVectorBase<Scalar> > f_bar;
00779   if (outArgs.supports(MEB::OUT_ARG_f)) {
00780     f_bar = rcp_dynamic_cast<ProductVectorBase<Scalar> >(
00781       outArgs.get_f(), true );
00782   }
00783 
00784   Array<MEB::Derivative<Scalar> > DfDp_bar(Np);
00785   Array<RCP<ProductMultiVectorBase<Scalar> > > DfDp_bar_mv(Np);
00786   for ( int l = 0; l < Np; ++l ) {
00787     if (!outArgs.supports(MEB::OUT_ARG_DfDp,l).none()) {
00788       MEB::Derivative<Scalar>
00789         DfDp_bar_l = outArgs.get_DfDp(l);
00790       DfDp_bar[l] = DfDp_bar_l;
00791       DfDp_bar_mv[l] = rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(
00792         DfDp_bar_l.getMultiVector(), true );
00793       TEST_FOR_EXCEPTION(
00794         (
00795           !DfDp_bar_l.isEmpty()
00796           &&
00797           (
00798             is_null(DfDp_bar_mv[l])
00799             ||
00800             DfDp_bar_l.getMultiVectorOrientation() != MEB::DERIV_MV_BY_COL
00801             )
00802           ),
00803         std::logic_error,
00804         "Error, we currently can only handle DfDp as an column-based multi-vector!"
00805         );
00806     }
00807   }
00808 
00809   RCP<BlockedLinearOpBase<Scalar> > W_op_bar;
00810   if (outArgs.supports(MEB::OUT_ARG_W_op)) {
00811     W_op_bar = rcp_dynamic_cast<BlockedLinearOpBase<Scalar> >(
00812       outArgs.get_W_op(), true
00813       );
00814   }
00815 
00816   RCP<VectorBase<Scalar> >
00817     g_bar = outArgs.get_g(0);
00818 
00819   MEB::Derivative<Scalar> DgDx_dot_bar;
00820   RCP<ProductMultiVectorBase<Scalar> > DgDx_dot_bar_mv;
00821   if (!outArgs.supports(MEB::OUT_ARG_DgDx_dot,0).none()) {
00822     DgDx_dot_bar = outArgs.get_DgDx_dot(0);
00823     DgDx_dot_bar_mv = rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(
00824       DgDx_dot_bar.getMultiVector(), true );
00825     TEST_FOR_EXCEPTION(
00826       (
00827         !DgDx_dot_bar.isEmpty()
00828         &&
00829         (
00830           is_null(DgDx_dot_bar_mv)
00831           ||
00832           DgDx_dot_bar.getMultiVectorOrientation() != MEB::DERIV_TRANS_MV_BY_ROW
00833           )
00834         ),
00835       std::logic_error,
00836       "Error, we currently can only handle DgDx_dot as an row-based multi-vector!"
00837       );
00838   }
00839 
00840   MEB::Derivative<Scalar> DgDx_bar;
00841   RCP<ProductMultiVectorBase<Scalar> > DgDx_bar_mv;
00842   if (!outArgs.supports(MEB::OUT_ARG_DgDx,0).none()) {
00843     DgDx_bar = outArgs.get_DgDx(0);
00844     DgDx_bar_mv = rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(
00845       DgDx_bar.getMultiVector(), true );
00846     TEST_FOR_EXCEPTION(
00847       (
00848         !DgDx_bar.isEmpty()
00849         &&
00850         (
00851           is_null(DgDx_bar_mv)
00852           ||
00853           DgDx_bar.getMultiVectorOrientation() != MEB::DERIV_TRANS_MV_BY_ROW
00854           )
00855         ),
00856       std::logic_error,
00857       "Error, we currently can only handle DgDx as an row-based multi-vector!"
00858       );
00859   }
00860 
00861   Array<MEB::Derivative<Scalar> > DgDp_bar(Np);
00862   Array<RCP<MultiVectorBase<Scalar> > > DgDp_bar_mv(Np);
00863   for ( int l = 0; l < Np; ++l ) {
00864     if (!outArgs.supports(MEB::OUT_ARG_DgDp,0,l).none()) {
00865       MEB::Derivative<Scalar>
00866         DgDp_bar_l = outArgs.get_DgDp(0,l);
00867       DgDp_bar[l] = DgDp_bar_l;
00868       DgDp_bar_mv[l] = DgDp_bar_l.getMultiVector();
00869       TEST_FOR_EXCEPTION(
00870         !DgDp_bar_l.isEmpty() && is_null(DgDp_bar_mv[l]),
00871         std::logic_error,
00872         "Error, we currently can only handle DgDp as some type of multi-vector!"
00873         );
00874     }
00875   }
00876 
00877   //
00878   // C) Evaluate the model
00879   //
00880 
00881   // C.1) Set up storage and do some initializations
00882 
00883   MEB::InArgs<Scalar>
00884     periodInArgs = periodModel_->createInArgs();
00885   // ToDo: The above will have to change if you allow different structures for
00886   // each period model!
00887   
00888   // Set all of the parameters that will just be passed through 
00889   for ( int l = 0; l < Np; ++l )
00890     periodInArgs.set_p( period_l(l), inArgs.get_p(l) ); // Can be null just fine
00891   
00892   MEB::OutArgs<Scalar>
00893     periodOutArgs = periodModel_->createOutArgs();
00894   // ToDo: The above will have to change if you allow different structures for
00895   // each period model!
00896   
00897   // Create storage for period g's that will be summed into global g_bar
00898   periodOutArgs.set_g(
00899     g_index_, createMember<Scalar>( periodModel_->get_g_space(g_index_) ) );
00900 
00901   // Zero out global g_bar that will be summed into below
00902   if (!is_null(g_bar) )
00903     assign( &*g_bar, ST::zero() );
00904   
00905   // Set up storage for peroid DgDp[l] objects that will be summed into global
00906   // DgDp_bar[l] and zero out DgDp_bar[l] that will be summed into.
00907   for ( int l = 0; l < Np; ++l ) {
00908     if ( !is_null(DgDp_bar_mv[l]) ) {
00909       assign(&*DgDp_bar_mv[l],ST::zero());
00910       periodOutArgs.set_DgDp(
00911         g_index_, period_l(l),
00912         create_DgDp_mv(
00913           *periodModel_, g_index_, period_l(l),
00914           DgDp_bar[l].getMultiVectorOrientation()
00915           )
00916         );
00917     }
00918   }
00919   
00920   // C.2) Loop over periods and assemble the model
00921 
00922   for ( int i = 0; i < N; ++i ) {
00923 
00924     VOTSME thyraModel_outputTempState(periodModels_[i],out,verbLevel);
00925 
00926     // C.2.a) Set period-speicific InArgs and OutArgs
00927 
00928     for ( int k = 0; k < numPeriodZs(); ++k )
00929       periodInArgs.set_p( z_indexes_[k], z_[i][k] );
00930 
00931     if (!is_null(x_bar))
00932       periodInArgs.set_x(x_bar->getVectorBlock(i));
00933 
00934     if (!is_null(f_bar))
00935       periodOutArgs.set_f(f_bar->getNonconstVectorBlock(i)); // Updated in place!
00936 
00937     if ( !is_null(W_op_bar) )
00938       periodOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,i));
00939 
00940     for ( int l = 0; l < Np; ++l ) {
00941       if ( !is_null(DfDp_bar_mv[l]) ) {
00942         periodOutArgs.set_DfDp(
00943           period_l(l),
00944           MEB::Derivative<Scalar>(
00945             DfDp_bar_mv[l]->getNonconstMultiVectorBlock(i),
00946             MEB::DERIV_MV_BY_COL
00947             )
00948           );
00949       }
00950     }
00951     
00952     if ( !is_null(DgDx_dot_bar_mv) ) {
00953       periodOutArgs.set_DgDx_dot(
00954         g_index_,
00955         MEB::Derivative<Scalar>(
00956           DgDx_dot_bar_mv->getNonconstMultiVectorBlock(i),
00957           MEB::DERIV_TRANS_MV_BY_ROW
00958           )
00959         );
00960     }
00961     
00962     if ( !is_null(DgDx_bar_mv) ) {
00963       periodOutArgs.set_DgDx(
00964         g_index_,
00965         MEB::Derivative<Scalar>(
00966           DgDx_bar_mv->getNonconstMultiVectorBlock(i),
00967           MEB::DERIV_TRANS_MV_BY_ROW
00968           )
00969         );
00970     }
00971 
00972     // C.2.b) Evaluate the period model
00973 
00974     periodModels_[i]->evalModel( periodInArgs, periodOutArgs );
00975 
00976     // C.2.c) Process output arguments that need processed
00977 
00978     // Sum into global g_bar
00979     if (!is_null(g_bar)) {
00980       Vp_StV( &*g_bar, g_weights_[i], *periodOutArgs.get_g(g_index_) );
00981     }
00982 
00983     // Sum into global DgDp_bar
00984     for ( int l = 0; l < Np; ++l ) {
00985       if ( !is_null(DgDp_bar_mv[l]) ) {
00986         update(
00987           g_weights_[i],
00988           *periodOutArgs.get_DgDp(g_index_,period_l(l)).getMultiVector(),
00989           &*DgDp_bar_mv[l]
00990           );
00991       }
00992     }
00993     
00994     // Scale DgDx_dot_bar_mv[i]
00995     if ( !is_null(DgDx_dot_bar_mv) ) {
00996       scale( g_weights_[i], &*DgDx_dot_bar_mv->getNonconstMultiVectorBlock(i) );
00997     }
00998     
00999     // Scale DgDx_bar_mv[i]
01000     if ( !is_null(DgDx_bar_mv) ) {
01001       scale( g_weights_[i], &*DgDx_bar_mv->getNonconstMultiVectorBlock(i) );
01002     }
01003 
01004   }
01005 
01006   // ToDo: We need to do some type of global sum of g_bar and DgDp_bar to
01007   // account for other clusters of processes.  I might do this with a separate
01008   // non-ANA class.
01009 
01010   // Once we get here, all of the quantities should be updated and we should
01011   // be all done!
01012 
01013   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
01014 
01015 }
01016 
01017 
01018 // private
01019 
01020 
01021 template<class Scalar>
01022 void DefaultMultiPeriodModelEvaluator<Scalar>::set_z_indexes_and_create_period_l_map(
01023   const Array<int> &z_indexes
01024   )
01025 {
01026 #ifdef TEUCHOS_DEBUG
01027   TEST_FOR_EXCEPT( Np_ <= 0 && "Error, Np must be set!" );
01028 #endif
01029   z_indexes_ = z_indexes;
01030   period_l_map_.resize(0);
01031   const int numTotalParams = Np_ + z_indexes_.size();
01032   Array<int>::const_iterator
01033     z_indexes_itr = z_indexes_.begin(),
01034     z_indexes_end = z_indexes_.end();
01035   int last_z_index = -1;
01036   for ( int k = 0; k < numTotalParams; ++k ) {
01037     if ( z_indexes_itr == z_indexes_end || k < *z_indexes_itr ) {
01038       // This is a "free" parameter subvector
01039       period_l_map_.push_back(k);
01040     }
01041     else {
01042       // This is a "fixed" period parameter subvector so increment
01043       // z_indexes iterator.
01044 #ifdef TEUCHOS_DEBUG
01045       TEST_FOR_EXCEPT( k != *z_indexes_itr && "This should never happen!" );
01046 #endif
01047       const int tmp_last_z_index = *z_indexes_itr;
01048       ++z_indexes_itr;
01049       if ( z_indexes_itr != z_indexes_end ) {
01050 #ifdef TEUCHOS_DEBUG
01051         if ( last_z_index >= 0 ) {
01052           TEST_FOR_EXCEPTION(
01053             *z_indexes_itr <= last_z_index, std::logic_error,
01054             "Error, the z_indexes array = " << toString(z_indexes_)
01055             << " is not sorted or contains duplicate entries!"
01056             );
01057         }
01058 #endif
01059         last_z_index = tmp_last_z_index;
01060       }
01061     }
01062   }
01063 }
01064 
01065 
01066 template<class Scalar>
01067 void DefaultMultiPeriodModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
01068 {
01069 
01070   using Teuchos::rcp_dynamic_cast;
01071   typedef ModelEvaluatorBase MEB;
01072 
01073   nominalValues_ = this->createInArgs();
01074   lowerBounds_ = this->createInArgs();
01075   upperBounds_ = this->createInArgs();
01076 
01077   const MEB::InArgs<Scalar>
01078     periodNominalValues = periodModel_->getNominalValues(),
01079     periodLowerBounds = periodModel_->getLowerBounds(),
01080     periodUpperBounds = periodModel_->getUpperBounds();
01081   
01082   if (periodNominalValues.supports(MEB::IN_ARG_x)) {
01083 
01084     if( !is_null(periodNominalValues.get_x()) ) {
01085       // If the first peroid model has nominal values for x, then all of them
01086       // must also!
01087       RCP<Thyra::ProductVectorBase<Scalar> >
01088         x_bar_init = createProductVector(x_bar_space_);
01089       const int N = this->N();
01090       for ( int i = 0; i < N; ++i ) {
01091         assign(
01092           &*x_bar_init->getNonconstVectorBlock(i),
01093           *periodModels_[i]->getNominalValues().get_x()
01094           );
01095       }
01096       nominalValues_.set_x(x_bar_init);
01097     }
01098       
01099     if( !is_null(periodLowerBounds.get_x()) ) {
01100       // If the first peroid model has lower bounds for for x, then all of
01101       // them must also!
01102       RCP<Thyra::ProductVectorBase<Scalar> >
01103         x_bar_l = createProductVector(x_bar_space_);
01104       const int N = this->N();
01105       for ( int i = 0; i < N; ++i ) {
01106         assign(
01107           &*x_bar_l->getNonconstVectorBlock(i),
01108           *periodModels_[i]->getLowerBounds().get_x()
01109           );
01110       }
01111       lowerBounds_.set_x(x_bar_l);
01112     }
01113       
01114     if( !is_null(periodUpperBounds.get_x()) ) {
01115       // If the first peroid model has upper bounds for for x, then all of
01116       // them must also!
01117       RCP<Thyra::ProductVectorBase<Scalar> >
01118         x_bar_u = createProductVector(x_bar_space_);
01119       const int N = this->N();
01120       for ( int i = 0; i < N; ++i ) {
01121         assign(
01122           &*x_bar_u->getNonconstVectorBlock(i),
01123           *periodModels_[i]->getUpperBounds().get_x()
01124           );
01125       }
01126       upperBounds_.set_x(x_bar_u);
01127     }
01128 
01129   }
01130 
01131   // There can only be one set of nominal values for the non-period parameters
01132   // so just take them from the first period!
01133   for ( int l = 0; l < Np_; ++l ) {
01134     const int period_l = this->period_l(l);
01135     nominalValues_.set_p(l,periodNominalValues.get_p(period_l));
01136     lowerBounds_.set_p(l,periodLowerBounds.get_p(period_l));
01137     upperBounds_.set_p(l,periodUpperBounds.get_p(period_l));
01138   }
01139   
01140 }
01141 
01142 
01143 template<class Scalar>
01144 RCP<ProductVectorBase<Scalar> >
01145 DefaultMultiPeriodModelEvaluator<Scalar>::createProductVector(
01146   const RCP<const ProductVectorSpaceBase<Scalar> > &prodVecSpc
01147   )
01148 {
01149   return Teuchos::rcp_dynamic_cast<ProductVectorBase<Scalar> >(
01150     Thyra::createMember<Scalar>(prodVecSpc), true );
01151 }
01152 
01153 
01154 } // namespace Thyra
01155 
01156 
01157 #endif // THYRA_DEFAULT_MULTI_PERIOD_MODEL_EVALUATOR_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines