Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultProductVectorSpace_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#ifndef THYRA_DEFAULT_PRODUCT_VECTOR_SPACE_HPP
43#define THYRA_DEFAULT_PRODUCT_VECTOR_SPACE_HPP
44
45
46#include "Thyra_DefaultProductVectorSpace_decl.hpp"
47#include "Thyra_DefaultProductVector.hpp"
48#include "Thyra_DefaultProductMultiVector.hpp"
49#include "Thyra_ProductMultiVectorBase.hpp"
50#include "Teuchos_Workspace.hpp"
51#include "Teuchos_dyn_cast.hpp"
52
53
54namespace Thyra {
55
56
57// Constructors/initializers/accessors
58
59
60template<class Scalar>
62 : numBlocks_(-1), dim_(-1), isInCore_(false)
63{}
64
65
66template<class Scalar>
68 const ArrayView<const RCP<const VectorSpaceBase<Scalar> > > &vecSpaces_in
69 )
70 : numBlocks_(-1), dim_(-1)
71{
72 initialize(vecSpaces_in);
73}
74
75
76template<class Scalar>
78 const ArrayView<const RCP<const VectorSpaceBase<Scalar> > > &vecSpaces_in
79 )
80{
81
82 //
83 // Check preconditions and compute cached quantities
84 //
85 const int nBlocks = vecSpaces_in.size();
86#ifdef TEUCHOS_DEBUG
87 TEUCHOS_TEST_FOR_EXCEPT( nBlocks == 0 );
88#endif
89 bool overallHasInCoreView = true;
90 for (int k = 0; k < nBlocks; ++k) {
91#ifdef TEUCHOS_DEBUG
93 vecSpaces_in[k].get() == NULL, std::invalid_argument
94 ,"Error, the smart pointer vecSpaces["<<k<<"] can not be NULL!"
95 );
96#endif
97 if (!vecSpaces_in[k]->hasInCoreView()) overallHasInCoreView = false;
98 }
99
100 //
101 // Setup private data members (should not throw an exception from here)
102 //
103 numBlocks_ = nBlocks;
104 vecSpaces_ = Teuchos::rcp(new vecSpaces_t);
105 *vecSpaces_ = vecSpaces_in;
106 vecSpacesOffsets_ = Teuchos::rcp(new vecSpacesOffsets_t(nBlocks+1));
107 (*vecSpacesOffsets_)[0] = 0;
108 dim_ = 0;
109 for( int k = 1; k <= nBlocks; ++k ) {
110 const Ordinal dim_km1 = vecSpaces_in[k-1]->dim();
111 (*vecSpacesOffsets_)[k] = (*vecSpacesOffsets_)[k-1] + dim_km1;
112 dim_ += dim_km1;
113 }
114 isInCore_ = overallHasInCoreView;
115
116}
117
118
119template<class Scalar>
121 const ArrayView<RCP<const VectorSpaceBase<Scalar> > > &vecSpaces_in
122 )
123{
124 TEUCHOS_TEST_FOR_EXCEPT(!is_null(vecSpaces_in)); // ToDo: Implement!
125 vecSpaces_ = Teuchos::null;
126 vecSpacesOffsets_ = Teuchos::null;
127 numBlocks_ = -1;
128 dim_ = -1;
129 isInCore_ = false;
130}
131
132
133template<class Scalar>
135 Ordinal i, int* kth_vector_space, Ordinal* kth_global_offset
136 ) const
137{
138 // Validate the preconditions
139#ifdef TEUCHOS_DEBUG
141 !(0 <= i && i < this->dim()), std::out_of_range
142 ,"VectorSpaceBlocked::get_vector_space_position(...): Error, i = "
143 << i << " is not in range [0,"<<(this->dim()-1)<<"]"
144 );
145#endif
146 *kth_vector_space = 0;
147 *kth_global_offset = 0;
148 while( *kth_vector_space < numBlocks_ ) {
149 const Ordinal off_kp1 = (*vecSpacesOffsets_)[*kth_vector_space+1];
150 if( off_kp1 > i ) {
151 *kth_global_offset = (*vecSpacesOffsets_)[*kth_vector_space];
152 break;
153 }
154 ++(*kth_vector_space);
155 }
156 TEUCHOS_TEST_FOR_EXCEPT( !(*kth_vector_space < numBlocks_) );
157}
158
159
160// Overridden from DefaultProductVectorSpace
161
162
163template<class Scalar>
165{
166 return numBlocks_;
167}
168
169
170template<class Scalar>
173{
174 TEUCHOS_TEST_FOR_EXCEPT( k < 0 || numBlocks_ < k );
175 return (*vecSpaces_)[k];
176}
177
178
179// Overridden from VectorSpaceBase
180
181
182template<class Scalar>
184{
185 return dim_;
186}
187
188
189template<class Scalar>
191 const VectorSpaceBase<Scalar>& vecSpc ) const
192{
193
194 using Teuchos::ptrFromRef;
195 using Teuchos::ptr_dynamic_cast;
196
197 const int nBlocks = this->numBlocks();
198
199 // Check for product vector interface
201 ptr_dynamic_cast<const ProductVectorSpaceBase<Scalar> >(ptrFromRef(vecSpc));
202
203 if (nonnull(pvsb)) {
204 // Validate that constituent vector spaces are compatible
205 if( nBlocks != pvsb->numBlocks() )
206 return false;
207 for( int i = 0; i < nBlocks; ++i ) {
208 if( !this->getBlock(i)->isCompatible(*pvsb->getBlock(i)) )
209 return false;
210 }
211 return true;
212 }
213
214 // Check for a single vector single vector space
215 if (nBlocks == 1) {
216 return this->getBlock(0)->isCompatible(vecSpc);
217 }
218
219 // If we get here, the RHS is not a product vector space and/or this is not
220 // a single block VS so we can assume the spaces are *not* compatible!
221 return false;
222
223}
224
225
226template<class Scalar>
229{
230 return defaultProductVector<Scalar>(Teuchos::rcpFromRef(*this));
231}
232
233
234template<class Scalar>
236 const VectorBase<Scalar> &x_in,
237 const VectorBase<Scalar> &y_in
238 ) const
239{
240 const int nBlocks = this->numBlocks();
244#ifdef TEUCHOS_DEBUG
246 nBlocks!=x.productSpace()->numBlocks()
247 || nBlocks!=y.productSpace()->numBlocks()
248 );
249#endif
250 Scalar scalarProd_rtn = Teuchos::ScalarTraits<Scalar>::zero();
251 for( int k = 0; k < nBlocks; ++k )
252 scalarProd_rtn += (*vecSpaces_)[k]->scalarProd(
253 *x.getVectorBlock(k),*y.getVectorBlock(k)
254 );
255 return scalarProd_rtn;
256}
257
258
259template<class Scalar>
261 const MultiVectorBase<Scalar> &X_in,
262 const MultiVectorBase<Scalar> &Y_in,
263 const ArrayView<Scalar> &scalarProds_out
264 ) const
265{
266 using Teuchos::as;
267 using Teuchos::Workspace;
268 const VectorSpaceBase<Scalar> &domain = *X_in.domain();
269 const Ordinal m = domain.dim();
270#ifdef TEUCHOS_DEBUG
271 TEUCHOS_TEST_FOR_EXCEPT(is_null(scalarProds_out));
272 TEUCHOS_TEST_FOR_EXCEPT( !domain.isCompatible(*Y_in.domain()) );
273 TEUCHOS_ASSERT_EQUALITY( as<Ordinal>(scalarProds_out.size()),
274 as<Ordinal>(m) )
275#endif
276 if(m==1) {
277 scalarProds_out[0] = this->scalarProd(*X_in.col(0),*Y_in.col(0));
278 return;
279 // ToDo: Remove this if(...) block once we have a DefaultProductMultiVector implementation!
280 }
282 const int nBlocks = this->numBlocks();
286#ifdef TEUCHOS_DEBUG
287 TEUCHOS_TEST_FOR_EXCEPT( nBlocks!=X.productSpace()->numBlocks() || nBlocks!=Y.productSpace()->numBlocks() );
288#endif
289 Workspace<Scalar> _scalarProds_out(wss, m, false);
290 std::fill( scalarProds_out.begin(), scalarProds_out.end(),
292 for( int k = 0; k < nBlocks; ++k ) {
293 (*vecSpaces_)[k]->scalarProds(
294 *X.getMultiVectorBlock(k), *Y.getMultiVectorBlock(k), _scalarProds_out());
295 for( int j = 0; j < m; ++j )
296 scalarProds_out[j] += _scalarProds_out[j];
297 }
298}
299
300
301template<class Scalar>
302bool DefaultProductVectorSpace<Scalar>::hasInCoreView(const Range1D& rng_in, const EViewType viewType, const EStrideType strideType) const
303{
304 const Range1D rng = full_range(rng_in,0,dim_-1);
305 // First see if rng fits in a single constituent vector
306 int kth_vector_space = -1;
307 Ordinal kth_global_offset = 0;
308 this->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
309#ifdef TEUCHOS_DEBUG
310 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
311#endif
312 if( rng.lbound() + rng.size() <= kth_global_offset + (*vecSpaces_)[kth_vector_space]->dim() ) {
313 return (*vecSpaces_)[kth_vector_space]->hasInCoreView(rng_in-kth_global_offset,viewType,strideType);
314 }
315 // If we get here, rng does not fit in a single constituent vector which
316 // also means that numBlocks_ > 1 must also be true!
317 //
318 // Next, if the client is asking for a direct view then we have to return
319 // false since this range spans more than one constituent vector.
320 if( viewType == VIEW_TYPE_DIRECT )
321 return false;
322 // If we get here then hasDirectView==false and therefore we are allowed to
323 // create a copy. Therefore, if all of the constituent vectors are "in
324 // core" then we can return true.
325 if(isInCore_)
326 return true;
327 // Finally, loop through all of the constituent vectors spaned by rng and
328 // see if they are each in core.
329 //
330 // Todo: Implement this if you have to!
331 //
332 // We must give up and return false
333 return false;
334}
335
336
337template<class Scalar>
340{
341 if (dim_)
342 return (*vecSpaces_)[0]->smallVecSpcFcty(); // They should all be compatible?
343 return Teuchos::null;
344}
345
346
347template<class Scalar>
350{
351 return defaultProductMultiVector<Scalar>(Teuchos::rcpFromRef(*this),
352 numMembers);
353}
354
355
356template<class Scalar>
359{
360 // Warning! If the client uninitialized this object then changes the
361 // constituent vector spaces then we are in trouble! The client is warned
362 // in documentation!
364 pvs = productVectorSpace<Scalar>();
365 pvs->numBlocks_ = numBlocks_;
366 pvs->vecSpaces_ = vecSpaces_;
367 pvs->vecSpacesOffsets_ = vecSpacesOffsets_;
368 pvs->dim_ = dim_;
369 pvs->isInCore_ = isInCore_;
370 return pvs;
371}
372
373
374// Overridden from Teuchos::Describable
375
376
377template<class Scalar>
379{
380 std::ostringstream oss;
381 oss
383 << "dim="<<dim_
384 << ",numBlocks="<<numBlocks_
385 << "}";
386 return oss.str();
387}
388
389
390template<class Scalar>
392 Teuchos::FancyOStream &out_arg
393 ,const Teuchos::EVerbosityLevel verbLevel
394 ) const
395{
397 using Teuchos::OSTab;
398 RCP<FancyOStream> out = rcpFromRef(out_arg);
399 OSTab tab(out);
400 if (includesVerbLevel(verbLevel, Teuchos::VERB_LOW, true)) {
401 *out << this->description() << std::endl;
402 }
403 if (includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM) && numBlocks_ > 0) {
404 OSTab tab2(out);
405 *out << "Constituent vector spaces V[0], V[1], ... V[numBlocks-1]:\n";
406 OSTab tab3(out);
407 for( int k = 0; k < numBlocks_; ++k ) {
408 *out << "V["<<k<<"] = " << Teuchos::describe(*(*vecSpaces_)[k],verbLevel);
409 }
410 }
411}
412
413
414} // namespace Thyra
415
416
417#endif // THYRA_DEFAULT_PRODUCT_VECTOR_SPACE_HPP
iterator end() const
iterator begin() const
size_type size() const
virtual std::string description() const
Ordinal size() const
Ordinal lbound() const
DefaultProductVectorSpace()
Default construct to uninitialized.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Prints the details about the constituent vector spaces.
RCP< MultiVectorBase< Scalar > > createMembers(int numMembers) const
Returns a DefaultProductMultiVector object.
RCP< const VectorSpaceBase< Scalar > > getBlock(const int k) const
RCP< const VectorSpaceBase< Scalar > > clone() const
Clones the object as promised.
bool hasInCoreView(const Range1D &rng, const EViewType viewType, const EStrideType strideType) const
Returns true if all of the constituent vector spaces return true.
void getVecSpcPoss(Ordinal i, int *kth_vector_space, Ordinal *kth_global_offset) const
Get the position of the vector space object and its offset into a composite vector that owns the ith ...
std::string description() const
Prints just the name DefaultProductVectorSpace along with the overall dimension and the number of blo...
RCP< VectorBase< Scalar > > createMember() const
Returns a DefaultProductVector object.
virtual void initialize(const ArrayView< const RCP< const VectorSpaceBase< Scalar > > > &vecSpaces)
Initialize with a list of constituent vector spaces.
virtual void uninitialize(const ArrayView< RCP< const VectorSpaceBase< Scalar > > > &vecSpaces=Teuchos::null)
Uninitialize.
RCP< const VectorSpaceFactoryBase< Scalar > > smallVecSpcFcty() const
Returns getBlock(0)->smallVecSpcFcty().
bool isCompatible(const VectorSpaceBase< Scalar > &vecSpc) const
Returns true only if also a product vector space and all constituent vectors are compatible.
Ordinal dim() const
Returns the summation of the constituent vector spaces.
void scalarProdsImpl(const MultiVectorBase< Scalar > &X, const MultiVectorBase< Scalar > &Y, const ArrayView< Scalar > &scalarProds) const
Returns the sum of the scalar products of each of the columns of the constituent multi-vectors.
Scalar scalarProd(const VectorBase< Scalar > &x, const VectorBase< Scalar > &y) const
Returns the sum of the scalar products of the constituent vectors.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Interface for a collection of column vectors called a multi-vector.
RCP< const VectorBase< Scalar > > col(Ordinal j) const
Calls colImpl().
Base interface for product multi-vectors.
virtual Teuchos::RCP< const ProductVectorSpaceBase< Scalar > > productSpace() const =0
Returns the associated product vector space that represents the range.
virtual Teuchos::RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const =0
Returns a non-persisting const view of the (zero-based) kth block multi-vector.
Base interface for product vectors.
virtual RCP< const VectorBase< Scalar > > getVectorBlock(const int k) const =0
Returns a non-persisting const view of the (zero-based) kth block vector.
Abstract interface for finite-dimensional dense vectors.
Abstract interface for objects that represent a space for vectors.
virtual Ordinal dim() const =0
Return the dimension of the vector space.
virtual bool isCompatible(const VectorSpaceBase< Scalar > &vecSpc) const =0
Compare the compatibility of two vector spaces.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
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. `*.
@ VIEW_TYPE_DIRECT
The view is a direct view of data and no copies are made.
TypeTo as(const TypeFrom &t)
T_To & dyn_cast(T_From &from)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
TEUCHOSCORE_LIB_DLL_EXPORT bool includesVerbLevel(const EVerbosityLevel verbLevel, const EVerbosityLevel requestedVerbLevel, const bool isDefaultLevel=false)