Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_TpetraVectorSpace_def.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 
43 #ifndef THYRA_TPETRA_VECTOR_SPACE_HPP
44 #define THYRA_TPETRA_VECTOR_SPACE_HPP
45 
46 
47 #include "Thyra_TpetraVectorSpace_decl.hpp"
48 #include "Thyra_TpetraThyraWrappers.hpp"
49 #include "Thyra_TpetraVector.hpp"
50 #include "Thyra_TpetraMultiVector.hpp"
51 #include "Thyra_TpetraEuclideanScalarProd.hpp"
52 #include "Tpetra_Details_StaticView.hpp"
53 
54 namespace Thyra {
55 
56 
57 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
58 RCP<TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
60 {
61  const RCP<this_t> vs(new this_t);
62  vs->weakSelfPtr_ = vs.create_weak();
63  return vs;
64 }
65 
66 
67 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
69  const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > &tpetraMap
70  )
71 {
72  comm_ = convertTpetraToThyraComm(tpetraMap->getComm());
73  tpetraMap_ = tpetraMap;
74  this->updateState(tpetraMap->getGlobalNumElements(),
75  !tpetraMap->isDistributed());
76  this->setScalarProd(tpetraEuclideanScalarProd<Scalar,LocalOrdinal,GlobalOrdinal,Node>());
77 }
78 
79 
80 // Utility functions
81 
82 
83 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
87 {
88  return tpetraVectorSpace<Scalar>(
89  Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
90  size, tpetraMap_->getComm() ) );
91 }
92 
93 
94 // Overridden from VectorSpace
95 
96 
97 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
100 {
101  return tpetraVector<Scalar>(
102  weakSelfPtr_.create_strong().getConst(),
103  Teuchos::rcp(
105  )
106  );
107 }
108 
109 
110 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
113 {
114  return tpetraMultiVector<Scalar>(
115  weakSelfPtr_.create_strong().getConst(),
116  this->createLocallyReplicatedVectorSpace(numMembers),
117  Teuchos::rcp(
119  tpetraMap_, numMembers, false)
120  )
121  );
122 }
123 
124 
125 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
126 class CopyTpetraMultiVectorViewBack {
127 public:
128  CopyTpetraMultiVectorViewBack( RCP<MultiVectorBase<Scalar> > mv, const RTOpPack::SubMultiVectorView<Scalar> &raw_mv )
129  :mv_(mv), raw_mv_(raw_mv)
130  {
132  bool inUse = Teuchos::get_extra_data<bool>(tmv,"inUse");
134  std::runtime_error,
135  "Cannot use the cached vector simultaneously more than once.");
136  inUse = true;
137  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tmv), Teuchos::POST_DESTROY, false);
138  }
139  ~CopyTpetraMultiVectorViewBack()
140  {
142  mv_->acquireDetachedView(Range1D(),Range1D(),&smv);
143  RTOpPack::assign_entries<Scalar>( Teuchos::outArg(raw_mv_), smv );
144  mv_->releaseDetachedView(&smv);
145  bool inUse = false;
146  RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmv = Teuchos::rcp_dynamic_cast<TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(mv_,true)->getTpetraMultiVector();
147  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tmv), Teuchos::POST_DESTROY, false);
148  }
149 private:
150  RCP<MultiVectorBase<Scalar> > mv_;
152 };
153 
154 
155 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
156 RCP< MultiVectorBase<Scalar> >
159  const bool initialize) const
160 {
161 #ifdef TEUCHOS_DEBUG
162  TEUCHOS_TEST_FOR_EXCEPT( raw_mv.subDim() != this->dim() );
163 #endif
164 
165  // Create a multi-vector
167  if (!tpetraMap_->isDistributed()) {
168 
169  if (tpetraMV_.is_null() || (tpetraMV_->getNumVectors() != size_t (raw_mv.numSubCols()))) {
170  if (!tpetraMV_.is_null())
171  // The MV is already allocated. If we are still using it, then very bad things can happen.
172  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::get_extra_data<bool>(tpetraMV_,"inUse"),
173  std::runtime_error,
174  "Cannot use the cached vector simultaneously more than once.");
177  auto dv = ::Tpetra::Details::getStatic2dDualView<IST, DT> (tpetraMap_->getGlobalNumElements(), raw_mv.numSubCols());
179  bool inUse = false;
180  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tpetraMV_));
181  }
182 
183  if (tpetraDomainSpace_.is_null() || raw_mv.numSubCols() != tpetraDomainSpace_->localSubDim())
184  tpetraDomainSpace_ = tpetraVectorSpace<Scalar>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(raw_mv.numSubCols(), tpetraMap_->getComm()));
185 
186  mv = tpetraMultiVector<Scalar>(weakSelfPtr_.create_strong().getConst(), tpetraDomainSpace_, tpetraMV_);
187  } else {
188  mv = this->createMembers(raw_mv.numSubCols());
189  bool inUse = false;
191  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tmv));
192  }
193  if (initialize) {
194  // Copy initial values in raw_mv into multi-vector
196  mv->acquireDetachedView(Range1D(),Range1D(),&smv);
197  RTOpPack::assign_entries<Scalar>(
198  Ptr<const RTOpPack::SubMultiVectorView<Scalar> >(Teuchos::outArg(smv)),
199  raw_mv
200  );
201  mv->commitDetachedView(&smv);
202  }
203  // Setup smart pointer to multi-vector to copy view back out just before multi-vector is destroyed
204  Teuchos::set_extra_data(
205  // We create a duplicate of the RCP, otherwise the ref count does not go to zero.
206  Teuchos::rcp(new CopyTpetraMultiVectorViewBack<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Teuchos::rcpFromRef(*mv),raw_mv)),
207  "CopyTpetraMultiVectorViewBack",
208  Teuchos::outArg(mv),
209  Teuchos::PRE_DESTROY
210  );
211  return mv;
212 }
213 
214 
215 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
218  const RTOpPack::ConstSubMultiVectorView<Scalar> &raw_mv ) const
219 {
220 #ifdef TEUCHOS_DEBUG
221  TEUCHOS_TEST_FOR_EXCEPT( raw_mv.subDim() != this->dim() );
222 #endif
223  // Create a multi-vector
225  if (!tpetraMap_->isDistributed()) {
226  if (tpetraMV_.is_null() || (tpetraMV_->getNumVectors() != size_t (raw_mv.numSubCols()))) {
227  if (!tpetraMV_.is_null())
228  // The MV is already allocated. If we are still using it, then very bad things can happen.
229  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::get_extra_data<bool>(tpetraMV_,"inUse"),
230  std::runtime_error,
231  "Cannot use the cached vector simultaneously more than once.");
234  auto dv = ::Tpetra::Details::getStatic2dDualView<IST, DT> (tpetraMap_->getGlobalNumElements(), raw_mv.numSubCols());
236  bool inUse = false;
237  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tpetraMV_));
238  }
239 
240  if (tpetraDomainSpace_.is_null() || raw_mv.numSubCols() != tpetraDomainSpace_->localSubDim())
241  tpetraDomainSpace_ = tpetraVectorSpace<Scalar>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(raw_mv.numSubCols(), tpetraMap_->getComm()));
242 
243  mv = tpetraMultiVector<Scalar>(weakSelfPtr_.create_strong().getConst(), tpetraDomainSpace_, tpetraMV_);
244  } else {
245  mv = this->createMembers(raw_mv.numSubCols());
246  bool inUse = false;
248  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tmv));
249  }
250  // Copy values in raw_mv into multi-vector
252  mv->acquireDetachedView(Range1D(),Range1D(),&smv);
253  RTOpPack::assign_entries<Scalar>(
254  Ptr<const RTOpPack::SubMultiVectorView<Scalar> >(Teuchos::outArg(smv)),
255  raw_mv );
256  mv->commitDetachedView(&smv);
257  return mv;
258 }
259 
260 
261 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
263  const Range1D& rng_in, const EViewType viewType, const EStrideType strideType
264  ) const
265 {
266  const Range1D rng = full_range(rng_in,0,this->dim()-1);
267  const Ordinal l_localOffset = this->localOffset();
268 
269  const Ordinal myLocalSubDim = tpetraMap_.is_null () ?
270  static_cast<Ordinal> (0) : tpetraMap_->getLocalNumElements ();
271 
272  return ( l_localOffset<=rng.lbound() && rng.ubound()<l_localOffset+myLocalSubDim );
273 }
274 
275 
276 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
279 {
280  return tpetraVectorSpace<Scalar>(tpetraMap_);
281 }
282 
283 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
286 {
287  return tpetraMap_;
288 }
289 
290 // Overridden from SpmdVectorSpaceDefaultBase
291 
292 
293 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
296 {
297  return comm_;
298 }
299 
300 
301 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
303 {
304  return tpetraMap_.is_null () ? static_cast<Ordinal> (0) :
305  static_cast<Ordinal> (tpetraMap_->getLocalNumElements ());
306 }
307 
308 // private
309 
310 
311 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
313 {
314  // The base classes should automatically default initialize to a safe
315  // uninitialized state.
316 }
317 
318 
319 } // end namespace Thyra
320 
321 
322 #endif // THYRA_TPETRA_VECTOR_SPACE_HPP
Concrete implementation of Thyra::MultiVector in terms of Tpetra::MultiVector.
Concrete implementation of an SPMD vector space for Tpetra.
RCP< const VectorSpaceBase< Scalar > > clone() const
RCP< T > create_weak() const
typename map_type::device_type device_type
RCP< TpetraVectorSpace< Scalar, LocalOrdinal, GlobalOrdinal, Node > > createLocallyReplicatedVectorSpace(int size) const
Create Tpetra locally replicated vector space.
static RCP< TpetraVectorSpace< Scalar, LocalOrdinal, GlobalOrdinal, Node > > create()
Create with weak ownership to self.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void initialize(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &tpetraMap)
Initialize a serial space.
EViewType
Determines if a view is a direct view of data or a detached copy of data.
Ordinal ubound() const
typename Kokkos::ArithTraits< Scalar >::val_type impl_scalar_type
RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getTpetraMap() const
Get the embedded Tpetra::Map.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Interface for a collection of column vectors called a multi-vector.
bool hasInCoreView(const Range1D &rng, const EViewType viewType, const EStrideType strideType) const
Returns true if all the elements in rng are in this process.
RCP< const Teuchos::Comm< Ordinal > > convertTpetraToThyraComm(const RCP< const Teuchos::Comm< int > > &tpetraComm)
Given an Tpetra Teuchos::Comm&lt;int&gt; object, return an equivalent Teuchos::Comm&lt;Ordinal&gt; object...
RCP< VectorBase< Scalar > > createMember() const
Ordinal lbound() const
RCP< MultiVectorBase< Scalar > > createMembers(int numMembers) const
EStrideType
Determine if data is unit stride or non-unit stride.
RCP< const Teuchos::Comm< Ordinal > > getComm() const
RCP< MultiVectorBase< Scalar > > createCachedMembersView(const RTOpPack::SubMultiVectorView< Scalar > &raw_mv, bool initialize=true) const
Create a (possibly) cached multi-vector member that is a non-const view of raw multi-vector data...
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Teuchos::Range1D Range1D