Thyra  Version of the Day
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> >
59 TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::create()
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>
68 void TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::initialize(
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 // Overridden from VectorSpace
81 
82 
83 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84 RCP<VectorBase<Scalar> >
85 TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createMember() const
86 {
87  return tpetraVector<Scalar>(
88  weakSelfPtr_.create_strong().getConst(),
89  Teuchos::rcp(
90  new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetraMap_, false)
91  )
92  );
93 }
94 
95 
96 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
97 RCP< MultiVectorBase<Scalar> >
98 TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createMembers(int numMembers) const
99 {
100  return tpetraMultiVector<Scalar>(
101  weakSelfPtr_.create_strong().getConst(),
102  tpetraVectorSpace<Scalar>(
103  Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(
104  numMembers, tpetraMap_->getComm()
105  )
106  ),
107  Teuchos::rcp(
108  new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(
109  tpetraMap_, numMembers, false)
110  )
111  );
112  // ToDo: Create wrapper function to create locally replicated vector space
113  // and use it.
114 }
115 
116 
117 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
118 class CopyTpetraMultiVectorViewBack {
119 public:
120  CopyTpetraMultiVectorViewBack( RCP<MultiVectorBase<Scalar> > mv, const RTOpPack::SubMultiVectorView<Scalar> &raw_mv )
121  :mv_(mv), raw_mv_(raw_mv)
122  {
123  RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmv = Teuchos::rcp_dynamic_cast<TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(mv_,true)->getTpetraMultiVector();
124  bool inUse = Teuchos::get_extra_data<bool>(tmv,"inUse");
125  TEUCHOS_TEST_FOR_EXCEPTION(inUse,
126  std::runtime_error,
127  "Cannot use the cached vector simultaneously more than once.");
128  inUse = true;
129  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tmv), Teuchos::POST_DESTROY, false);
130  }
131  ~CopyTpetraMultiVectorViewBack()
132  {
133  RTOpPack::ConstSubMultiVectorView<Scalar> smv;
134  mv_->acquireDetachedView(Range1D(),Range1D(),&smv);
135  RTOpPack::assign_entries<Scalar>( Teuchos::outArg(raw_mv_), smv );
136  mv_->releaseDetachedView(&smv);
137  bool inUse = false;
138  RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmv = Teuchos::rcp_dynamic_cast<TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(mv_,true)->getTpetraMultiVector();
139  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tmv), Teuchos::POST_DESTROY, false);
140  }
141 private:
142  RCP<MultiVectorBase<Scalar> > mv_;
143  const RTOpPack::SubMultiVectorView<Scalar> raw_mv_;
144 };
145 
146 
147 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
148 RCP< MultiVectorBase<Scalar> >
149 TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createCachedMembersView(
150  const RTOpPack::SubMultiVectorView<Scalar> &raw_mv ) const
151 {
152 #ifdef TEUCHOS_DEBUG
153  TEUCHOS_TEST_FOR_EXCEPT( raw_mv.subDim() != this->dim() );
154 #endif
155 
156  // Create a multi-vector
157  RCP< MultiVectorBase<Scalar> > mv;
158  if (!tpetraMap_->isDistributed()) {
159 
160  if (tpetraMV_.is_null() || (tpetraMV_->getNumVectors() != size_t (raw_mv.numSubCols()))) {
161  if (!tpetraMV_.is_null())
162  // The MV is already allocated. If we are still using it, then very bad things can happen.
163  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::get_extra_data<bool>(tpetraMV_,"inUse"),
164  std::runtime_error,
165  "Cannot use the cached vector simultaneously more than once.");
166  using IST = typename Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::impl_scalar_type;
167  using DT = typename Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::device_type;
168  auto dv = ::Tpetra::Details::getStatic2dDualView<IST, DT> (tpetraMap_->getGlobalNumElements(), raw_mv.numSubCols());
169  tpetraMV_ = Teuchos::rcp(new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetraMap_, dv));
170  bool inUse = false;
171  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tpetraMV_));
172  }
173 
174  if (tpetraDomainSpace_.is_null() || raw_mv.numSubCols() != tpetraDomainSpace_->localSubDim())
175  tpetraDomainSpace_ = tpetraVectorSpace<Scalar>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(raw_mv.numSubCols(), tpetraMap_->getComm()));
176 
177  mv = tpetraMultiVector<Scalar>(weakSelfPtr_.create_strong().getConst(), tpetraDomainSpace_, tpetraMV_);
178  } else {
179  mv = this->createMembers(raw_mv.numSubCols());
180  bool inUse = false;
181  RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmv = Teuchos::rcp_dynamic_cast<TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(mv,true)->getTpetraMultiVector();
182  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tmv));
183  }
184  // Copy initial values in raw_mv into multi-vector
185  RTOpPack::SubMultiVectorView<Scalar> smv;
186  mv->acquireDetachedView(Range1D(),Range1D(),&smv);
187  RTOpPack::assign_entries<Scalar>(
188  Ptr<const RTOpPack::SubMultiVectorView<Scalar> >(Teuchos::outArg(smv)),
189  raw_mv
190  );
191  mv->commitDetachedView(&smv);
192  // Setup smart pointer to multi-vector to copy view back out just before multi-vector is destroyed
193  Teuchos::set_extra_data(
194  // We create a duplicate of the RCP, otherwise the ref count does not go to zero.
195  Teuchos::rcp(new CopyTpetraMultiVectorViewBack<Scalar,LocalOrdinal,GlobalOrdinal,Node>(Teuchos::rcpFromRef(*mv),raw_mv)),
196  "CopyTpetraMultiVectorViewBack",
197  Teuchos::outArg(mv),
198  Teuchos::PRE_DESTROY
199  );
200  return mv;
201 }
202 
203 
204 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
205 RCP<const MultiVectorBase<Scalar> >
206 TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createCachedMembersView(
207  const RTOpPack::ConstSubMultiVectorView<Scalar> &raw_mv ) const
208 {
209 #ifdef TEUCHOS_DEBUG
210  TEUCHOS_TEST_FOR_EXCEPT( raw_mv.subDim() != this->dim() );
211 #endif
212  // Create a multi-vector
213  RCP< MultiVectorBase<Scalar> > mv;
214  if (!tpetraMap_->isDistributed()) {
215  if (tpetraMV_.is_null() || (tpetraMV_->getNumVectors() != size_t (raw_mv.numSubCols()))) {
216  if (!tpetraMV_.is_null())
217  // The MV is already allocated. If we are still using it, then very bad things can happen.
218  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::get_extra_data<bool>(tpetraMV_,"inUse"),
219  std::runtime_error,
220  "Cannot use the cached vector simultaneously more than once.");
221  using IST = typename Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::impl_scalar_type;
222  using DT = typename Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::device_type;
223  auto dv = ::Tpetra::Details::getStatic2dDualView<IST, DT> (tpetraMap_->getGlobalNumElements(), raw_mv.numSubCols());
224  tpetraMV_ = Teuchos::rcp(new Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetraMap_, dv));
225  bool inUse = false;
226  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tpetraMV_));
227  }
228 
229  if (tpetraDomainSpace_.is_null() || raw_mv.numSubCols() != tpetraDomainSpace_->localSubDim())
230  tpetraDomainSpace_ = tpetraVectorSpace<Scalar>(Tpetra::createLocalMapWithNode<LocalOrdinal, GlobalOrdinal, Node>(raw_mv.numSubCols(), tpetraMap_->getComm()));
231 
232  mv = tpetraMultiVector<Scalar>(weakSelfPtr_.create_strong().getConst(), tpetraDomainSpace_, tpetraMV_);
233  } else {
234  mv = this->createMembers(raw_mv.numSubCols());
235  bool inUse = false;
236  RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmv = Teuchos::rcp_dynamic_cast<TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(mv,true)->getTpetraMultiVector();
237  Teuchos::set_extra_data(inUse,"inUse",Teuchos::outArg(tmv));
238  }
239  // Copy values in raw_mv into multi-vector
240  RTOpPack::SubMultiVectorView<Scalar> smv;
241  mv->acquireDetachedView(Range1D(),Range1D(),&smv);
242  RTOpPack::assign_entries<Scalar>(
243  Ptr<const RTOpPack::SubMultiVectorView<Scalar> >(Teuchos::outArg(smv)),
244  raw_mv );
245  mv->commitDetachedView(&smv);
246  return mv;
247 }
248 
249 
250 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
251 bool TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::hasInCoreView(
252  const Range1D& rng_in, const EViewType viewType, const EStrideType strideType
253  ) const
254 {
255  const Range1D rng = full_range(rng_in,0,this->dim()-1);
256  const Ordinal l_localOffset = this->localOffset();
257 
258  const Ordinal myLocalSubDim = tpetraMap_.is_null () ?
259  static_cast<Ordinal> (0) : tpetraMap_->getNodeNumElements ();
260 
261  return ( l_localOffset<=rng.lbound() && rng.ubound()<l_localOffset+myLocalSubDim );
262 }
263 
264 
265 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
266 RCP< const VectorSpaceBase<Scalar> >
267 TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::clone() const
268 {
269  return tpetraVectorSpace<Scalar>(tpetraMap_);
270 }
271 
272 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
273 RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >
274 TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getTpetraMap() const
275 {
276  return tpetraMap_;
277 }
278 
279 // Overridden from SpmdVectorSpaceDefaultBase
280 
281 
282 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
283 RCP<const Teuchos::Comm<Ordinal> >
284 TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getComm() const
285 {
286  return comm_;
287 }
288 
289 
290 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
291 Ordinal TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::localSubDim() const
292 {
293  return tpetraMap_.is_null () ? static_cast<Ordinal> (0) :
294  static_cast<Ordinal> (tpetraMap_->getNodeNumElements ());
295 }
296 
297 // private
298 
299 
300 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
301 TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TpetraVectorSpace()
302 {
303  // The base classes should automatically default initialize to a safe
304  // uninitialized state.
305 }
306 
307 
308 } // end namespace Thyra
309 
310 
311 #endif // THYRA_TPETRA_VECTOR_SPACE_HPP
Teuchos::Range1D::lbound
Ordinal lbound() const
Thyra::TpetraVectorSpace::getTpetraMap
RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getTpetraMap() const
Get the embedded Tpetra::Map.
Definition: Thyra_TpetraVectorSpace_def.hpp:274
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Teuchos::Range1D::ubound
Ordinal ubound() const
RTOpPack::ConstSubMultiVectorView::subDim
Ordinal subDim() const
Thyra::convertTpetraToThyraComm
RCP< const Teuchos::Comm< Ordinal > > convertTpetraToThyraComm(const RCP< const Teuchos::Comm< int > > &tpetraComm)
Given an Tpetra Teuchos::Comm<int> object, return an equivalent Teuchos::Comm<Ordinal> object.
Definition: Thyra_TpetraThyraWrappers.cpp:55
Thyra::TpetraVectorSpace::getComm
RCP< const Teuchos::Comm< Ordinal > > getComm() const
Definition: Thyra_TpetraVectorSpace_def.hpp:284
Thyra::EViewType
EViewType
Determines if a view is a direct view of data or a detached copy of data.
Definition: Thyra_OperatorVectorTypes.hpp:309
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
RTOpPack::SubMultiVectorView< Scalar >
Teuchos::RCP
Thyra::TpetraVectorSpace::localSubDim
Ordinal localSubDim() const
Definition: Thyra_TpetraVectorSpace_def.hpp:291
Thyra::Ordinal
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Definition: Thyra_OperatorVectorTypes.hpp:127
Thyra::TpetraVectorSpace::create
static RCP< TpetraVectorSpace< Scalar, LocalOrdinal, GlobalOrdinal, Node > > create()
Create with weak ownership to self.
Definition: Thyra_TpetraVectorSpace_def.hpp:59
Thyra::TpetraVectorSpace::hasInCoreView
bool hasInCoreView(const Range1D &rng, const EViewType viewType, const EStrideType strideType) const
Returns true if all the elements in rng are in this process.
Definition: Thyra_TpetraVectorSpace_def.hpp:251
Thyra::TpetraVectorSpace::initialize
void initialize(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &tpetraMap)
Initialize a serial space.
Definition: Thyra_TpetraVectorSpace_def.hpp:68
Thyra::MultiVectorBase
Interface for a collection of column vectors called a multi-vector.
Definition: Thyra_MultiVectorBase_decl.hpp:496
RTOpPack::ConstSubMultiVectorView::numSubCols
Ordinal numSubCols() const
Thyra::TpetraVectorSpace::createMembers
RCP< MultiVectorBase< Scalar > > createMembers(int numMembers) const
Definition: Thyra_TpetraVectorSpace_def.hpp:98
Thyra::EStrideType
EStrideType
Determine if data is unit stride or non-unit stride.
Definition: Thyra_OperatorVectorTypes.hpp:319
Teuchos::RCP::create_weak
RCP< T > create_weak() const
Thyra::TpetraVectorSpace
Concrete implementation of an SPMD vector space for Tpetra.
Definition: Thyra_TpetraVectorSpace_decl.hpp:61
Thyra::TpetraVectorSpace::createMember
RCP< VectorBase< Scalar > > createMember() const
Definition: Thyra_TpetraVectorSpace_def.hpp:85
RTOpPack::ConstSubMultiVectorView< Scalar >
Thyra::TpetraVectorSpace::clone
RCP< const VectorSpaceBase< Scalar > > clone() const
Definition: Thyra_TpetraVectorSpace_def.hpp:267
Teuchos::Range1D
Thyra::Range1D
Teuchos::Range1D Range1D
Definition: Thyra_OperatorVectorTypes.hpp:97
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Generated by doxygen 1.8.20