Thyra_DefaultProductVector.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_HPP
00030 #define THYRA_PRODUCT_VECTOR_HPP
00031 
00032 #include "Thyra_DefaultProductVectorDecl.hpp"
00033 #include "Thyra_DefaultProductVectorSpace.hpp"
00034 #include "Teuchos_Workspace.hpp"
00035 
00036 namespace Thyra {
00037 
00038 // Constructors/initializers/accessors
00039 
00040 template <class Scalar>
00041 DefaultProductVector<Scalar>::DefaultProductVector()
00042 {
00043   uninitialize();
00044 }
00045 
00046 template <class Scalar>
00047 DefaultProductVector<Scalar>::DefaultProductVector(
00048   const Teuchos::RefCountPtr<const DefaultProductVectorSpace<Scalar> >  &productSpace
00049   )
00050 {
00051   initialize(productSpace);
00052 }
00053 
00054 template <class Scalar>
00055 DefaultProductVector<Scalar>::DefaultProductVector(
00056   const Teuchos::RefCountPtr<const DefaultProductVectorSpace<Scalar> >  &productSpace
00057   ,const Teuchos::RefCountPtr<VectorBase<Scalar> >                      vecs[]
00058   )
00059 {
00060   initialize(productSpace,vecs);
00061 }
00062 
00063 template <class Scalar>
00064 DefaultProductVector<Scalar>::DefaultProductVector(
00065   const Teuchos::RefCountPtr<const DefaultProductVectorSpace<Scalar> >  &productSpace
00066   ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >                vecs[]
00067   )
00068 {
00069   initialize(productSpace,vecs);
00070 }
00071 
00072 template <class Scalar>
00073 void DefaultProductVector<Scalar>::initialize(
00074   const Teuchos::RefCountPtr<const DefaultProductVectorSpace<Scalar> >  &productSpace
00075   )
00076 {
00077   // ToDo: Validate input!
00078   numBlocks_ = productSpace->numBlocks();
00079   productSpace_ = productSpace;
00080   vecs_.resize(numBlocks_);
00081   for( int k = 0; k < numBlocks_; ++k )
00082     vecs_[k].initialize(createMember(productSpace->getBlock(k)));
00083 }
00084 
00085 template <class Scalar>
00086 void DefaultProductVector<Scalar>::initialize(
00087   const Teuchos::RefCountPtr<const DefaultProductVectorSpace<Scalar> >  &productSpace
00088   ,const Teuchos::RefCountPtr<VectorBase<Scalar> >                      vecs[]
00089   )
00090 {
00091   // ToDo: Validate input!
00092   numBlocks_ = productSpace->numBlocks();
00093   productSpace_ = productSpace;
00094   vecs_.resize(numBlocks_);
00095   for( int k = 0; k < numBlocks_; ++k )
00096     vecs_[k].initialize(vecs[k]);
00097 }
00098 
00099 template <class Scalar>
00100 void DefaultProductVector<Scalar>::initialize(
00101   const Teuchos::RefCountPtr<const DefaultProductVectorSpace<Scalar> >  &productSpace
00102   ,const Teuchos::RefCountPtr<const VectorBase<Scalar> >                vecs[]
00103   )
00104 {
00105   // ToDo: Validate input!
00106   numBlocks_ = productSpace->numBlocks();
00107   productSpace_ = productSpace;
00108   vecs_.resize(numBlocks_);
00109   for( int k = 0; k < numBlocks_; ++k )
00110     vecs_[k].initialize(vecs[k]);
00111 }
00112 
00113 template <class Scalar>
00114 void DefaultProductVector<Scalar>::uninitialize()
00115 {
00116   productSpace_ = Teuchos::null;
00117   vecs_.resize(0);
00118   numBlocks_ = 0;
00119 }
00120 
00121 // Overridden from ProductVectorBase
00122 
00123 template <class Scalar>
00124 Teuchos::RefCountPtr<VectorBase<Scalar> >
00125 DefaultProductVector<Scalar>::getNonconstVectorBlock(const int k)
00126 {
00127   TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
00128   return vecs_[k].getNonconstObj();
00129 }
00130 
00131 template <class Scalar>
00132 Teuchos::RefCountPtr<const VectorBase<Scalar> >
00133 DefaultProductVector<Scalar>::getVectorBlock(const int k) const
00134 {
00135   TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
00136   return vecs_[k].getConstObj();
00137 }
00138 
00139 // Overridden from ProductMultiVectorBase
00140 
00141 template <class Scalar>
00142 Teuchos::RefCountPtr<const ProductVectorSpaceBase<Scalar> >
00143 DefaultProductVector<Scalar>::productSpace() const
00144 {
00145   return productSpace_;
00146 }
00147 
00148 template <class Scalar>
00149 bool DefaultProductVector<Scalar>::blockIsConst(const int k) const
00150 {
00151   TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
00152   return vecs_[k].isConst();
00153 }
00154 
00155 template <class Scalar>
00156 Teuchos::RefCountPtr<MultiVectorBase<Scalar> >
00157 DefaultProductVector<Scalar>::getNonconstMultiVectorBlock(const int k)
00158 {
00159   return getNonconstVectorBlock(k);
00160 }
00161 
00162 template <class Scalar>
00163 void DefaultProductVector<Scalar>
00164 ::setBlock(int i, 
00165            const Teuchos::RefCountPtr<const VectorBase<Scalar> >& b)
00166 {
00167   TEST_FOR_EXCEPT(i < 0 || i >= numBlocks_);
00168   TEST_FOR_EXCEPT(!productSpace_->getBlock(i)->isCompatible(*(b->space())));
00169   vecs_[i] = b;
00170 }
00171 
00172 template <class Scalar>
00173 void DefaultProductVector<Scalar>
00174 ::setNonconstBlock(int i, 
00175                    const Teuchos::RefCountPtr<VectorBase<Scalar> >& b)
00176 {
00177   TEST_FOR_EXCEPT(i < 0 || i >= numBlocks_);
00178   TEST_FOR_EXCEPT(!productSpace_->getBlock(i)->isCompatible(*(b->space())));
00179   vecs_[i] = b;
00180 }
00181 
00182 template <class Scalar>
00183 Teuchos::RefCountPtr<const MultiVectorBase<Scalar> >
00184 DefaultProductVector<Scalar>::getMultiVectorBlock(const int k) const
00185 {
00186   return getVectorBlock(k);
00187 }
00188 
00189 // Overridden from VectorBase
00190 
00191 template <class Scalar>
00192 Teuchos::RefCountPtr< const VectorSpaceBase<Scalar> >
00193 DefaultProductVector<Scalar>::space() const
00194 {
00195   return productSpace_;
00196 }
00197 
00198 template <class Scalar>
00199 void DefaultProductVector<Scalar>::applyOp(
00200   const RTOpPack::RTOpT<Scalar>    &op
00201   ,const int                       num_vecs
00202   ,const VectorBase<Scalar>*const  vecs[]
00203   ,const int                       num_targ_vecs
00204   ,VectorBase<Scalar>*const        targ_vecs[]
00205   ,RTOpPack::ReductTarget          *reduct_obj
00206   ,const Index                     first_ele_offset_in
00207   ,const Index                     sub_dim_in
00208   ,const Index                     global_offset_in
00209   ) const
00210 {
00211   using Teuchos::Workspace;
00212   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00213   //
00214   const Index n = productSpace_->dim();
00215   // Validate the compatibility of the vectors!
00216 #ifdef TEUCHOS_DEBUG
00217   TEST_FOR_EXCEPTION(
00218     !(0 <= first_ele_offset_in && first_ele_offset_in < n), std::out_of_range
00219     ,"DefaultProductVector::applyOp(...): Error, "
00220     "first_ele_offset_in = " << first_ele_offset_in << " is not in range [0,"<<(n-1)<<"]" );
00221   TEST_FOR_EXCEPTION(
00222     global_offset_in < 0, std::invalid_argument
00223     ,"DefaultProductVector::applyOp(...): Error, "
00224     "global_offset_in = " << global_offset_in << " is not acceptable" );
00225   TEST_FOR_EXCEPTION(
00226     sub_dim_in >= 0 && sub_dim_in - first_ele_offset_in > n, std::length_error
00227     ,"DefaultProductVector::applyOp(...): Error, "
00228     "global_offset_in = " << global_offset_in << ", sub_dim_in = " << sub_dim_in
00229     << "first_ele_offset_in = " << first_ele_offset_in << " and n = " << n
00230     << " are not compatible" );
00231   bool test_failed;
00232   for(int k = 0; k < num_vecs; ++k) {
00233     test_failed = !this->space()->isCompatible(*vecs[k]->space());
00234     TEST_FOR_EXCEPTION(
00235       test_failed, Exceptions::IncompatibleVectorSpaces
00236       ,"DefaultProductVector::applyOp(...): Error vecs["<<k<<"]->space() "
00237       <<"of type \'"<<typeid(*vecs[k]->space()).name()<<"\' is not compatible with this "
00238       <<"\'VectorSpaceBlocked\' vector space!"
00239       );
00240   }
00241   for(int k = 0; k < num_targ_vecs; ++k) {
00242     test_failed = !this->space()->isCompatible(*targ_vecs[k]->space());
00243     TEST_FOR_EXCEPTION(
00244       test_failed, Exceptions::IncompatibleVectorSpaces
00245       ,"DefaultProductVector::applyOp(...): Error targ_vecs["<<k<<"]->space() "
00246       <<"of type \'"<<typeid(*vecs[k]->space()).name()<<"\' is not compatible with this "
00247       <<"\'VectorSpaceBlocked\' vector space!"
00248       );
00249   }
00250 #endif
00251   //
00252   // The first thing that we do is to see if all of the vectors involved are
00253   // incore, serial vectors.  In this case, the vectors should be compatible.
00254   // To accomplish this we will pick the continguous vector to implement
00255   // the applyOp(...) function.
00256   //
00257   // Get the index of an incore-only input vector and input/output vector?
00258   const bool this_isInCore = productSpace_->hasInCoreView(
00259     Range1D(),VIEW_TYPE_DETACHED,STRIDE_TYPE_NONUNIT
00260     );
00261   int incore_vec_k = -1, incore_targ_vec_k = -1;
00262   // Dynamic cast the pointers for the vector arguments
00263   Workspace<const DefaultProductVector<Scalar>*>
00264     vecs_args(wss,num_vecs,false);
00265   for(int k = 0; k < num_vecs; ++k) {
00266     vecs_args[k] = dynamic_cast<const DefaultProductVector<Scalar>*>(vecs[k]);
00267     if( vecs_args[k] == NULL ) {
00268       const bool isInCore_k = vecs[k]->space()->hasInCoreView();
00269       if( this_isInCore && isInCore_k ) {
00270         incore_vec_k = k;
00271         break;
00272       }
00273       TEST_FOR_EXCEPTION(
00274         !this_isInCore || (this_isInCore && !isInCore_k), Exceptions::IncompatibleVectorSpaces
00275         ,"DefaultProductVector::applyOp(...): Error vecs["<<k<<"] "
00276         <<"of type \'"<<typeid(*vecs[k]).name()<<"\' does not support the "
00277         <<"\'DefaultProductVector<Scalar>\' interface and is not an incore vector or this is not an incore vector!"
00278         );
00279     }
00280   }
00281   Workspace<DefaultProductVector<Scalar>*>
00282     targ_vecs_args(wss,num_targ_vecs,false);
00283   for(int k = 0; k < num_targ_vecs; ++k) {
00284     targ_vecs_args[k] = dynamic_cast<DefaultProductVector<Scalar>*>(targ_vecs[k]);
00285     if( targ_vecs_args[k] == NULL ) {
00286       const bool isInCore_k = targ_vecs[k]->space()->hasInCoreView();
00287       if( this_isInCore && isInCore_k ) {
00288         incore_targ_vec_k = k;
00289         break;
00290       }
00291       TEST_FOR_EXCEPTION(
00292         !this_isInCore || (this_isInCore && !isInCore_k), Exceptions::IncompatibleVectorSpaces
00293         ,"DefaultProductVector::applyOp(...): Error targ_vecs["<<k<<"] "
00294         <<"of type \'"<<typeid(*targ_vecs[k]).name()<<"\' does not support the "
00295         <<"\'DefaultProductVector<Scalar>\' interface and is not an incore vector or this is not an incore vector!"
00296         );
00297     }
00298   }
00299   // Let a incore-only vector with a contiguous view handle this through
00300   // explicit vector access?
00301   if( incore_vec_k >= 0 ) {
00302     vecs[incore_vec_k]->applyOp(
00303       op,num_vecs,vecs,num_targ_vecs,targ_vecs,reduct_obj
00304       ,first_ele_offset_in,sub_dim_in,global_offset_in
00305       );
00306     return;
00307   }
00308   else if ( incore_targ_vec_k >= 0 ) {
00309     targ_vecs[incore_targ_vec_k]->applyOp(
00310       op,num_vecs,vecs,num_targ_vecs,targ_vecs,reduct_obj
00311       ,first_ele_offset_in,sub_dim_in,global_offset_in
00312       );
00313     return;
00314   }
00315   //
00316   // If we get here, then we will implement the applyOp(...)  one vector block
00317   // at a time.
00318   //
00319   const Index this_dim = n;
00320   const Index sub_dim  = ( sub_dim_in < 0
00321                            ? this_dim - first_ele_offset_in
00322                            : sub_dim_in );
00323   Index num_elements_remaining = sub_dim;
00324   const int  numBlocks = productSpace_->numBlocks();
00325   Workspace<const VectorBase<Scalar>*>
00326     sub_vecs(wss,num_vecs,false);
00327   Workspace<VectorBase<Scalar>*>
00328     sub_targ_vecs(wss,num_targ_vecs,false);
00329   Index g_off = -first_ele_offset_in;
00330   for(int k = 0; k < numBlocks; ++k) {
00331     const Index local_dim = productSpace_->getBlock(k)->dim();
00332     if( g_off < 0 && -g_off+1 > local_dim ) {
00333       g_off += local_dim;
00334       continue;
00335     }
00336     const Index
00337       local_sub_dim
00338       = ( g_off >= 0
00339           ? std::min( local_dim, num_elements_remaining )
00340           : std::min( local_dim + g_off, num_elements_remaining ) );
00341     if( local_sub_dim <= 0 )
00342       break;
00343     for( int i = 0; i < num_vecs; ++i )       // Fill constituent vectors for block k
00344       sub_vecs[i] = &*vecs_args[i]->vecs_[k].getConstObj();
00345     for( int j = 0; j < num_targ_vecs; ++j )  // Fill constituent target vectors for block k
00346       sub_targ_vecs[j] = &*targ_vecs_args[j]->vecs_[k].getNonconstObj();
00347     Thyra::applyOp( 
00348       op
00349       ,num_vecs,num_vecs?&sub_vecs[0]:NULL
00350       ,num_targ_vecs,num_targ_vecs?&sub_targ_vecs[0]:NULL
00351       ,reduct_obj
00352       ,g_off < 0 ? -g_off : 0                                    // first_ele_offset
00353       ,local_sub_dim                                             // sub_dim
00354       ,g_off < 0 ? global_offset_in : global_offset_in + g_off   // global_offset
00355       );
00356     g_off += local_dim;
00357     num_elements_remaining -= local_sub_dim;
00358   }
00359   TEST_FOR_EXCEPT(!(num_elements_remaining==0));
00360 }
00361 
00362 template <class Scalar>
00363 void DefaultProductVector<Scalar>::acquireDetachedView(
00364   const Range1D& rng_in, RTOpPack::ConstSubVectorView<Scalar>* sub_vec
00365   ) const
00366 {
00367   const Range1D
00368     rng = rng_in.full_range() ? Range1D(0,productSpace_->dim()-1) : rng_in;
00369   int    kth_vector_space  = -1;
00370   Index  kth_global_offset = 0;
00371   productSpace_->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
00372 #ifdef TEUCHOS_DEBUG
00373   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00374 #endif
00375   if(
00376     rng.lbound() + rng.size()
00377     <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
00378     )
00379   {
00380     // This involves only one sub-vector so just return it.
00381     const_cast<const VectorBase<Scalar>*>(
00382       &*vecs_[kth_vector_space].getConstObj()
00383       )->acquireDetachedView( rng - kth_global_offset, sub_vec );
00384     sub_vec->setGlobalOffset( sub_vec->globalOffset() + kth_global_offset );
00385   }
00386   else {
00387     // Just let the default implementation handle this.  ToDo: In the future
00388     // we could manually construct an explicit sub-vector that spanned
00389     // two or more constituent vectors but this would be a lot of work.
00390     // However, this would require the use of temporary memory but
00391     // so what.
00392     VectorDefaultBase<Scalar>::acquireDetachedView(rng_in,sub_vec);
00393   }
00394 }
00395 
00396 template <class Scalar>
00397 void DefaultProductVector<Scalar>::releaseDetachedView(
00398   RTOpPack::ConstSubVectorView<Scalar>* sub_vec
00399   ) const
00400 {
00401   if( sub_vec->values() == NULL ) return;
00402   int    kth_vector_space  = -1;
00403   Index  kth_global_offset = 0;
00404   productSpace_->getVecSpcPoss(sub_vec->globalOffset(),&kth_vector_space,&kth_global_offset);
00405 #ifdef TEUCHOS_DEBUG
00406   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00407 #endif
00408   if(
00409     sub_vec->globalOffset() + sub_vec->subDim()
00410     <= kth_global_offset +  vecs_[kth_vector_space].getConstObj()->space()->dim()
00411     )
00412   {
00413     // This sub_vec was extracted from a single constituent vector
00414     sub_vec->setGlobalOffset( sub_vec->globalOffset() - kth_global_offset );
00415     vecs_[kth_vector_space].getConstObj()->releaseDetachedView(sub_vec);
00416   }
00417   else {
00418     // This sub_vec was created by the default implementation!
00419     VectorDefaultBase<Scalar>::releaseDetachedView(sub_vec);
00420   }
00421 }
00422 
00423 template <class Scalar>
00424 void DefaultProductVector<Scalar>::acquireDetachedView(
00425   const Range1D& rng_in, RTOpPack::SubVectorView<Scalar>* sub_vec
00426   )
00427 {
00428   const Range1D
00429     rng = rng_in.full_range() ? Range1D(0,productSpace_->dim()-1) : rng_in;
00430   int    kth_vector_space  = -1;
00431   Index  kth_global_offset = 0;
00432   productSpace_->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
00433 #ifdef TEUCHOS_DEBUG
00434   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00435 #endif
00436   if(
00437     rng.lbound() + rng.size()
00438     <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
00439     )
00440   {
00441     // This involves only one sub-vector so just return it.
00442     vecs_[kth_vector_space].getConstObj()->acquireDetachedView(
00443       rng - kth_global_offset, sub_vec
00444       );
00445     sub_vec->setGlobalOffset( sub_vec->globalOffset() + kth_global_offset );
00446   }
00447   else {
00448     // Just let the default implementation handle this.  ToDo: In the future
00449     // we could manually construct an explicit sub-vector that spanned
00450     // two or more constituent vectors but this would be a lot of work.
00451     // However, this would require the use of temporary memory but
00452     // so what.
00453     VectorDefaultBase<Scalar>::acquireDetachedView(rng_in,sub_vec);
00454   }
00455 }
00456 
00457 template <class Scalar>
00458 void DefaultProductVector<Scalar>::commitDetachedView(
00459   RTOpPack::SubVectorView<Scalar>* sub_vec
00460   )
00461 {
00462   if( sub_vec->values() == NULL ) return;
00463   int    kth_vector_space  = -1;
00464   Index  kth_global_offset = 0;
00465   productSpace_->getVecSpcPoss(sub_vec->globalOffset(),&kth_vector_space,&kth_global_offset);
00466 #ifdef TEUCHOS_DEBUG
00467   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00468 #endif
00469   if(
00470     sub_vec->globalOffset() + sub_vec->subDim()
00471     <= kth_global_offset +  vecs_[kth_vector_space].getConstObj()->space()->dim()
00472     )
00473   {
00474     // This sub_vec was extracted from a single constituent vector
00475     sub_vec->setGlobalOffset( sub_vec->globalOffset() - kth_global_offset );
00476     vecs_[kth_vector_space].getNonconstObj()->commitDetachedView(sub_vec);
00477   }
00478   else {
00479     // This sub_vec was created by the default implementation!
00480     VectorDefaultBase<Scalar>::commitDetachedView(sub_vec);
00481   }
00482 }
00483 
00484 template <class Scalar>
00485 void DefaultProductVector<Scalar>::setSubVector(
00486   const RTOpPack::SparseSubVectorT<Scalar>& sub_vec
00487   )
00488 {
00489   int    kth_vector_space  = -1;
00490   Index  kth_global_offset = 0;
00491   productSpace_->getVecSpcPoss(sub_vec.globalOffset(),&kth_vector_space,&kth_global_offset);
00492 #ifdef TEUCHOS_DEBUG
00493   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00494 #endif
00495   if(
00496     sub_vec.globalOffset() + sub_vec.subDim()
00497     <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
00498     )
00499   {
00500     // This sub-vector fits into a single constituent vector
00501     RTOpPack::SparseSubVectorT<Scalar> sub_vec_g = sub_vec;
00502     sub_vec_g.setGlobalOffset( sub_vec_g.globalOffset() - kth_global_offset );
00503     vecs_[kth_vector_space].getNonconstObj()->setSubVector(sub_vec_g);
00504   }
00505   else {
00506     // Let the default implementation take care of this.  ToDo: In the future
00507     // it would be possible to manually set the relevant constituent
00508     // vectors with no temp memory allocations.
00509     VectorDefaultBase<Scalar>::setSubVector(sub_vec);
00510   }
00511 }
00512 
00513 } // namespace Thyra
00514 
00515 #endif // THYRA_PRODUCT_VECTOR_HPP

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