Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultColumnwiseMultiVector_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_COLUMNWISE_MULTI_VECTOR_DEF_HPP
43#define THYRA_DEFAULT_COLUMNWISE_MULTI_VECTOR_DEF_HPP
44
45#include "Thyra_DefaultColumnwiseMultiVector_decl.hpp"
46#include "Thyra_MultiVectorDefaultBase.hpp"
47#include "Thyra_VectorSpaceBase.hpp"
48#include "Thyra_VectorBase.hpp"
49#include "Thyra_MultiVectorBase.hpp"
50#include "Thyra_VectorStdOps.hpp"
51#include "Thyra_VectorSpaceFactoryBase.hpp"
52#include "Thyra_AssertOp.hpp"
53#include "Teuchos_Assert.hpp"
54#include "Teuchos_as.hpp"
55
56namespace Thyra {
57
58
59// Constructors/Initializers
60
61
62template<class Scalar>
65
66
67template<class Scalar>
69 const RCP<VectorBase<Scalar> > &col_vec
70 )
71{
72 this->initialize(col_vec);
73}
74
75
76template<class Scalar>
78 const RCP<const VectorSpaceBase<Scalar> > &range_in,
79 const RCP<const VectorSpaceBase<Scalar> > &domain_in,
80 const ArrayView<const RCP<VectorBase<Scalar> > > &col_vecs_in
81 )
82{
83 this->initialize(range_in, domain_in, col_vecs_in);
84}
85
86
87template<class Scalar>
89 const RCP<VectorBase<Scalar> > &col_vec
90 )
91{
92#ifdef TEUCHOS_DEBUG
93 const std::string err_msg =
94 "DefaultColumnwiseMultiVector<Scalar>::initialize(...): Error!";
95 TEUCHOS_TEST_FOR_EXCEPT_MSG( is_null(col_vec), err_msg );
96 TEUCHOS_TEST_FOR_EXCEPT_MSG( is_null(col_vec->space()), err_msg );
97#endif
98 range_ = col_vec->space();
99 domain_ = range_->smallVecSpcFcty()->createVecSpc(1);
100 col_vecs_.resize(1);
101 col_vecs_[0] = col_vec;
102}
103
104
105template<class Scalar>
107 const RCP<const VectorSpaceBase<Scalar> > &range_in,
108 const RCP<const VectorSpaceBase<Scalar> > &domain_in,
109 const ArrayView<const RCP<VectorBase<Scalar> > > &col_vecs
110 )
111{
112#ifdef TEUCHOS_DEBUG
113 const std::string err_msg =
114 "DefaultColumnwiseMultiVector<Scalar>::initialize(...): Error!";
115 TEUCHOS_TEST_FOR_EXCEPT_MSG( is_null(range_in), err_msg );
116 TEUCHOS_TEST_FOR_EXCEPT_MSG( is_null(domain_in), err_msg );
117 TEUCHOS_TEST_FOR_EXCEPT_MSG( range_in->dim() == 0, err_msg );
118 TEUCHOS_TEST_FOR_EXCEPT_MSG( domain_in->dim() == 0, err_msg );
119 // ToDo: Check the compatibility of the vectors in col_vecs!
120#endif
121 const int domainDim = domain_in->dim();
122 range_ = range_in;
123 domain_ = domain_in;
124 col_vecs_.clear();
125 col_vecs_.reserve(domainDim);
126 if (col_vecs.size()) {
127 for( Ordinal j = 0; j < domainDim; ++j )
128 col_vecs_.push_back(col_vecs[j]);
129 }
130 else {
131 for( Ordinal j = 0; j < domainDim; ++j )
132 col_vecs_.push_back(createMember(range_));
133 }
134}
135
136
137template<class Scalar>
139{
140 col_vecs_.resize(0);
141 range_ = Teuchos::null;
142 domain_ = Teuchos::null;
143}
144
145
146// Overridden from OpBase
147
148
149template<class Scalar>
152{
153 return range_;
154}
155
156
157template<class Scalar>
160{
161 return domain_;
162}
163
164
165// Overridden protected functions from MultiVectorBase
166
167
168template<class Scalar>
170{
171 const Ordinal m = col_vecs_.size();
172 for (Ordinal col_j = 0; col_j < m; ++col_j) {
173 col_vecs_[col_j]->assign(alpha);
174 }
175}
176
177
178template<class Scalar>
181 )
182{
183#ifdef TEUCHOS_DEBUG
184 TEUCHOS_ASSERT_EQUALITY(col_vecs_.size(), mv.domain()->dim());
185 TEUCHOS_ASSERT(range_->isCompatible(*mv.range()));
186#endif
187 for (Ordinal col_j = 0; col_j < col_vecs_.size(); ++col_j) {
188 col_vecs_[col_j]->assign(*mv.col(col_j));
189 }
190}
191
192
193template<class Scalar>
195{
196 for (Ordinal col_j = 0; col_j < col_vecs_.size(); ++col_j) {
197 col_vecs_[col_j]->scale(alpha);
198 }
199}
200
201
202template<class Scalar>
204 Scalar alpha,
206 )
207{
208#ifdef TEUCHOS_DEBUG
209 TEUCHOS_ASSERT_EQUALITY(col_vecs_.size(), mv.domain()->dim());
210 TEUCHOS_ASSERT(range_->isCompatible(*mv.range()));
211#endif
212 for (Ordinal col_j = 0; col_j < col_vecs_.size(); ++col_j) {
213 col_vecs_[col_j]->update(alpha, *mv.col(col_j));
214 }
215}
216
217
218template<class Scalar>
220 const ArrayView<const Scalar>& alpha,
221 const ArrayView<const Ptr<const MultiVectorBase<Scalar> > >& mv,
222 const Scalar& beta
223 )
224{
225#ifdef TEUCHOS_DEBUG
226 TEUCHOS_ASSERT_EQUALITY(alpha.size(), mv.size());
227 for (Ordinal i = 0; i < mv.size(); ++i) {
228 TEUCHOS_ASSERT_EQUALITY(col_vecs_.size(), mv[i]->domain()->dim());
229 TEUCHOS_ASSERT(range_->isCompatible(*mv[i]->range()));
230 }
231#endif
232 Array<RCP<const VectorBase<Scalar> > > v_rcp(mv.size());
233 Array<Ptr<const VectorBase<Scalar> > > v(mv.size());
234 for (Ordinal col_j = 0; col_j < col_vecs_.size(); ++col_j) {
235 for (Ordinal i = 0; i < mv.size(); ++i) {
236 v_rcp[i] = mv[i]->col(col_j);
237 v[i] = v_rcp[i].ptr();
238 }
239 col_vecs_[col_j]->linear_combination(alpha, v(), beta);
240 }
241}
242
243
244template<class Scalar>
246 const MultiVectorBase<Scalar>& mv,
247 const ArrayView<Scalar>& prods
248 ) const
249{
250#ifdef TEUCHOS_DEBUG
251 TEUCHOS_ASSERT_EQUALITY(col_vecs_.size(), mv.domain()->dim());
252 TEUCHOS_ASSERT_EQUALITY(col_vecs_.size(), prods.size());
253 TEUCHOS_ASSERT(range_->isCompatible(*mv.range()));
254#endif
255 for (Ordinal col_j = 0; col_j < col_vecs_.size(); ++col_j) {
256 prods[col_j] = col_vecs_[col_j]->dot(*mv.col(col_j));
257 }
258}
259
260
261template<class Scalar>
264 ) const
265{
266#ifdef TEUCHOS_DEBUG
267 TEUCHOS_ASSERT_EQUALITY(col_vecs_.size(), norms.size());
268#endif
269 for (Ordinal col_j = 0; col_j < col_vecs_.size(); ++col_j) {
270 norms[col_j] = col_vecs_[col_j]->norm_1();
271 }
272}
273
274
275template<class Scalar>
278 ) const
279{
280#ifdef TEUCHOS_DEBUG
281 TEUCHOS_ASSERT_EQUALITY(col_vecs_.size(), norms.size());
282#endif
283 for (Ordinal col_j = 0; col_j < col_vecs_.size(); ++col_j) {
284 norms[col_j] = col_vecs_[col_j]->norm_2();
285 }
286}
287
288
289template<class Scalar>
292 ) const
293{
294#ifdef TEUCHOS_DEBUG
295 TEUCHOS_ASSERT_EQUALITY(col_vecs_.size(), norms.size());
296#endif
297 for (Ordinal col_j = 0; col_j < col_vecs_.size(); ++col_j) {
298 norms[col_j] = col_vecs_[col_j]->norm_inf();
299 }
300}
301
302
303// Overridden protected functions from LinearOpBase
304
305
306
307template<class Scalar>
309{
311 return ( ST::isComplex ? ( M_trans==NOTRANS || M_trans==CONJTRANS ) : true );
312}
313
314
315template<class Scalar>
317 const EOpTransp M_trans,
319 const Ptr<MultiVectorBase<Scalar> > &Y,
320 const Scalar alpha,
321 const Scalar beta
322 ) const
323{
324#ifdef TEUCHOS_DEBUG
326 "MultiVectorBase<Scalar>::apply()", *this, M_trans, X, &*Y);
327#endif
328 const Ordinal nc = this->domain()->dim();
329 const Ordinal m = X.domain()->dim();
330 for (Ordinal col_j = 0; col_j < m; ++col_j) {
331 const RCP<const VectorBase<Scalar> > x_j = X.col(col_j);
332 const RCP<VectorBase<Scalar> > y_j = Y->col(col_j);
333 // y_j *= beta
334 Vt_S(y_j.ptr(), beta);
335 // y_j += alpha*op(M)*x_j
336 if(M_trans == NOTRANS) {
337 //
338 // y_j += alpha*M*x_j = alpha*M.col(0)*x_j(0) + ... + alpha*M.col(nc-1)*x_j(nc-1)
339 //
340 // Extract an explicit view of x_j
342 x_j->acquireDetachedView(Range1D(), &x_sub_vec);
343 // Loop through and add the multiple of each column
344 for (Ordinal j = 0; j < nc; ++j )
345 Vp_StV( y_j.ptr(), Scalar(alpha*x_sub_vec(j)), *this->col(j) );
346 // Release the view of x
347 x_j->releaseDetachedView(&x_sub_vec);
348 }
349 else {
350 //
351 // [ alpha*dot(M.col(0),x_j) ]
352 // y_j += alpha*M^T*x_j = [ alpha*dot(M.col(1),x_j) ]
353 // [ ... ]
354 // [ alpha*dot(M.col(nc-1),x_j) ]
355 //
356 // Extract an explicit view of y_j
358 y_j->acquireDetachedView(Range1D(), &y_sub_vec);
359 // Loop through and add to each element in y_j
360 for (Ordinal j = 0; j < nc; ++j )
361 y_sub_vec(j) += alpha*dot(*this->col(j), *x_j);
362 // Commit explicit view of y
363 y_j->commitDetachedView(&y_sub_vec);
364 }
365 }
366}
367
368
369// Overridden from MultiVectorBase
370
371
372template<class Scalar>
375{
376 using Teuchos::as;
377 const int num_cols = col_vecs_.size();
379 !( 0 <= j && j < num_cols ), std::logic_error
380 ,"Error, j = " << j << " does not fall in the range [0,"<<(num_cols-1)<< "]!"
381 );
382 return col_vecs_[j];
383}
384
385
386template<class Scalar>
389 const Range1D& col_rng_in
390 )
391{
392 const Ordinal numCols = domain_->dim();
393 const Range1D col_rng = Teuchos::full_range(col_rng_in,0,numCols-1);
394#ifdef TEUCHOS_DEBUG
396 !( col_rng.ubound() < numCols ), std::logic_error
397 ,"DefaultColumnwiseMultiVector<Scalar>::subView(col_rng):"
398 "Error, the input range col_rng = ["<<col_rng.lbound()
399 <<","<<col_rng.ubound()<<"] "
400 "is not in the range [0,"<<(numCols-1)<<"]!"
401 );
402#endif
403 return Teuchos::rcp(
405 range_,
406 domain_->smallVecSpcFcty()->createVecSpc(col_rng.size()),
407 col_vecs_(col_rng.lbound(),col_rng.size())
408 )
409 );
410}
411
412
413} // end namespace Thyra
414
415
416#endif // THYRA_DEFAULT_COLUMNWISE_MULTI_VECTOR_DEF_HPP
size_type size() const
Ptr< T > ptr() const
Ordinal size() const
Ordinal lbound() const
Ordinal ubound() const
Default subclass for MultiVectorBase implemented using columns of separate abstract vectors.
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
This function is implemented in terms of the multi-vector applyOp() function.
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
RCP< MultiVectorBase< Scalar > > nonconstContigSubViewImpl(const Range1D &col_rng)
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
virtual void norms2Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
virtual void normsInfImpl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
virtual void norms1Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
RCP< const VectorSpaceBase< Scalar > > domain() const
void initialize(const RCP< VectorBase< Scalar > > &col_vec)
Initialize given a single vector object.
RCP< const VectorSpaceBase< Scalar > > range() const
bool opSupportedImpl(EOpTransp M_trans) const
For complex Scalar types returns true for NOTRANS and CONJTRANS and for real types returns true for a...
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
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().
Abstract interface for finite-dimensional dense vectors.
Abstract interface for objects that represent a space for vectors.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
#define THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(FUNC_NAME, M, M_T, X, Y)
This is a very useful macro that should be used to validate that the spaces for the multi-vector vers...
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Teuchos::Range1D Range1D
@ NOTRANS
Use the non-transposed operator.
@ CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types).
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)