Thyra_DefaultProductVectorSpace.hpp

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_PRODUCT_VECTOR_SPACE_STD_HPP
00030 #define THYRA_PRODUCT_VECTOR_SPACE_STD_HPP
00031 
00032 #include "Thyra_DefaultProductVectorSpaceDecl.hpp"
00033 #include "Thyra_DefaultProductVector.hpp"
00034 #include "Thyra_DefaultProductMultiVector.hpp"
00035 #include "Thyra_ProductMultiVectorBase.hpp"
00036 #include "Teuchos_dyn_cast.hpp"
00037 
00038 namespace Thyra {
00039 
00040 // Constructors/initializers/accessors
00041 
00042 template<class Scalar>
00043 DefaultProductVectorSpace<Scalar>::DefaultProductVectorSpace(
00044   const int                                                       numBlocks
00045   ,const Teuchos::RCP<const VectorSpaceBase<Scalar> >     vecSpaces[]
00046   )
00047 {
00048   initialize(numBlocks,vecSpaces);
00049 }
00050 
00051 template<class Scalar>
00052 void DefaultProductVectorSpace<Scalar>::initialize(
00053   const int                                                       numBlocks
00054   ,const Teuchos::RCP<const VectorSpaceBase<Scalar> >     vecSpaces[]
00055   )
00056 {
00057   //
00058   // Check preconditions and compute cached quantities
00059   //
00060   TEST_FOR_EXCEPT( numBlocks < 0 );
00061   TEST_FOR_EXCEPT( vecSpaces == NULL );
00062   bool  hasInCoreView = true;
00063   for( int k = 0; k < numBlocks; ++k ) {
00064     TEST_FOR_EXCEPTION(
00065       vecSpaces[k].get() == NULL, std::invalid_argument
00066       ,"Error, the smart pointer vecSpaces["<<k<<"] can not be NULL!"
00067       );
00068     if( !vecSpaces[k]->hasInCoreView() ) hasInCoreView = false;
00069   }
00070   //
00071   // Setup private data members (should not throw an exception from here)
00072   //
00073   numBlocks_ = numBlocks;
00074   vecSpaces_ = Teuchos::rcp(new vecSpaces_t(numBlocks));
00075   std::copy( vecSpaces, vecSpaces+numBlocks, &(*vecSpaces_)[0] );
00076   vecSpacesOffsets_ = Teuchos::rcp(new vecSpacesOffsets_t(numBlocks+1));
00077   (*vecSpacesOffsets_)[0] = 0;
00078   dim_ = 0;
00079   for( int k = 1; k <= numBlocks; ++k ) {
00080     const Index dim_km1 = vecSpaces[k-1]->dim();
00081     (*vecSpacesOffsets_)[k] = (*vecSpacesOffsets_)[k-1] + dim_km1;
00082     dim_ += dim_km1;
00083   }
00084   isInCore_ = hasInCoreView;
00085 }
00086 
00087 template<class Scalar>
00088 void DefaultProductVectorSpace<Scalar>::uninitialize(
00089   int                                                             *numBlocks
00090   ,Teuchos::RCP<const VectorSpaceBase<Scalar> >           vecSpaces[]
00091   )
00092 {
00093   vecSpaces_        = Teuchos::null;
00094   vecSpacesOffsets_ = Teuchos::null;
00095   dim_              = 0;
00096   isInCore_         = false;
00097 }
00098 
00099 template<class Scalar>
00100 void DefaultProductVectorSpace<Scalar>::getVecSpcPoss(
00101   Index i, int* kth_vector_space, Index* kth_global_offset
00102   ) const
00103 {
00104   // Validate the preconditions
00105 #ifdef TEUCHOS_DEBUG
00106   TEST_FOR_EXCEPTION(
00107     !(0 <= i && i < this->dim()), std::out_of_range
00108     ,"VectorSpaceBlocked::get_vector_space_position(...): Error, i = "
00109     << i << " is not in range [0,"<<(this->dim()-1)<<"]"
00110     );
00111 #endif
00112   *kth_vector_space  = 0;
00113   *kth_global_offset = 0;
00114   while( *kth_vector_space < numBlocks_ ) {
00115     const Index off_kp1 = (*vecSpacesOffsets_)[*kth_vector_space+1];
00116     if( off_kp1 > i ) {
00117       *kth_global_offset = (*vecSpacesOffsets_)[*kth_vector_space];
00118       break;
00119     }
00120     ++(*kth_vector_space);
00121   }
00122   TEST_FOR_EXCEPT( !(*kth_vector_space < numBlocks_) );
00123 }
00124 
00125 // Overridden from DefaultProductVectorSpace
00126 
00127 template<class Scalar>
00128 int DefaultProductVectorSpace<Scalar>::numBlocks() const
00129 {
00130   return numBlocks_;
00131 }
00132 
00133 template<class Scalar>
00134 Teuchos::RCP<const VectorSpaceBase<Scalar> >
00135 DefaultProductVectorSpace<Scalar>::getBlock(const int k) const
00136 {
00137   TEST_FOR_EXCEPT( k < 0 || numBlocks_ < k );
00138   return (*vecSpaces_)[k];
00139 }
00140 
00141 // Overridden from VectorSpaceBase
00142 
00143 template<class Scalar>
00144 Index DefaultProductVectorSpace<Scalar>::dim() const
00145 {
00146   return dim_;
00147 }
00148 
00149 template<class Scalar>
00150 bool DefaultProductVectorSpace<Scalar>::isCompatible( const VectorSpaceBase<Scalar>& vecSpc ) const
00151 {
00152   // Check for in-core
00153   if( this->hasInCoreView(Range1D(),VIEW_TYPE_DETACHED,STRIDE_TYPE_NONUNIT) && vecSpc.hasInCoreView() && ( this->dim() == vecSpc.dim() ) )
00154     return true;
00155   // Check for product vector interface
00156   const ProductVectorSpaceBase<Scalar> *pvsb = dynamic_cast<const ProductVectorSpaceBase<Scalar>*>(&vecSpc);
00157   if( !pvsb )
00158     return false;
00159   // Validate that constituent vector spaces are compatible
00160   const int numBlocks = this->numBlocks(); 
00161   if( numBlocks != pvsb->numBlocks() )
00162     return false;
00163   for( int i = 0; i < numBlocks; ++i ) {
00164     if( !this->getBlock(i)->isCompatible(*pvsb->getBlock(i)) )
00165       return false;
00166   }
00167   return true;
00168 }
00169 
00170 template<class Scalar>
00171 Teuchos::RCP< VectorBase<Scalar> >
00172 DefaultProductVectorSpace<Scalar>::createMember() const
00173 {
00174   return Teuchos::rcp(
00175     new DefaultProductVector<Scalar>(Teuchos::rcp(this,false))
00176     );
00177 }
00178 
00179 template<class Scalar>
00180 Scalar DefaultProductVectorSpace<Scalar>::scalarProd(
00181   const VectorBase<Scalar>   &x_in
00182   ,const VectorBase<Scalar>  &y_in
00183   ) const
00184 {
00185   const int numBlocks = this->numBlocks(); 
00186   const ProductVectorBase<Scalar>
00187     &x = Teuchos::dyn_cast<const ProductVectorBase<Scalar> >(x_in),
00188     &y = Teuchos::dyn_cast<const ProductVectorBase<Scalar> >(y_in);
00189 #ifdef TEUCHOS_DEBUG
00190   TEST_FOR_EXCEPT(
00191     numBlocks!=x.productSpace()->numBlocks()
00192     || numBlocks!=y.productSpace()->numBlocks()
00193     );
00194 #endif
00195   Scalar scalarProd = Teuchos::ScalarTraits<Scalar>::zero();
00196   for( int k = 0; k < numBlocks; ++k )
00197     scalarProd += (*vecSpaces_)[k]->scalarProd(
00198       *x.getVectorBlock(k),*y.getVectorBlock(k)
00199       );
00200   return scalarProd;
00201 }
00202 
00203 template<class Scalar>
00204 void DefaultProductVectorSpace<Scalar>::scalarProds(
00205   const MultiVectorBase<Scalar>    &X_in
00206   ,const MultiVectorBase<Scalar>   &Y_in
00207   ,Scalar                          scalar_prods[]
00208   ) const
00209 {
00210   using Teuchos::Workspace;
00211   const VectorSpaceBase<Scalar> &domain = *X_in.domain();
00212   const Index m = domain.dim();
00213 #ifdef TEUCHOS_DEBUG
00214   TEST_FOR_EXCEPT(scalar_prods==NULL);
00215   TEST_FOR_EXCEPT( !domain.isCompatible(*Y_in.domain()) );
00216 #endif
00217   if(m==1) {
00218     scalar_prods[0] = this->scalarProd(*X_in.col(0),*Y_in.col(0));
00219     return;
00220     // ToDo: Remove this if(...) block once we have a DefaultProductMultiVector implementation!
00221   }
00222   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00223   const int numBlocks = this->numBlocks(); 
00224   const ProductMultiVectorBase<Scalar>
00225     &X = Teuchos::dyn_cast<const ProductMultiVectorBase<Scalar> >(X_in),
00226     &Y = Teuchos::dyn_cast<const ProductMultiVectorBase<Scalar> >(Y_in);
00227 #ifdef TEUCHOS_DEBUG
00228   TEST_FOR_EXCEPT( numBlocks!=X.productSpace()->numBlocks() || numBlocks!=Y.productSpace()->numBlocks() );
00229 #endif
00230   Workspace<Scalar> _scalar_prods(wss,m,false);
00231   std::fill_n( scalar_prods, m, Teuchos::ScalarTraits<Scalar>::zero() );
00232   for( int k = 0; k < numBlocks; ++k ) {
00233     (*vecSpaces_)[k]->scalarProds(*X.getMultiVectorBlock(k),*Y.getMultiVectorBlock(k),&_scalar_prods[0]);
00234     for( int j = 0; j < m; ++j ) scalar_prods[j] += _scalar_prods[j];
00235   }
00236 }
00237 
00238 template<class Scalar>
00239 bool DefaultProductVectorSpace<Scalar>::hasInCoreView(const Range1D& rng_in, const EViewType viewType, const EStrideType strideType) const
00240 {
00241   const Range1D rng = full_range(rng_in,0,dim_-1);
00242   // First see if rng fits in a single constituent vector
00243   int    kth_vector_space  = -1;
00244   Index  kth_global_offset = 0;
00245   this->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
00246 #ifdef TEUCHOS_DEBUG
00247   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00248 #endif
00249   if( rng.lbound() + rng.size() <= kth_global_offset + (*vecSpaces_)[kth_vector_space]->dim() ) {
00250     return (*vecSpaces_)[kth_vector_space]->hasInCoreView(rng_in-kth_global_offset,viewType,strideType);
00251   }
00252   // If we get here, rng does not fit in a single constituent vector which
00253   // also means that numBlocks_ > 1 must also be true!
00254   //
00255   // Next, if the client is asking for a direct view then we have to return
00256   // false since this range spans more than one constituent vector.
00257   if( viewType == VIEW_TYPE_DIRECT )
00258     return false;
00259   // If we get here then hasDirectView==false and therefore we are allowed to
00260   // create a copy.  Therefore, if all of the constituent vectors are "in
00261   // core" then we can return true.
00262   if(isInCore_)
00263     return true;
00264   // Finally, loop through all of the constituent vectors spaned by rng and
00265   // see if they are each in core.
00266   //
00267   // Todo: Implement this if you have to!
00268   //
00269   // We must give up and return false
00270   return false;
00271 }
00272 
00273 template<class Scalar>
00274 Teuchos::RCP< const VectorSpaceFactoryBase<Scalar> >
00275 DefaultProductVectorSpace<Scalar>::smallVecSpcFcty() const
00276 {
00277   if( dim_ ) return (*vecSpaces_)[0]->smallVecSpcFcty(); // They should all be compatible!
00278   return Teuchos::null;
00279 }
00280 
00281 template<class Scalar>
00282 Teuchos::RCP< MultiVectorBase<Scalar> >
00283 DefaultProductVectorSpace<Scalar>::createMembers(int numMembers) const
00284 {
00285   return defaultProductMultiVector<Scalar>(Teuchos::rcp(this,false),numMembers);
00286 /*
00287   if(numMembers==1)
00288     return createMember();
00289   return VectorSpaceDefaultBase<Scalar>::createMembers(numMembers);
00290   // ToDo: Specialize to return DefaultProductMultiVector when needed!
00291   */
00292 }
00293 
00294 template<class Scalar>
00295 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00296 DefaultProductVectorSpace<Scalar>::clone() const
00297 {
00298   // Warning! If the client uninitialized this object then changes the
00299   // constituent vector spaces then we are in trouble!  The client is warned
00300   // in documentation!
00301   Teuchos::RCP<DefaultProductVectorSpace<Scalar> >
00302     pvs = Teuchos::rcp(new DefaultProductVectorSpace<Scalar>());
00303   pvs->numBlocks_          = numBlocks_;
00304   pvs->vecSpaces_          = vecSpaces_;
00305   pvs->vecSpacesOffsets_   = vecSpacesOffsets_;
00306   pvs->dim_                = dim_;
00307   pvs->isInCore_           = isInCore_;
00308   return pvs;
00309 }
00310 
00311 
00312 // Overridden from Teuchos::Describable
00313 
00314                                                 
00315 template<class Scalar>
00316 std::string DefaultProductVectorSpace<Scalar>::description() const
00317 {
00318   assertInitialized();
00319   std::ostringstream oss;
00320   oss
00321     << Teuchos::Describable::description() << "{"
00322     << "dim="<<dim_
00323     << ",numBlocks="<<numBlocks_
00324     << "}";
00325   return oss.str();
00326 }
00327 
00328 
00329 template<class Scalar>
00330 void DefaultProductVectorSpace<Scalar>::describe(
00331   Teuchos::FancyOStream                &out_arg
00332   ,const Teuchos::EVerbosityLevel      verbLevel
00333   ) const
00334 {
00335   typedef Teuchos::ScalarTraits<Scalar>  ST;
00336   using Teuchos::RCP;
00337   using Teuchos::FancyOStream;
00338   using Teuchos::OSTab;
00339   assertInitialized();
00340   RCP<FancyOStream> out = rcp(&out_arg,false);
00341   OSTab tab(out);
00342   switch(verbLevel) {
00343     case Teuchos::VERB_DEFAULT:
00344     case Teuchos::VERB_LOW:
00345       *out << this->description() << std::endl;
00346       break;
00347     case Teuchos::VERB_MEDIUM:
00348     case Teuchos::VERB_HIGH:
00349     case Teuchos::VERB_EXTREME:
00350     {
00351       *out
00352         << this->description() << std::endl;
00353       OSTab tab(out);
00354       *out
00355         <<  "Constituent vector spaces V[0], V[1], ... V[numBlocks-1]:\n";
00356       tab.incrTab();
00357       for( int k = 0; k < numBlocks_; ++k ) {
00358         *out << "V["<<k<<"] = " << Teuchos::describe(*(*vecSpaces_)[k],verbLevel);
00359       }
00360       break;
00361     }
00362     default:
00363       TEST_FOR_EXCEPT(true); // Should never get here!
00364   }
00365 }
00366 
00367 
00368 } // namespace Thyra
00369 
00370 
00371 #endif // THYRA_PRODUCT_VECTOR_SPACE_STD_HPP

Generated on Wed May 12 21:42:27 2010 for Thyra Operator/Vector Support by  doxygen 1.4.7