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

Generated on Thu Sep 18 12:33:02 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1