Thyra Version of the Day
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef THYRA_DEFAULT_PRODUCT_VECTOR_SPACE_HPP
00043 #define THYRA_DEFAULT_PRODUCT_VECTOR_SPACE_HPP
00044 
00045 
00046 #include "Thyra_DefaultProductVectorSpace_decl.hpp"
00047 #include "Thyra_DefaultProductVector.hpp"
00048 #include "Thyra_DefaultProductMultiVector.hpp"
00049 #include "Thyra_ProductMultiVectorBase.hpp"
00050 #include "Teuchos_Workspace.hpp"
00051 #include "Teuchos_dyn_cast.hpp"
00052 
00053 
00054 namespace Thyra {
00055 
00056 
00057 // Constructors/initializers/accessors
00058 
00059 
00060 template<class Scalar>
00061 DefaultProductVectorSpace<Scalar>::DefaultProductVectorSpace()
00062   : numBlocks_(-1), dim_(-1)
00063 {}
00064 
00065 
00066 template<class Scalar>
00067 DefaultProductVectorSpace<Scalar>::DefaultProductVectorSpace(
00068   const ArrayView<const RCP<const VectorSpaceBase<Scalar> > > &vecSpaces_in
00069   )
00070   : numBlocks_(-1), dim_(-1)
00071 {
00072   initialize(vecSpaces_in);
00073 }
00074 
00075 
00076 template<class Scalar>
00077 void DefaultProductVectorSpace<Scalar>::initialize(
00078   const ArrayView<const RCP<const VectorSpaceBase<Scalar> > > &vecSpaces_in
00079   )
00080 {
00081 
00082   //
00083   // Check preconditions and compute cached quantities
00084   //
00085   const int nBlocks = vecSpaces_in.size();
00086 #ifdef TEUCHOS_DEBUG
00087   TEUCHOS_TEST_FOR_EXCEPT( nBlocks == 0 );
00088 #endif
00089   bool overallHasInCoreView = true;
00090   for (int k = 0; k < nBlocks; ++k) {
00091 #ifdef TEUCHOS_DEBUG
00092     TEUCHOS_TEST_FOR_EXCEPTION(
00093       vecSpaces_in[k].get() == NULL, std::invalid_argument
00094       ,"Error, the smart pointer vecSpaces["<<k<<"] can not be NULL!"
00095       );
00096 #endif
00097     if (!vecSpaces_in[k]->hasInCoreView()) overallHasInCoreView = false;
00098   }
00099 
00100   //
00101   // Setup private data members (should not throw an exception from here)
00102   //
00103   numBlocks_ = nBlocks;
00104   vecSpaces_ = Teuchos::rcp(new vecSpaces_t);
00105   *vecSpaces_ = vecSpaces_in;
00106   vecSpacesOffsets_ = Teuchos::rcp(new vecSpacesOffsets_t(nBlocks+1));
00107   (*vecSpacesOffsets_)[0] = 0;
00108   dim_ = 0;
00109   for( int k = 1; k <= nBlocks; ++k ) {
00110     const Ordinal dim_km1 = vecSpaces_in[k-1]->dim();
00111     (*vecSpacesOffsets_)[k] = (*vecSpacesOffsets_)[k-1] + dim_km1;
00112     dim_ += dim_km1;
00113   }
00114   isInCore_ = overallHasInCoreView;
00115 
00116 }
00117 
00118 
00119 template<class Scalar>
00120 void DefaultProductVectorSpace<Scalar>::uninitialize(
00121   const ArrayView<RCP<const VectorSpaceBase<Scalar> > > &vecSpaces_in
00122   )
00123 {
00124   TEUCHOS_TEST_FOR_EXCEPT(!is_null(vecSpaces_in)); // ToDo: Implement!
00125   vecSpaces_ = Teuchos::null;
00126   vecSpacesOffsets_ = Teuchos::null;
00127   numBlocks_ = -1;
00128   dim_ = -1;
00129   isInCore_ = false;
00130 }
00131 
00132 
00133 template<class Scalar>
00134 void DefaultProductVectorSpace<Scalar>::getVecSpcPoss(
00135   Ordinal i, int* kth_vector_space, Ordinal* kth_global_offset
00136   ) const
00137 {
00138   // Validate the preconditions
00139 #ifdef TEUCHOS_DEBUG
00140   TEUCHOS_TEST_FOR_EXCEPTION(
00141     !(0 <= i && i < this->dim()), std::out_of_range
00142     ,"VectorSpaceBlocked::get_vector_space_position(...): Error, i = "
00143     << i << " is not in range [0,"<<(this->dim()-1)<<"]"
00144     );
00145 #endif
00146   *kth_vector_space  = 0;
00147   *kth_global_offset = 0;
00148   while( *kth_vector_space < numBlocks_ ) {
00149     const Ordinal off_kp1 = (*vecSpacesOffsets_)[*kth_vector_space+1];
00150     if( off_kp1 > i ) {
00151       *kth_global_offset = (*vecSpacesOffsets_)[*kth_vector_space];
00152       break;
00153     }
00154     ++(*kth_vector_space);
00155   }
00156   TEUCHOS_TEST_FOR_EXCEPT( !(*kth_vector_space < numBlocks_) );
00157 }
00158 
00159 
00160 // Overridden from DefaultProductVectorSpace
00161 
00162 
00163 template<class Scalar>
00164 int DefaultProductVectorSpace<Scalar>::numBlocks() const
00165 {
00166   return numBlocks_;
00167 }
00168 
00169 
00170 template<class Scalar>
00171 Teuchos::RCP<const VectorSpaceBase<Scalar> >
00172 DefaultProductVectorSpace<Scalar>::getBlock(const int k) const
00173 {
00174   TEUCHOS_TEST_FOR_EXCEPT( k < 0 || numBlocks_ < k );
00175   return (*vecSpaces_)[k];
00176 }
00177 
00178 
00179 // Overridden from VectorSpaceBase
00180 
00181 
00182 template<class Scalar>
00183 Ordinal DefaultProductVectorSpace<Scalar>::dim() const
00184 {
00185   return dim_;
00186 }
00187 
00188 
00189 template<class Scalar>
00190 bool DefaultProductVectorSpace<Scalar>::isCompatible(
00191   const VectorSpaceBase<Scalar>& vecSpc ) const
00192 {
00193 
00194   using Teuchos::ptrFromRef;
00195   using Teuchos::ptr_dynamic_cast;
00196 
00197   const int nBlocks = this->numBlocks();
00198 
00199   // Check for product vector interface
00200   const Ptr<const ProductVectorSpaceBase<Scalar> > pvsb =
00201     ptr_dynamic_cast<const ProductVectorSpaceBase<Scalar> >(ptrFromRef(vecSpc));
00202 
00203   if (nonnull(pvsb)) {
00204     // Validate that constituent vector spaces are compatible
00205     if( nBlocks != pvsb->numBlocks() )
00206       return false;
00207     for( int i = 0; i < nBlocks; ++i ) {
00208       if( !this->getBlock(i)->isCompatible(*pvsb->getBlock(i)) )
00209         return false;
00210     }
00211     return true;
00212   }
00213 
00214   // Check for a single vector single vector space
00215   if (nBlocks == 1) {
00216     return this->getBlock(0)->isCompatible(vecSpc);
00217   }
00218 
00219   // If we get here, the RHS is not a product vector space and/or this is not
00220   // a single block VS so we can assume the spaces are *not* compatible!
00221   return false;
00222 
00223 }
00224 
00225 
00226 template<class Scalar>
00227 Teuchos::RCP< VectorBase<Scalar> >
00228 DefaultProductVectorSpace<Scalar>::createMember() const
00229 {
00230   return defaultProductVector<Scalar>(Teuchos::rcpFromRef(*this));
00231 }
00232 
00233 
00234 template<class Scalar>
00235 Scalar DefaultProductVectorSpace<Scalar>::scalarProd(
00236   const VectorBase<Scalar> &x_in,
00237   const VectorBase<Scalar> &y_in
00238   ) const
00239 {
00240   const int nBlocks = this->numBlocks(); 
00241   const ProductVectorBase<Scalar>
00242     &x = Teuchos::dyn_cast<const ProductVectorBase<Scalar> >(x_in),
00243     &y = Teuchos::dyn_cast<const ProductVectorBase<Scalar> >(y_in);
00244 #ifdef TEUCHOS_DEBUG
00245   TEUCHOS_TEST_FOR_EXCEPT(
00246     nBlocks!=x.productSpace()->numBlocks()
00247     || nBlocks!=y.productSpace()->numBlocks()
00248     );
00249 #endif
00250   Scalar scalarProd_rtn = Teuchos::ScalarTraits<Scalar>::zero();
00251   for( int k = 0; k < nBlocks; ++k )
00252     scalarProd_rtn += (*vecSpaces_)[k]->scalarProd(
00253       *x.getVectorBlock(k),*y.getVectorBlock(k)
00254       );
00255   return scalarProd_rtn;
00256 }
00257 
00258 
00259 template<class Scalar>
00260 void DefaultProductVectorSpace<Scalar>::scalarProdsImpl(
00261   const MultiVectorBase<Scalar> &X_in,
00262  const MultiVectorBase<Scalar> &Y_in,
00263   const ArrayView<Scalar> &scalarProds_out
00264   ) const
00265 {
00266   using Teuchos::as;
00267   using Teuchos::Workspace;
00268   const VectorSpaceBase<Scalar> &domain = *X_in.domain();
00269   const Ordinal m = domain.dim();
00270 #ifdef TEUCHOS_DEBUG
00271   TEUCHOS_TEST_FOR_EXCEPT(is_null(scalarProds_out));
00272   TEUCHOS_TEST_FOR_EXCEPT( !domain.isCompatible(*Y_in.domain()) );
00273   TEUCHOS_ASSERT_EQUALITY( as<Ordinal>(scalarProds_out.size()),
00274     as<Ordinal>(m) )
00275 #endif
00276   if(m==1) {
00277     scalarProds_out[0] = this->scalarProd(*X_in.col(0),*Y_in.col(0));
00278     return;
00279     // ToDo: Remove this if(...) block once we have a DefaultProductMultiVector implementation!
00280   }
00281   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00282   const int nBlocks = this->numBlocks(); 
00283   const ProductMultiVectorBase<Scalar>
00284     &X = Teuchos::dyn_cast<const ProductMultiVectorBase<Scalar> >(X_in),
00285     &Y = Teuchos::dyn_cast<const ProductMultiVectorBase<Scalar> >(Y_in);
00286 #ifdef TEUCHOS_DEBUG
00287   TEUCHOS_TEST_FOR_EXCEPT( nBlocks!=X.productSpace()->numBlocks() || nBlocks!=Y.productSpace()->numBlocks() );
00288 #endif
00289   Workspace<Scalar> _scalarProds_out(wss, m, false);
00290   std::fill( scalarProds_out.begin(), scalarProds_out.end(),
00291     ScalarTraits<Scalar>::zero() );
00292   for( int k = 0; k < nBlocks; ++k ) {
00293     (*vecSpaces_)[k]->scalarProds(
00294       *X.getMultiVectorBlock(k), *Y.getMultiVectorBlock(k), _scalarProds_out());
00295     for( int j = 0; j < m; ++j )
00296       scalarProds_out[j] += _scalarProds_out[j];
00297   }
00298 }
00299 
00300 
00301 template<class Scalar>
00302 bool DefaultProductVectorSpace<Scalar>::hasInCoreView(const Range1D& rng_in, const EViewType viewType, const EStrideType strideType) const
00303 {
00304   const Range1D rng = full_range(rng_in,0,dim_-1);
00305   // First see if rng fits in a single constituent vector
00306   int    kth_vector_space  = -1;
00307   Ordinal  kth_global_offset = 0;
00308   this->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
00309 #ifdef TEUCHOS_DEBUG
00310   TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00311 #endif
00312   if( rng.lbound() + rng.size() <= kth_global_offset + (*vecSpaces_)[kth_vector_space]->dim() ) {
00313     return (*vecSpaces_)[kth_vector_space]->hasInCoreView(rng_in-kth_global_offset,viewType,strideType);
00314   }
00315   // If we get here, rng does not fit in a single constituent vector which
00316   // also means that numBlocks_ > 1 must also be true!
00317   //
00318   // Next, if the client is asking for a direct view then we have to return
00319   // false since this range spans more than one constituent vector.
00320   if( viewType == VIEW_TYPE_DIRECT )
00321     return false;
00322   // If we get here then hasDirectView==false and therefore we are allowed to
00323   // create a copy.  Therefore, if all of the constituent vectors are "in
00324   // core" then we can return true.
00325   if(isInCore_)
00326     return true;
00327   // Finally, loop through all of the constituent vectors spaned by rng and
00328   // see if they are each in core.
00329   //
00330   // Todo: Implement this if you have to!
00331   //
00332   // We must give up and return false
00333   return false;
00334 }
00335 
00336 
00337 template<class Scalar>
00338 Teuchos::RCP< const VectorSpaceFactoryBase<Scalar> >
00339 DefaultProductVectorSpace<Scalar>::smallVecSpcFcty() const
00340 {
00341   if (dim_)
00342     return (*vecSpaces_)[0]->smallVecSpcFcty(); // They should all be compatible?
00343   return Teuchos::null;
00344 }
00345 
00346 
00347 template<class Scalar>
00348 Teuchos::RCP< MultiVectorBase<Scalar> >
00349 DefaultProductVectorSpace<Scalar>::createMembers(int numMembers) const
00350 {
00351   return defaultProductMultiVector<Scalar>(Teuchos::rcpFromRef(*this),
00352     numMembers);
00353 }
00354 
00355 
00356 template<class Scalar>
00357 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00358 DefaultProductVectorSpace<Scalar>::clone() const
00359 {
00360   // Warning! If the client uninitialized this object then changes the
00361   // constituent vector spaces then we are in trouble!  The client is warned
00362   // in documentation!
00363   Teuchos::RCP<DefaultProductVectorSpace<Scalar> >
00364     pvs = productVectorSpace<Scalar>();
00365   pvs->numBlocks_          = numBlocks_;
00366   pvs->vecSpaces_          = vecSpaces_;
00367   pvs->vecSpacesOffsets_   = vecSpacesOffsets_;
00368   pvs->dim_                = dim_;
00369   pvs->isInCore_           = isInCore_;
00370   return pvs;
00371 }
00372 
00373 
00374 // Overridden from Teuchos::Describable
00375 
00376                                                 
00377 template<class Scalar>
00378 std::string DefaultProductVectorSpace<Scalar>::description() const
00379 {
00380   std::ostringstream oss;
00381   oss
00382     << Teuchos::Describable::description() << "{"
00383     << "dim="<<dim_
00384     << ",numBlocks="<<numBlocks_
00385     << "}";
00386   return oss.str();
00387 }
00388 
00389 
00390 template<class Scalar>
00391 void DefaultProductVectorSpace<Scalar>::describe(
00392   Teuchos::FancyOStream                &out_arg
00393   ,const Teuchos::EVerbosityLevel      verbLevel
00394   ) const
00395 {
00396   typedef Teuchos::ScalarTraits<Scalar>  ST;
00397   using Teuchos::includesVerbLevel;
00398   using Teuchos::OSTab;
00399   RCP<FancyOStream> out = rcpFromRef(out_arg);
00400   OSTab tab(out);
00401   if (includesVerbLevel(verbLevel, Teuchos::VERB_LOW, true)) {
00402     *out << this->description() << std::endl;
00403   }
00404   if (includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM) && numBlocks_ > 0) {
00405     OSTab tab2(out);
00406     *out <<  "Constituent vector spaces V[0], V[1], ... V[numBlocks-1]:\n";
00407     OSTab tab3(out);
00408     for( int k = 0; k < numBlocks_; ++k ) {
00409       *out << "V["<<k<<"] = " << Teuchos::describe(*(*vecSpaces_)[k],verbLevel);
00410     }
00411   }
00412 }
00413 
00414 
00415 } // namespace Thyra
00416 
00417 
00418 #endif // THYRA_DEFAULT_PRODUCT_VECTOR_SPACE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines