Thyra Version of the Day
Loading...
Searching...
No Matches
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
54namespace Thyra {
55
56
57template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
58RCP<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
67template <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// Overridden from VectorSpace
81
82
83template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
86{
87 return tpetraVector<Scalar>(
88 weakSelfPtr_.create_strong().getConst(),
90 new Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetraMap_, false)
91 )
92 );
93}
94
95
96template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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 ),
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
117template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
118class CopyTpetraMultiVectorViewBack {
119public:
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");
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 {
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 }
141private:
142 RCP<MultiVectorBase<Scalar> > mv_;
144};
145
146
147template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
148RCP< MultiVectorBase<Scalar> >
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
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
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
204template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
208{
209#ifdef TEUCHOS_DEBUG
210 TEUCHOS_TEST_FOR_EXCEPT( raw_mv.subDim() != this->dim() );
211#endif
212 // Create a multi-vector
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
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
250template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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_->getLocalNumElements ();
260
261 return ( l_localOffset<=rng.lbound() && rng.ubound()<l_localOffset+myLocalSubDim );
262}
263
264
265template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
268{
269 return tpetraVectorSpace<Scalar>(tpetraMap_);
270}
271
272template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
275{
276 return tpetraMap_;
277}
278
279// Overridden from SpmdVectorSpaceDefaultBase
280
281
282template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
285{
286 return comm_;
287}
288
289
290template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
292{
293 return tpetraMap_.is_null () ? static_cast<Ordinal> (0) :
294 static_cast<Ordinal> (tpetraMap_->getLocalNumElements ());
295}
296
297// private
298
299
300template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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
RCP< T > create_weak() const
Ordinal lbound() const
Ordinal ubound() const
Interface for a collection of column vectors called a multi-vector.
Concrete implementation of an SPMD vector space for Tpetra.
RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getTpetraMap() const
Get the embedded Tpetra::Map.
bool hasInCoreView(const Range1D &rng, const EViewType viewType, const EStrideType strideType) const
Returns true if all the elements in rng are in this process.
void initialize(const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &tpetraMap)
Initialize a serial space.
RCP< MultiVectorBase< Scalar > > createCachedMembersView(const RTOpPack::SubMultiVectorView< Scalar > &raw_mv) const
Create a (possibly) cached multi-vector member that is a non-const view of raw multi-vector data....
RCP< const VectorSpaceBase< Scalar > > clone() const
RCP< VectorBase< Scalar > > createMember() const
RCP< MultiVectorBase< Scalar > > createMembers(int numMembers) const
RCP< const Teuchos::Comm< Ordinal > > getComm() const
static RCP< TpetraVectorSpace< Scalar, LocalOrdinal, GlobalOrdinal, Node > > create()
Create with weak ownership to self.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
EStrideType
Determine if data is unit stride or non-unit stride.
EViewType
Determines if a view is a direct view of data or a detached copy of data.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Teuchos::Range1D Range1D
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)