Thyra Version of the Day
Thyra_DefaultSpmdMultiVector_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_SPMD_MULTI_VECTOR_DEF_HPP
00043 #define THYRA_DEFAULT_SPMD_MULTI_VECTOR_DEF_HPP
00044 
00045 // Define to make some verbose output
00046 //#define THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
00047 
00048 #include "Thyra_DefaultSpmdMultiVector_decl.hpp"
00049 #include "Thyra_SpmdMultiVectorDefaultBase.hpp"
00050 #include "Thyra_VectorSpaceFactoryBase.hpp"
00051 #include "Thyra_DefaultSpmdVector.hpp"
00052 #include "Teuchos_Assert.hpp"
00053 
00054 
00055 namespace Thyra {
00056 
00057 
00058 //
00059 // Simple utility class to copy back multi-vector entries
00060 //
00061 // ToDo: Refactor above to not use raw pointers!
00062 //
00063 template<class Scalar>
00064 class CopyBackSpmdMultiVectorEntries {
00065 public:
00066   CopyBackSpmdMultiVectorEntries(
00067     const ArrayView<const int> &cols,
00068     const ArrayRCP<const Scalar> &localValuesView, const Ordinal localSubDim,
00069     const ArrayRCP<Scalar> &localValues, const Ordinal leadingDim
00070     )
00071     : cols_(cols), localValuesView_(localValuesView), localSubDim_(localSubDim),
00072       localValues_(localValues), leadingDim_(leadingDim)
00073     {}
00074   ~CopyBackSpmdMultiVectorEntries()
00075     {
00076       typedef typename ArrayRCP<const Scalar>::const_iterator const_itr_t;
00077       typedef typename ArrayRCP<Scalar>::iterator itr_t;
00078       // Copy from contiguous storage column by column
00079       if (localValues_.strong_count()) {
00080         const int numCols = cols_.size();
00081         const const_itr_t lvv = localValuesView_.begin();
00082         const itr_t lv = localValues_.begin();
00083         for (int k = 0; k < numCols; ++k) {
00084           const int col_k = cols_[k];
00085           const const_itr_t lvv_k = lvv + localSubDim_*k;
00086           const itr_t lv_k = lv + leadingDim_*col_k;
00087           std::copy( lvv_k, lvv_k + localSubDim_, lv_k );
00088         }
00089       }
00090 #ifdef THYRA_DEBUG
00091       else {
00092         ++DefaultSpmdMultiVector<Scalar>::numSkipCopyBack;
00093       }
00094 #endif // THYRA_DEBUG
00095     }
00096 private:
00097   Array<int> cols_;
00098   ArrayRCP<const Scalar> localValuesView_;
00099   Ordinal localSubDim_;
00100   ArrayRCP<Scalar> localValues_;
00101   Ordinal leadingDim_;
00102   // Not defined and not to be called
00103   CopyBackSpmdMultiVectorEntries();
00104   CopyBackSpmdMultiVectorEntries(const CopyBackSpmdMultiVectorEntries&);
00105   CopyBackSpmdMultiVectorEntries& operator=(const CopyBackSpmdMultiVectorEntries&);
00106 };
00107 
00108 
00109 template<class Scalar>
00110 RCP<CopyBackSpmdMultiVectorEntries<Scalar> >
00111 copyBackSpmdMultiVectorEntries(
00112   const ArrayView<const int> &cols,
00113   const ArrayRCP<const Scalar> &localValuesView, const Ordinal localSubDim,
00114   const ArrayRCP<Scalar> &localValues, const Ordinal leadingDim
00115   )
00116 {
00117   return Teuchos::rcp(
00118     new CopyBackSpmdMultiVectorEntries<Scalar>(
00119       cols, localValuesView, localSubDim, localValues, leadingDim
00120       )
00121     );
00122 }
00123 
00124 
00125 //
00126 // DefaultSpmdMultiVector
00127 //
00128 
00129 
00130 #ifdef THYRA_DEBUG
00131 template<class Scalar>
00132 int DefaultSpmdMultiVector<Scalar>::numSkipCopyBack(0);
00133 #endif
00134 
00135 
00136 // Constructors/initializers/accessors
00137 
00138 
00139 template<class Scalar>
00140 DefaultSpmdMultiVector<Scalar>::DefaultSpmdMultiVector()
00141   :leadingDim_(0)
00142 {}
00143 
00144 
00145 template<class Scalar>
00146 DefaultSpmdMultiVector<Scalar>::DefaultSpmdMultiVector(
00147   const RCP<const SpmdVectorSpaceBase<Scalar> > &spmdRangeSpace
00148   ,const RCP<const ScalarProdVectorSpaceBase<Scalar> > &domainSpace
00149   )
00150 {
00151   initialize(spmdRangeSpace,domainSpace);
00152 }
00153 
00154 
00155 template<class Scalar>
00156 DefaultSpmdMultiVector<Scalar>::DefaultSpmdMultiVector(
00157   const RCP<const SpmdVectorSpaceBase<Scalar> > &spmdRangeSpace
00158   ,const RCP<const ScalarProdVectorSpaceBase<Scalar> > &domainSpace
00159   ,const ArrayRCP<Scalar> &localValues
00160   ,const Ordinal leadingDim
00161   )
00162 {
00163   initialize(spmdRangeSpace,domainSpace,localValues,leadingDim);
00164 }
00165 
00166 
00167 template<class Scalar>
00168 void DefaultSpmdMultiVector<Scalar>::initialize(
00169   const RCP<const SpmdVectorSpaceBase<Scalar> > &spmdRangeSpace
00170   ,const RCP<const ScalarProdVectorSpaceBase<Scalar> > &domainSpace
00171   )
00172 {
00173   const Ordinal localSubDim = spmdRangeSpace->localSubDim();
00174   ArrayRCP<Scalar> values;
00175   if (localSubDim)
00176     values = Teuchos::arcp<Scalar>(localSubDim * domainSpace->dim());
00177   initialize(spmdRangeSpace, domainSpace, values, localSubDim);
00178 }
00179 
00180 
00181 template<class Scalar>
00182 void DefaultSpmdMultiVector<Scalar>::initialize(
00183   const RCP<const SpmdVectorSpaceBase<Scalar> > &spmdRangeSpace,
00184   const RCP<const ScalarProdVectorSpaceBase<Scalar> > &domainSpace,
00185   const ArrayRCP<Scalar> &localValues,
00186   const Ordinal leadingDim_in
00187   )
00188 {
00189   const Ordinal localSubDim = spmdRangeSpace->localSubDim();
00190   const Ordinal leadingDim = (leadingDim_in >= 0 ? leadingDim_in : localSubDim);
00191 #ifdef TEUCHOS_DEBUG
00192   TEUCHOS_ASSERT(!is_null(spmdRangeSpace));
00193   TEUCHOS_ASSERT(!is_null(domainSpace));
00194   if (spmdRangeSpace->dim() && localSubDim) {
00195     TEUCHOS_ASSERT(nonnull(localValues));
00196   }
00197   TEUCHOS_ASSERT_INEQUALITY(leadingDim, >=, spmdRangeSpace->localSubDim());
00198 #endif
00199   spmdRangeSpace_ = spmdRangeSpace;
00200   domainSpace_ = domainSpace;
00201   localValues_ = localValues;
00202   leadingDim_ = leadingDim;
00203   this->updateSpmdSpace();
00204 }
00205 
00206 
00207 template<class Scalar>
00208 void DefaultSpmdMultiVector<Scalar>::uninitialize(
00209   RCP<const SpmdVectorSpaceBase<Scalar> > *spmdRangeSpace
00210   ,RCP<const ScalarProdVectorSpaceBase<Scalar> > *domainSpace
00211   ,ArrayRCP<Scalar> *localValues
00212   ,Ordinal *leadingDim
00213   )
00214 {
00215   if(spmdRangeSpace) *spmdRangeSpace = spmdRangeSpace_;
00216   if(domainSpace) *domainSpace = domainSpace_;
00217   if(localValues) *localValues = localValues_;
00218   if(leadingDim) *leadingDim = leadingDim_;
00219 
00220   spmdRangeSpace_ = Teuchos::null;
00221   domainSpace_ = Teuchos::null;
00222   localValues_ = Teuchos::null;
00223   leadingDim_ = 0;
00224 
00225   this->updateSpmdSpace();
00226 }
00227 
00228 
00229 template<class Scalar>
00230 RCP< const ScalarProdVectorSpaceBase<Scalar> >
00231 DefaultSpmdMultiVector<Scalar>::domainScalarProdVecSpc() const
00232 {
00233 #ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
00234   std::cerr << "\nSpmdMultiVectorStd<Scalar>::domainScalarProdVecSpc() const called!\n";
00235 #endif
00236   return domainSpace_;
00237 }
00238 
00239 
00240 // Overridden protected functions from MultiVectorBase
00241 
00242 
00243 template<class Scalar>
00244 RCP<VectorBase<Scalar> >
00245 DefaultSpmdMultiVector<Scalar>::nonconstColImpl(Ordinal j)
00246 {
00247 #ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
00248   std::cerr << "\nSpmdMultiVectorStd<Scalar>::col() called!\n";
00249 #endif
00250 #ifdef TEUCHOS_DEBUG
00251   TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= j && j < this->domain()->dim() ) );
00252 #endif
00253   return Teuchos::rcp(
00254     new DefaultSpmdVector<Scalar>(
00255       spmdRangeSpace_,
00256       localValues_.persistingView(j*leadingDim_,spmdRangeSpace_->localSubDim()),
00257       1
00258       )
00259     );
00260   //return Teuchos::rcp(new DefaultVectorMultiVector<Scalar>(subView(Range1D(j,j))));
00261 }
00262 
00263 
00264 template<class Scalar>
00265 RCP<const MultiVectorBase<Scalar> >
00266 DefaultSpmdMultiVector<Scalar>::contigSubViewImpl(
00267   const Range1D& col_rng_in
00268   ) const
00269 {
00270   const Range1D colRng = this->validateColRange(col_rng_in);
00271   return Teuchos::rcp(
00272     new DefaultSpmdMultiVector<Scalar>(
00273       spmdRangeSpace_,
00274       Teuchos::rcp_dynamic_cast<const ScalarProdVectorSpaceBase<Scalar> >(
00275         spmdRangeSpace_->smallVecSpcFcty()->createVecSpc(colRng.size())
00276         ,true
00277         ),
00278       localValues_.persistingView(colRng.lbound()*leadingDim_,colRng.size()*spmdRangeSpace_->localSubDim()),
00279       leadingDim_
00280       )
00281     );
00282 }
00283 
00284 
00285 template<class Scalar>
00286 RCP<MultiVectorBase<Scalar> >
00287 DefaultSpmdMultiVector<Scalar>::nonconstContigSubViewImpl(
00288   const Range1D& col_rng_in
00289   )
00290 {
00291   return Teuchos::rcp_const_cast<MultiVectorBase<Scalar> >(
00292     this->contigSubViewImpl(col_rng_in));
00293   // Have the nonconst version call the const version.  Note that in this case
00294   // we just need to take the const off of the returned MultiVectorBase object
00295   // because the localValues is already handled as nonconst.  This is the
00296   // perfect instance where the advice in Item 3 in "Effective C++ 3rd
00297   // edition" where Scott Meyers recommends having the nonconst version call
00298   // the const version.
00299 }
00300 
00301 
00302 template<class Scalar>
00303 RCP<const MultiVectorBase<Scalar> >
00304 DefaultSpmdMultiVector<Scalar>::nonContigSubViewImpl(
00305   const ArrayView<const int> &cols
00306   ) const
00307 {
00308   THYRA_DEBUG_ASSERT_MV_COLS("nonContigSubViewImpl(cols)", cols);
00309   const int numCols = cols.size();
00310   const ArrayRCP<Scalar> localValuesView = createContiguousCopy(cols);
00311   return defaultSpmdMultiVector<Scalar>(
00312     spmdRangeSpace_,
00313     createSmallScalarProdVectorSpaceBase<Scalar>(spmdRangeSpace_, numCols),
00314     localValuesView
00315     );
00316 }
00317 
00318 
00319 template<class Scalar>
00320 RCP<MultiVectorBase<Scalar> >
00321 DefaultSpmdMultiVector<Scalar>::nonconstNonContigSubViewImpl(
00322   const ArrayView<const int> &cols )
00323 {
00324   THYRA_DEBUG_ASSERT_MV_COLS("nonContigSubViewImpl(cols)", cols);
00325   const int numCols = cols.size();
00326   const ArrayRCP<Scalar> localValuesView = createContiguousCopy(cols);
00327   const Ordinal localSubDim = spmdRangeSpace_->localSubDim();
00328   RCP<CopyBackSpmdMultiVectorEntries<Scalar> > copyBackView =
00329     copyBackSpmdMultiVectorEntries<Scalar>(cols, localValuesView.getConst(),
00330       localSubDim, localValues_.create_weak(), leadingDim_);
00331   return Teuchos::rcpWithEmbeddedObjPreDestroy(
00332     new DefaultSpmdMultiVector<Scalar>(
00333       spmdRangeSpace_,
00334       createSmallScalarProdVectorSpaceBase<Scalar>(spmdRangeSpace_, numCols),
00335       localValuesView),
00336     copyBackView
00337     );
00338 }
00339 
00340 
00341 // Overridden protected members from SpmdMultiVectorBase
00342 
00343 
00344 template<class Scalar>
00345 RCP<const SpmdVectorSpaceBase<Scalar> >
00346 DefaultSpmdMultiVector<Scalar>::spmdSpaceImpl() const
00347 {
00348 #ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
00349   std::cerr << "\nSpmdMultiVectorStd<Scalar>::spmdSpace() const called!\n";
00350 #endif
00351   return spmdRangeSpace_;
00352 }
00353 
00354 
00355 template<class Scalar>
00356 void DefaultSpmdMultiVector<Scalar>::getNonconstLocalMultiVectorDataImpl(
00357   const Ptr<ArrayRCP<Scalar> > &localValues, const Ptr<Ordinal> &leadingDim
00358   )
00359 {
00360 #ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
00361   std::cerr << "\nSpmdMultiVectorStd<Scalar>::getLocalMultiVectorDataImpl() called!\n";
00362 #endif
00363   *localValues = localValues_;
00364   *leadingDim = leadingDim_;
00365 }
00366 
00367 
00368 template<class Scalar>
00369 void DefaultSpmdMultiVector<Scalar>::getLocalMultiVectorDataImpl(
00370   const Ptr<ArrayRCP<const Scalar> > &localValues, const Ptr<Ordinal> &leadingDim
00371   ) const
00372 {
00373 #ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
00374   std::cerr << "\nSpmdMultiVectorStd<Scalar>::getLocalData() called!\n";
00375 #endif
00376   *localValues = localValues_;
00377   *leadingDim = leadingDim_;
00378 }
00379 
00380 
00381 // private
00382 
00383 
00384 template<class Scalar>
00385 ArrayRCP<Scalar>
00386 DefaultSpmdMultiVector<Scalar>::createContiguousCopy(
00387   const ArrayView<const int> &cols ) const
00388 {
00389   typedef typename ArrayRCP<Scalar>::const_iterator const_itr_t;
00390   typedef typename ArrayRCP<Scalar>::iterator itr_t;
00391   const int numCols = cols.size();
00392   const Ordinal localSubDim = spmdRangeSpace_->localSubDim();
00393   ArrayRCP<Scalar> localValuesView = Teuchos::arcp<Scalar>(numCols*localSubDim);
00394   // Copy to contiguous storage column by column
00395   const const_itr_t lv = localValues_.begin();
00396   const itr_t lvv = localValuesView.begin();
00397   for (int k = 0; k < numCols; ++k) {
00398     const int col_k = cols[k];
00399     const const_itr_t lv_k = lv + leadingDim_*col_k;
00400     const itr_t lvv_k = lvv + localSubDim*k;
00401     std::copy(lv_k, lv_k+localSubDim, lvv_k);
00402   }
00403   return localValuesView;
00404 }
00405 
00406 
00407 } // end namespace Thyra
00408 
00409 
00410 #endif // THYRA_DEFAULT_SPMD_MULTI_VECTOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines