Thyra  Version of the Day
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>
150 RCP<const VectorSpaceBase<Scalar> >
152 {
153  return range_;
154 }
155 
156 
157 template<class Scalar>
158 RCP<const VectorSpaceBase<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 // Overridden protected functions from LinearOpBase
179 
180 
181 
182 template<class Scalar>
184 {
185  typedef Teuchos::ScalarTraits<Scalar> ST;
186  return ( ST::isComplex ? ( M_trans==NOTRANS || M_trans==CONJTRANS ) : true );
187 }
188 
189 
190 template<class Scalar>
192  const EOpTransp M_trans,
193  const MultiVectorBase<Scalar> &X,
194  const Ptr<MultiVectorBase<Scalar> > &Y,
195  const Scalar alpha,
196  const Scalar beta
197  ) const
198 {
199 #ifdef TEUCHOS_DEBUG
201  "MultiVectorBase<Scalar>::apply()", *this, M_trans, X, &*Y);
202 #endif
203  const Ordinal nc = this->domain()->dim();
204  const Ordinal m = X.domain()->dim();
205  for (Ordinal col_j = 0; col_j < m; ++col_j) {
206  const RCP<const VectorBase<Scalar> > x_j = X.col(col_j);
207  const RCP<VectorBase<Scalar> > y_j = Y->col(col_j);
208  // y_j *= beta
209  Vt_S(y_j.ptr(), beta);
210  // y_j += alpha*op(M)*x_j
211  if(M_trans == NOTRANS) {
212  //
213  // y_j += alpha*M*x_j = alpha*M.col(0)*x_j(0) + ... + alpha*M.col(nc-1)*x_j(nc-1)
214  //
215  // Extract an explicit view of x_j
217  x_j->acquireDetachedView(Range1D(), &x_sub_vec);
218  // Loop through and add the multiple of each column
219  for (Ordinal j = 0; j < nc; ++j )
220  Vp_StV( y_j.ptr(), Scalar(alpha*x_sub_vec(j)), *this->col(j) );
221  // Release the view of x
222  x_j->releaseDetachedView(&x_sub_vec);
223  }
224  else {
225  //
226  // [ alpha*dot(M.col(0),x_j) ]
227  // y_j += alpha*M^T*x_j = [ alpha*dot(M.col(1),x_j) ]
228  // [ ... ]
229  // [ alpha*dot(M.col(nc-1),x_j) ]
230  //
231  // Extract an explicit view of y_j
233  y_j->acquireDetachedView(Range1D(), &y_sub_vec);
234  // Loop through and add to each element in y_j
235  for (Ordinal j = 0; j < nc; ++j )
236  y_sub_vec(j) += alpha*dot(*this->col(j), *x_j);
237  // Commit explicit view of y
238  y_j->commitDetachedView(&y_sub_vec);
239  }
240  }
241 }
242 
243 
244 // Overridden from MultiVectorBase
245 
246 
247 template<class Scalar>
248 RCP<VectorBase<Scalar> >
250 {
251  using Teuchos::as;
252  const int num_cols = col_vecs_.size();
253  TEUCHOS_TEST_FOR_EXCEPTION(
254  !( 0 <= j && j < num_cols ), std::logic_error
255  ,"Error, j = " << j << " does not fall in the range [0,"<<(num_cols-1)<< "]!"
256  );
257  return col_vecs_[j];
258 }
259 
260 
261 template<class Scalar>
262 RCP<MultiVectorBase<Scalar> >
264  const Range1D& col_rng_in
265  )
266 {
267  const Ordinal numCols = domain_->dim();
268  const Range1D col_rng = Teuchos::full_range(col_rng_in,0,numCols-1);
269 #ifdef TEUCHOS_DEBUG
270  TEUCHOS_TEST_FOR_EXCEPTION(
271  !( col_rng.ubound() < numCols ), std::logic_error
272  ,"DefaultColumnwiseMultiVector<Scalar>::subView(col_rng):"
273  "Error, the input range col_rng = ["<<col_rng.lbound()
274  <<","<<col_rng.ubound()<<"] "
275  "is not in the range [0,"<<(numCols-1)<<"]!"
276  );
277 #endif
278  return Teuchos::rcp(
280  range_,
281  domain_->smallVecSpcFcty()->createVecSpc(col_rng.size()),
282  col_vecs_(col_rng.lbound(),col_rng.size())
283  )
284  );
285 }
286 
287 
288 } // end namespace Thyra
289 
290 
291 #endif // THYRA_DEFAULT_COLUMNWISE_MULTI_VECTOR_DEF_HPP
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
#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...
Use the non-transposed operator.
RCP< const VectorSpaceBase< Scalar > > domain() 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.
bool opSupportedImpl(EOpTransp M_trans) const
For complex Scalar types returns true for NOTRANS and CONJTRANS and for real types returns true for a...
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.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Interface for a collection of column vectors called a multi-vector.
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...
RCP< MultiVectorBase< Scalar > > nonconstContigSubViewImpl(const Range1D &col_rng)
RCP< const VectorSpaceBase< Scalar > > range() const
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Teuchos::Range1D Range1D