Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
56 namespace Thyra {
57 
58 
59 // Constructors/Initializers
60 
61 
62 template<class Scalar>
64 {}
65 
66 
67 template<class Scalar>
69  const RCP<VectorBase<Scalar> > &col_vec
70  )
71 {
72  this->initialize(col_vec);
73 }
74 
75 
76 template<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 
87 template<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 
105 template<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 
137 template<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 
149 template<class Scalar>
152 {
153  return range_;
154 }
155 
156 
157 template<class Scalar>
160 {
161  return domain_;
162 }
163 
164 
165 // Overridden protected functions from MultiVectorBase
166 
167 
168 template<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 
178 template<class Scalar>
180  const MultiVectorBase<Scalar>& mv
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 
193 template<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 
202 template<class Scalar>
204  Scalar alpha,
205  const MultiVectorBase<Scalar>& mv
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 
218 template<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 
244 template<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 
261 template<class Scalar>
263  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
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 
275 template<class Scalar>
277  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
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 
289 template<class Scalar>
291  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
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 
307 template<class Scalar>
309 {
311  return ( ST::isComplex ? ( M_trans==NOTRANS || M_trans==CONJTRANS ) : true );
312 }
313 
314 
315 template<class Scalar>
317  const EOpTransp M_trans,
318  const MultiVectorBase<Scalar> &X,
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 
372 template<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 
386 template<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
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
RCP< const VectorSpaceBase< Scalar > > range() const
bool is_null(const boost::shared_ptr< T > &p)
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
virtual void norms1Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
#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...
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 norms2Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Use the non-transposed operator.
size_type size() const
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
Abstract interface for objects that represent a space for vectors.
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
Ordinal ubound() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
#define TEUCHOS_TEST_FOR_EXCEPT_MSG(throw_exception_test, msg)
Interface for a collection of column vectors called a multi-vector.
Ptr< T > ptr() const
RCP< const VectorBase< Scalar > > col(Ordinal j) const
Calls colImpl().
RCP< VectorBase< Scalar > > nonconstColImpl(Ordinal j)
Abstract interface for finite-dimensional dense vectors.
void initialize(const RCP< VectorBase< Scalar > > &col_vec)
Initialize given a single vector object.
Default subclass for MultiVectorBase implemented using columns of separate abstract vectors...
virtual void normsInfImpl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
RCP< MultiVectorBase< Scalar > > nonconstContigSubViewImpl(const Range1D &col_rng)
bool opSupportedImpl(EOpTransp M_trans) const
For complex Scalar types returns true for NOTRANS and CONJTRANS and for real types returns true for a...
Ordinal lbound() const
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
RCP< const VectorSpaceBase< Scalar > > domain() const
Ordinal size() const
TypeTo as(const TypeFrom &t)
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
Teuchos::Range1D Range1D