Belos Package Browser (Single Doxygen Collection)  Development
MyMultiVec.hpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ************************************************************************
4 //
5 // Belos: Block Linear Solvers Package
6 // Copyright 2004 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 */
43 
44 #ifndef MY_MULTIVECTOR_HPP
45 #define MY_MULTIVECTOR_HPP
46 
47 #include "BelosConfigDefs.hpp"
48 #include "BelosMultiVec.hpp"
49 
51 
61 template <class ScalarType>
62 class MyMultiVec : public Belos::MultiVec<ScalarType>
63 {
64 public:
65 
67  MyMultiVec (const ptrdiff_t Length, const int NumberVecs) :
68  Length_ (Length),
69  NumberVecs_ (NumberVecs)
70  {
71  Check();
72 
73  data_.resize(NumberVecs);
74  ownership_.resize(NumberVecs);
75 
76  // Allocates memory to store the vectors.
77  for (int v = 0; v < NumberVecs_; ++v) {
78  data_[v] = new ScalarType[Length];
79  ownership_[v] = true;
80  }
81 
82  // Initializes all elements to zero.
83  MvInit(0.0);
84  }
85 
87  MyMultiVec(const ptrdiff_t Length, const std::vector<ScalarType*>& rhs) :
88  Length_(Length),
89  NumberVecs_(rhs.size())
90  {
91  Check();
92 
93  data_.resize(NumberVecs_);
94  ownership_.resize(NumberVecs_);
95 
96  // Copies pointers from input array, set ownership so that
97  // this memory is not free'd in the destructor
98  for (int v = 0; v < NumberVecs_; ++v) {
99  data_[v] = rhs[v];
100  ownership_[v] = false;
101  }
102  }
103 
105  MyMultiVec(const MyMultiVec& rhs) :
106  Length_(rhs.GetGlobalLength()),
108  {
109  Check();
110 
111  data_.resize(NumberVecs_);
112  ownership_.resize(NumberVecs_);
113 
114  for (int v = 0 ; v < NumberVecs_ ; ++v)
115  {
116  data_[v] = new ScalarType[Length_];
117  ownership_[v] = true;
118  }
119 
120  for (int v = 0 ; v < NumberVecs_ ; ++v)
121  {
122  for (int i = 0 ; i < Length_ ; ++i)
123  (*this)(i, v) = rhs(i, v);
124  }
125  }
126 
129  {
130  for (int v = 0 ; v < NumberVecs_ ; ++v)
131  if (ownership_[v])
132  delete[] data_[v];
133  }
134 
136  MyMultiVec* Clone(const int NumberVecs) const
137  {
138  // FIXME
139  MyMultiVec* tmp = new MyMultiVec(Length_, NumberVecs);
140 
141  // for (int v = 0 ; v < NumberVecs ; ++v)
142  // for (int i = 0 ; i < Length_ ; ++i)
143  // (*tmp)(i, v) = (*this)(i, v);
144 
145  return(tmp);
146  }
147 
148  // Returns a clone of the corrent multi-std::vector.
150  {
151  return(new MyMultiVec(*this));
152  }
153 
155  MyMultiVec* CloneCopy(const std::vector< int > &index) const
156  {
157  int size = index.size();
158  MyMultiVec* tmp = new MyMultiVec(Length_, size);
159 
160  for (size_t v = 0 ; v < index.size() ; ++v) {
161  for (int i = 0 ; i < Length_ ; ++i) {
162  (*tmp)(i, v) = (*this)(i, index[v]);
163  }
164  }
165 
166  return(tmp);
167  }
168 
170  MyMultiVec* CloneViewNonConst(const std::vector< int > &index)
171  {
172  int size = index.size();
173  std::vector<ScalarType*> values(size);
174 
175  for (size_t v = 0 ; v < index.size() ; ++v)
176  values[v] = data_[index[v]];
177 
178  return(new MyMultiVec(Length_, values));
179  }
180 
182  const MyMultiVec* CloneView (const std::vector<int>& index) const
183  {
184  int size = index.size();
185  std::vector<ScalarType*> values (size);
186 
187  for (size_t v = 0; v < index.size (); ++v) {
188  values[v] = data_[index[v]];
189  }
190 
191  return(new MyMultiVec(Length_, values));
192  }
193 
194  ptrdiff_t GetGlobalLength () const
195  {
196  return Length_;
197  }
198 
199  int GetNumberVecs () const
200  {
201  return(NumberVecs_);
202  }
203 
204  // Update *this with alpha * A * B + beta * (*this).
205  void MvTimesMatAddMv (const ScalarType alpha, const Belos::MultiVec<ScalarType> &A,
207  const ScalarType beta)
208  {
209  assert (Length_ == A.GetGlobalLength());
210  assert (B.numRows() == A.GetNumberVecs());
211  assert (B.numCols() <= NumberVecs_);
212 
213  MyMultiVec* MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
214  TEUCHOS_ASSERT(MyA != NULL);
215 
216  if ((*this)[0] == (*MyA)[0]) {
217  // If this == A, then need additional storage ...
218  // This situation is a bit peculiar but it may be required by
219  // certain algorithms.
220 
221  std::vector<ScalarType> tmp(NumberVecs_);
222 
223  for (int i = 0 ; i < Length_ ; ++i) {
224  for (int v = 0; v < A.GetNumberVecs() ; ++v) {
225  tmp[v] = (*MyA)(i, v);
226  }
227 
228  for (int v = 0 ; v < B.numCols() ; ++v) {
229  (*this)(i, v) *= beta;
230  ScalarType res = Teuchos::ScalarTraits<ScalarType>::zero();
231 
232  for (int j = 0 ; j < A.GetNumberVecs() ; ++j) {
233  res += tmp[j] * B(j, v);
234  }
235 
236  (*this)(i, v) += alpha * res;
237  }
238  }
239  }
240  else {
241  for (int i = 0 ; i < Length_ ; ++i) {
242  for (int v = 0 ; v < B.numCols() ; ++v) {
243  (*this)(i, v) *= beta;
244  ScalarType res = 0.0;
245  for (int j = 0 ; j < A.GetNumberVecs() ; ++j) {
246  res += (*MyA)(i, j) * B(j, v);
247  }
248 
249  (*this)(i, v) += alpha * res;
250  }
251  }
252  }
253  }
254 
255  // Replace *this with alpha * A + beta * B.
256  void MvAddMv (const ScalarType alpha, const Belos::MultiVec<ScalarType>& A,
257  const ScalarType beta, const Belos::MultiVec<ScalarType>& B)
258  {
259  MyMultiVec* MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
260  TEUCHOS_ASSERT(MyA != NULL);
261 
262  MyMultiVec* MyB = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(B));
263  TEUCHOS_ASSERT(MyB != NULL);
264 
265  assert (NumberVecs_ == A.GetNumberVecs());
266  assert (NumberVecs_ == B.GetNumberVecs());
267 
268  assert (Length_ == A.GetGlobalLength());
269  assert (Length_ == B.GetGlobalLength());
270 
271  for (int v = 0 ; v < NumberVecs_ ; ++v) {
272  for (int i = 0 ; i < Length_ ; ++i) {
273  (*this)(i, v) = alpha * (*MyA)(i, v) + beta * (*MyB)(i, v);
274  }
275  }
276  }
277 
278  // Replace each element of the vectors in *this with (*this)*alpha.
279  void MvScale (const ScalarType alpha)
280  {
281  for (int v = 0 ; v < NumberVecs_ ; ++v) {
282  for (int i = 0 ; i < Length_ ; ++i) {
283  (*this)(i, v) *= alpha;
284  }
285  }
286  }
287 
288  // Replace each element of the vectors in *this with (*this)*alpha[i].
289  void MvScale (const std::vector<ScalarType>& alpha)
290  {
291  assert((int)alpha.size() >= NumberVecs_);
292  for (int v = 0 ; v < NumberVecs_ ; ++v) {
293  for (int i = 0 ; i < Length_ ; ++i) {
294  (*this)(i, v) *= alpha[v];
295  }
296  }
297  }
298 
299  // Compute a dense matrix B through the matrix-matrix multiply alpha * A^H * (*this).
300  void MvTransMv (const ScalarType alpha, const Belos::MultiVec<ScalarType>& A,
302  {
303  MyMultiVec* MyA;
304  MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
305  TEUCHOS_ASSERT(MyA != NULL);
306 
307  assert (A.GetGlobalLength() == Length_);
308  assert (NumberVecs_ <= B.numCols());
309  assert (A.GetNumberVecs() <= B.numRows());
310 
311  for (int v = 0 ; v < A.GetNumberVecs() ; ++v) {
312  for (int w = 0 ; w < NumberVecs_ ; ++w) {
313  ScalarType value = 0.0;
314  for (int i = 0 ; i < Length_ ; ++i) {
315  value += Teuchos::ScalarTraits<ScalarType>::conjugate((*MyA)(i, v)) * (*this)(i, w);
316  }
317  B(v, w) = alpha * value;
318  }
319  }
320  }
321 
322 
323  // Compute a std::vector b where the components are the individual dot-products, i.e.b[i] = A[i]^H*this[i] where A[i] is the i-th column of A.
324  void MvDot (const Belos::MultiVec<ScalarType>& A, std::vector<ScalarType>& b) const
325  {
326  MyMultiVec* MyA;
327  MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
328  TEUCHOS_ASSERT(MyA != NULL);
329 
330  assert (NumberVecs_ <= (int)b.size());
331  assert (NumberVecs_ == A.GetNumberVecs());
332  assert (Length_ == A.GetGlobalLength());
333 
334  for (int v = 0 ; v < NumberVecs_ ; ++v) {
335  ScalarType value = 0.0;
336  for (int i = 0 ; i < Length_ ; ++i) {
337  value += (*this)(i, v) * Teuchos::ScalarTraits<ScalarType>::conjugate((*MyA)(i, v));
338  }
339  b[v] = value;
340  }
341  }
342 
343 
344  void MvNorm ( std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec,
345  Belos::NormType type = Belos::TwoNorm ) const
346  {
347  assert (NumberVecs_ <= (int)normvec.size());
348 
349  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
350 
351  for (int v = 0 ; v < NumberVecs_ ; ++v) {
352  MagnitudeType value = Teuchos::ScalarTraits<MagnitudeType>::zero();
353  for (int i = 0 ; i < Length_ ; ++i) {
354  MagnitudeType val = Teuchos::ScalarTraits<ScalarType>::magnitude((*this)(i, v));
355  value += val * val;
356  }
358  }
359  }
360 
361  // Copy the vectors in A to a set of vectors in *this. The numvecs vectors in
362  // A are copied to a subset of vectors in *this indicated by the indices given
363  // in index.
364  // FIXME: not so clear what the size of A and index.size() are...
366  const std::vector<int> &index)
367  {
368  MyMultiVec* MyA;
369  MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
370  TEUCHOS_ASSERT(MyA != NULL);
371 
372  assert (A.GetNumberVecs() >= (int)index.size());
373  assert (A.GetGlobalLength() == Length_);
374 
375  for (unsigned int v = 0 ; v < index.size() ; ++v) {
376  for (int i = 0 ; i < Length_ ; ++i) {
377  (*this)(i, index[v]) = (*MyA)(i, v);
378  }
379  }
380  }
381 
382  // Fill the vectors in *this with random numbers.
383  void MvRandom ()
384  {
385  for (int v = 0 ; v < NumberVecs_ ; ++v) {
386  for (int i = 0 ; i < Length_ ; ++i) {
388  }
389  }
390  }
391 
392  // Replace each element of the vectors in *this with alpha.
393  void MvInit (const ScalarType alpha)
394  {
395  for (int v = 0 ; v < NumberVecs_ ; ++v) {
396  for (int i = 0 ; i < Length_ ; ++i) {
397  (*this)(i, v) = alpha;
398  }
399  }
400  }
401 
402  void MvPrint (std::ostream &os) const
403  {
404  os << "Object MyMultiVec" << std::endl;
405  os << "Number of rows = " << Length_ << std::endl;
406  os << "Number of vecs = " << NumberVecs_ << std::endl;
407 
408  for (int i = 0 ; i < Length_ ; ++i)
409  {
410  for (int v = 0 ; v < NumberVecs_ ; ++v)
411  os << (*this)(i, v) << " ";
412  os << std::endl;
413  }
414  }
415 
416  inline ScalarType& operator()(const int i, const int j)
417  {
418  if (j < 0 || j >= NumberVecs_) throw(-1);
419  if (i < 0 || i >= Length_) throw(-2);
420  //
421  return(data_[j][i]);
422  }
423 
424  inline const ScalarType& operator()(const int i, const int j) const
425  {
426  if (j < 0 || j >= NumberVecs_) throw(-1);
427  if (i < 0 || i >= Length_) throw(-2);
428  //
429  return(data_[j][i]);
430  }
431 
432  ScalarType* operator[](int v)
433  {
434  if (v < 0 || v >= NumberVecs_) throw(-1);
435  return(data_[v]);
436  }
437 
438  ScalarType* operator[](int v) const
439  {
440  return(data_[v]);
441  }
442 
443 private:
444  void Check()
445  {
446  if (Length_ <= 0)
447  throw("Length must be positive");
448 
449  if (NumberVecs_ <= 0)
450  throw("Number of vectors must be positive");
451  }
452 
454  const ptrdiff_t Length_;
456  const int NumberVecs_;
458  std::vector<ScalarType*> data_;
460  std::vector<bool> ownership_;
461 
462 };
463 
464 
465 #endif // MY_MULTIVECTOR_HPP
void Check()
Definition: MyMultiVec.hpp:444
MyMultiVec(const MyMultiVec &rhs)
Copy constructor, performs a deep copy.
Definition: MyMultiVec.hpp:105
std::vector< ScalarType * > data_
Pointers to the storage of the vectors.
Definition: MyMultiVec.hpp:458
MyMultiVec * CloneCopy() const
Create a new MultiVec and copy contents of *this into it (deep copy).
Definition: MyMultiVec.hpp:149
const ScalarType & operator()(const int i, const int j) const
Definition: MyMultiVec.hpp:424
ScalarType & operator()(const int i, const int j)
Definition: MyMultiVec.hpp:416
void MvNorm(std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, Belos::NormType type=Belos::TwoNorm) const
Compute the norm of each vector in *this.
Definition: MyMultiVec.hpp:344
static T squareroot(T x)
void MvPrint(std::ostream &os) const
Print *this multivector to the os output stream.
Definition: MyMultiVec.hpp:402
void MvScale(const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
Definition: MyMultiVec.hpp:279
~MyMultiVec()
Destructor.
Definition: MyMultiVec.hpp:128
MyMultiVec * CloneViewNonConst(const std::vector< int > &index)
Returns a view of current std::vector (shallow copy)
Definition: MyMultiVec.hpp:170
const MyMultiVec * CloneView(const std::vector< int > &index) const
Returns a view of current std::vector (shallow copy), const version.
Definition: MyMultiVec.hpp:182
void MvTimesMatAddMv(const ScalarType alpha, const Belos::MultiVec< ScalarType > &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta)
Update *this with alpha * A * B + beta * (*this).
Definition: MyMultiVec.hpp:205
Simple example of a user&#39;s defined Belos::MultiVec class.
Definition: MyMultiVec.hpp:62
static T conjugate(T a)
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
Definition: MyMultiVec.hpp:199
MyMultiVec(const ptrdiff_t Length, const std::vector< ScalarType *> &rhs)
Constructor with already allocated memory.
Definition: MyMultiVec.hpp:87
MyMultiVec * CloneCopy(const std::vector< int > &index) const
Returns a clone copy of specified vectors.
Definition: MyMultiVec.hpp:155
ScalarType * operator[](int v)
Definition: MyMultiVec.hpp:432
static magnitudeType magnitude(T a)
NormType
The type of vector norm to compute.
Definition: BelosTypes.hpp:97
void MvRandom()
Fill all the vectors in *this with random numbers.
Definition: MyMultiVec.hpp:383
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
Definition: MyMultiVec.hpp:194
const int NumberVecs_
Number of multi-vectors.
Definition: MyMultiVec.hpp:456
void MvTransMv(const ScalarType alpha, const Belos::MultiVec< ScalarType > &A, Teuchos::SerialDenseMatrix< int, ScalarType > &B) const
Compute a dense matrix B through the matrix-matrix multiply alpha * A^T * (*this).
Definition: MyMultiVec.hpp:300
Interface for multivectors used by Belos&#39; linear solvers.
#define TEUCHOS_ASSERT(assertion_test)
void MvDot(const Belos::MultiVec< ScalarType > &A, std::vector< ScalarType > &b) const
Compute the dot product of each column of *this with the corresponding column of A.
Definition: MyMultiVec.hpp:324
void MvScale(const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
Definition: MyMultiVec.hpp:289
const ptrdiff_t Length_
Length of the vectors.
Definition: MyMultiVec.hpp:454
Belos header file which uses auto-configuration information to include necessary C++ headers...
MyMultiVec * Clone(const int NumberVecs) const
Returns a clone of the current std::vector.
Definition: MyMultiVec.hpp:136
void SetBlock(const Belos::MultiVec< ScalarType > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
Definition: MyMultiVec.hpp:365
std::vector< bool > ownership_
If true, then this object owns the vectors and must free them in dtor.
Definition: MyMultiVec.hpp:460
MyMultiVec(const ptrdiff_t Length, const int NumberVecs)
Constructor for a NumberVecs vectors of length Length.
Definition: MyMultiVec.hpp:67
void MvAddMv(const ScalarType alpha, const Belos::MultiVec< ScalarType > &A, const ScalarType beta, const Belos::MultiVec< ScalarType > &B)
Replace *this with alpha * A + beta * B.
Definition: MyMultiVec.hpp:256
ScalarType * operator[](int v) const
Definition: MyMultiVec.hpp:438
void MvInit(const ScalarType alpha)
Replace each element of the vectors in *this with alpha.
Definition: MyMultiVec.hpp:393
Interface for multivectors used by Belos&#39; linear solvers.