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