Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
50#include "Teuchos_SerialDenseMatrix.hpp"
51#include "Teuchos_SerialDenseHelpers.hpp"
52
54
64template <class ScalarType>
65class MyMultiVec : public Belos::MultiVec<ScalarType>
66{
67public:
68
70 MyMultiVec (const ptrdiff_t Length, const int NumberVecs) :
71 Length_ (Length),
72 NumberVecs_ (NumberVecs)
73 {
74 Check();
75
76 data_.resize(NumberVecs);
77 ownership_.resize(NumberVecs);
78
79 // Allocates memory to store the vectors.
80 for (int v = 0; v < NumberVecs_; ++v) {
81 data_[v] = new ScalarType[Length];
82 ownership_[v] = true;
83 }
84
85 // Initializes all elements to zero.
86 MvInit(0.0);
87 }
88
90 MyMultiVec(const ptrdiff_t Length, const std::vector<ScalarType*>& rhs) :
91 Length_(Length),
92 NumberVecs_(rhs.size())
93 {
94 Check();
95
96 data_.resize(NumberVecs_);
97 ownership_.resize(NumberVecs_);
98
99 // Copies pointers from input array, set ownership so that
100 // this memory is not free'd in the destructor
101 for (int v = 0; v < NumberVecs_; ++v) {
102 data_[v] = rhs[v];
103 ownership_[v] = false;
104 }
105 }
106
108 MyMultiVec(const MyMultiVec& rhs) :
111 {
112 Check();
113
114 data_.resize(NumberVecs_);
115 ownership_.resize(NumberVecs_);
116
117 for (int v = 0 ; v < NumberVecs_ ; ++v)
118 {
119 data_[v] = new ScalarType[Length_];
120 ownership_[v] = true;
121 }
122
123 for (int v = 0 ; v < NumberVecs_ ; ++v)
124 {
125 for (int i = 0 ; i < Length_ ; ++i)
126 (*this)(i, v) = rhs(i, v);
127 }
128 }
129
132 {
133 for (int v = 0 ; v < NumberVecs_ ; ++v)
134 if (ownership_[v])
135 delete[] data_[v];
136 }
137
139 MyMultiVec* Clone(const int NumberVecs) const
140 {
141 // FIXME
142 MyMultiVec* tmp = new MyMultiVec(Length_, NumberVecs);
143
144 // for (int v = 0 ; v < NumberVecs ; ++v)
145 // for (int i = 0 ; i < Length_ ; ++i)
146 // (*tmp)(i, v) = (*this)(i, v);
147
148 return(tmp);
149 }
150
151 // Returns a clone of the corrent multi-std::vector.
153 {
154 return(new MyMultiVec(*this));
155 }
156
158 MyMultiVec* CloneCopy(const std::vector< int > &index) const
159 {
160 int size = index.size();
161 MyMultiVec* tmp = new MyMultiVec(Length_, size);
162
163 for (size_t v = 0 ; v < index.size() ; ++v) {
164 for (int i = 0 ; i < Length_ ; ++i) {
165 (*tmp)(i, v) = (*this)(i, index[v]);
166 }
167 }
168
169 return(tmp);
170 }
171
173 MyMultiVec* CloneViewNonConst(const std::vector< int > &index)
174 {
175 int size = index.size();
176 std::vector<ScalarType*> values(size);
177
178 for (size_t v = 0 ; v < index.size() ; ++v)
179 values[v] = data_[index[v]];
180
181 return(new MyMultiVec(Length_, values));
182 }
183
185 const MyMultiVec* CloneView (const std::vector<int>& index) const
186 {
187 int size = index.size();
188 std::vector<ScalarType*> values (size);
189
190 for (size_t v = 0; v < index.size (); ++v) {
191 values[v] = data_[index[v]];
192 }
193
194 return(new MyMultiVec(Length_, values));
195 }
196
197 ptrdiff_t GetGlobalLength () const
198 {
199 return Length_;
200 }
201
202 int GetNumberVecs () const
203 {
204 return(NumberVecs_);
205 }
206
207 // Update *this with alpha * A * B + beta * (*this).
208 void MvTimesMatAddMv (const ScalarType alpha, const Belos::MultiVec<ScalarType> &A,
209 const Teuchos::SerialDenseMatrix<int, ScalarType> &B,
210 const ScalarType beta)
211 {
212 assert (Length_ == A.GetGlobalLength());
213 assert (B.numRows() == A.GetNumberVecs());
214 assert (B.numCols() <= NumberVecs_);
215
216 MyMultiVec* MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
217 TEUCHOS_ASSERT(MyA != NULL);
218
219 if ((*this)[0] == (*MyA)[0]) {
220 // If this == A, then need additional storage ...
221 // This situation is a bit peculiar but it may be required by
222 // certain algorithms.
223
224 std::vector<ScalarType> tmp(NumberVecs_);
225
226 for (int i = 0 ; i < Length_ ; ++i) {
227 for (int v = 0; v < A.GetNumberVecs() ; ++v) {
228 tmp[v] = (*MyA)(i, v);
229 }
230
231 for (int v = 0 ; v < B.numCols() ; ++v) {
232 (*this)(i, v) *= beta;
233 ScalarType res = Teuchos::ScalarTraits<ScalarType>::zero();
234
235 for (int j = 0 ; j < A.GetNumberVecs() ; ++j) {
236 res += tmp[j] * B(j, v);
237 }
238
239 (*this)(i, v) += alpha * res;
240 }
241 }
242 }
243 else {
244 for (int i = 0 ; i < Length_ ; ++i) {
245 for (int v = 0 ; v < B.numCols() ; ++v) {
246 (*this)(i, v) *= beta;
247 ScalarType res = 0.0;
248 for (int j = 0 ; j < A.GetNumberVecs() ; ++j) {
249 res += (*MyA)(i, j) * B(j, v);
250 }
251
252 (*this)(i, v) += alpha * res;
253 }
254 }
255 }
256 }
257
258 // Replace *this with alpha * A + beta * B.
259 void MvAddMv (const ScalarType alpha, const Belos::MultiVec<ScalarType>& A,
260 const ScalarType beta, const Belos::MultiVec<ScalarType>& B)
261 {
262 MyMultiVec* MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
263 TEUCHOS_ASSERT(MyA != NULL);
264
265 MyMultiVec* MyB = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(B));
266 TEUCHOS_ASSERT(MyB != NULL);
267
268 assert (NumberVecs_ == A.GetNumberVecs());
269 assert (NumberVecs_ == B.GetNumberVecs());
270
271 assert (Length_ == A.GetGlobalLength());
272 assert (Length_ == B.GetGlobalLength());
273
274 for (int v = 0 ; v < NumberVecs_ ; ++v) {
275 for (int i = 0 ; i < Length_ ; ++i) {
276 (*this)(i, v) = alpha * (*MyA)(i, v) + beta * (*MyB)(i, v);
277 }
278 }
279 }
280
281 // Replace each element of the vectors in *this with (*this)*alpha.
282 void MvScale (const ScalarType alpha)
283 {
284 for (int v = 0 ; v < NumberVecs_ ; ++v) {
285 for (int i = 0 ; i < Length_ ; ++i) {
286 (*this)(i, v) *= alpha;
287 }
288 }
289 }
290
291 // Replace each element of the vectors in *this with (*this)*alpha[i].
292 void MvScale (const std::vector<ScalarType>& alpha)
293 {
294 assert((int)alpha.size() >= NumberVecs_);
295 for (int v = 0 ; v < NumberVecs_ ; ++v) {
296 for (int i = 0 ; i < Length_ ; ++i) {
297 (*this)(i, v) *= alpha[v];
298 }
299 }
300 }
301
302 // Compute a dense matrix B through the matrix-matrix multiply alpha * A^H * (*this).
303 void MvTransMv (const ScalarType alpha, const Belos::MultiVec<ScalarType>& A,
304 Teuchos::SerialDenseMatrix< int, ScalarType >& B) const
305 {
306 MyMultiVec* MyA;
307 MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
308 TEUCHOS_ASSERT(MyA != NULL);
309
310 assert (A.GetGlobalLength() == Length_);
311 assert (NumberVecs_ <= B.numCols());
312 assert (A.GetNumberVecs() <= B.numRows());
313
314 for (int v = 0 ; v < A.GetNumberVecs() ; ++v) {
315 for (int w = 0 ; w < NumberVecs_ ; ++w) {
316 ScalarType value = 0.0;
317 for (int i = 0 ; i < Length_ ; ++i) {
318 value += Teuchos::ScalarTraits<ScalarType>::conjugate((*MyA)(i, v)) * (*this)(i, w);
319 }
320 B(v, w) = alpha * value;
321 }
322 }
323 }
324
325
326 // 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.
327 void MvDot (const Belos::MultiVec<ScalarType>& A, std::vector<ScalarType>& b) const
328 {
329 MyMultiVec* MyA;
330 MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
331 TEUCHOS_ASSERT(MyA != NULL);
332
333 assert (NumberVecs_ <= (int)b.size());
334 assert (NumberVecs_ == A.GetNumberVecs());
335 assert (Length_ == A.GetGlobalLength());
336
337 for (int v = 0 ; v < NumberVecs_ ; ++v) {
338 ScalarType value = 0.0;
339 for (int i = 0 ; i < Length_ ; ++i) {
340 value += (*this)(i, v) * Teuchos::ScalarTraits<ScalarType>::conjugate((*MyA)(i, v));
341 }
342 b[v] = value;
343 }
344 }
345
346
347 void MvNorm ( std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec,
348 Belos::NormType type = Belos::TwoNorm ) const
349 {
350 assert (NumberVecs_ <= (int)normvec.size());
351
352 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
353
354 for (int v = 0 ; v < NumberVecs_ ; ++v) {
355 MagnitudeType value = Teuchos::ScalarTraits<MagnitudeType>::zero();
356 for (int i = 0 ; i < Length_ ; ++i) {
357 MagnitudeType val = Teuchos::ScalarTraits<ScalarType>::magnitude((*this)(i, v));
358 value += val * val;
359 }
360 normvec[v] = Teuchos::ScalarTraits<MagnitudeType>::squareroot(value);
361 }
362 }
363
364 // Copy the vectors in A to a set of vectors in *this. The numvecs vectors in
365 // A are copied to a subset of vectors in *this indicated by the indices given
366 // in index.
367 // FIXME: not so clear what the size of A and index.size() are...
369 const std::vector<int> &index)
370 {
371 MyMultiVec* MyA;
372 MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A));
373 TEUCHOS_ASSERT(MyA != NULL);
374
375 assert (A.GetNumberVecs() >= (int)index.size());
376 assert (A.GetGlobalLength() == Length_);
377
378 for (unsigned int v = 0 ; v < index.size() ; ++v) {
379 for (int i = 0 ; i < Length_ ; ++i) {
380 (*this)(i, index[v]) = (*MyA)(i, v);
381 }
382 }
383 }
384
385 // Fill the vectors in *this with random numbers.
386 void MvRandom ()
387 {
388 Teuchos::SerialDenseMatrix<int,ScalarType> R( Length_, NumberVecs_ );
389 Teuchos::randomSyncedMatrix( R );
390 for (int v = 0 ; v < NumberVecs_ ; ++v) {
391 for (int i = 0 ; i < Length_ ; ++i) {
392 (*this)(i, v) = R(i, v);
393 }
394 }
395 }
396
397 // Replace each element of the vectors in *this with alpha.
398 void MvInit (const ScalarType alpha)
399 {
400 for (int v = 0 ; v < NumberVecs_ ; ++v) {
401 for (int i = 0 ; i < Length_ ; ++i) {
402 (*this)(i, v) = alpha;
403 }
404 }
405 }
406
407 void MvPrint (std::ostream &os) const
408 {
409 os << "Object MyMultiVec" << std::endl;
410 os << "Number of rows = " << Length_ << std::endl;
411 os << "Number of vecs = " << NumberVecs_ << std::endl;
412
413 for (int i = 0 ; i < Length_ ; ++i)
414 {
415 for (int v = 0 ; v < NumberVecs_ ; ++v)
416 os << (*this)(i, v) << " ";
417 os << std::endl;
418 }
419 }
420
421 inline ScalarType& operator()(const int i, const int j)
422 {
423 if (j < 0 || j >= NumberVecs_) throw(-1);
424 if (i < 0 || i >= Length_) throw(-2);
425 //
426 return(data_[j][i]);
427 }
428
429 inline const ScalarType& operator()(const int i, const int j) const
430 {
431 if (j < 0 || j >= NumberVecs_) throw(-1);
432 if (i < 0 || i >= Length_) throw(-2);
433 //
434 return(data_[j][i]);
435 }
436
437 ScalarType* operator[](int v)
438 {
439 if (v < 0 || v >= NumberVecs_) throw(-1);
440 return(data_[v]);
441 }
442
443 ScalarType* operator[](int v) const
444 {
445 return(data_[v]);
446 }
447
448private:
449 void Check()
450 {
451 if (Length_ <= 0)
452 throw("Length must be positive");
453
454 if (NumberVecs_ <= 0)
455 throw("Number of vectors must be positive");
456 }
457
459 const ptrdiff_t Length_;
461 const int NumberVecs_;
463 std::vector<ScalarType*> data_;
465 std::vector<bool> ownership_;
466
467};
468
469
470#endif // MY_MULTIVECTOR_HPP
Belos header file which uses auto-configuration information to include necessary C++ headers.
Interface for multivectors used by Belos' linear solvers.
Interface for multivectors used by Belos' linear solvers.
virtual ptrdiff_t GetGlobalLength() const =0
The number of rows in the multivector.
virtual int GetNumberVecs() const =0
The number of vectors (i.e., columns) in the multivector.
Simple example of a user's defined Belos::MultiVec class.
void MvScale(const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
ScalarType & operator()(const int i, const int j)
ScalarType * operator[](int v) const
MyMultiVec(const MyMultiVec &rhs)
Copy constructor, performs a deep copy.
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
MyMultiVec(const ptrdiff_t Length, const int NumberVecs)
Constructor for a NumberVecs vectors of length Length.
MyMultiVec(const ptrdiff_t Length, const std::vector< ScalarType * > &rhs)
Constructor with already allocated memory.
ScalarType * operator[](int v)
void MvPrint(std::ostream &os) const
Print *this multivector to the os output stream.
const ptrdiff_t Length_
Length of the vectors.
void SetBlock(const Belos::MultiVec< ScalarType > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
const ScalarType & operator()(const int i, const int j) const
MyMultiVec * Clone(const int NumberVecs) const
Returns a clone of the current std::vector.
std::vector< ScalarType * > data_
Pointers to the storage of the vectors.
void MvRandom()
Fill all the vectors in *this with random numbers.
MyMultiVec * CloneCopy(const std::vector< int > &index) const
Returns a clone copy of specified vectors.
const int NumberVecs_
Number of multi-vectors.
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
void MvNorm(std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, Belos::NormType type=Belos::TwoNorm) const
Compute the norm of each vector in *this.
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).
void MvInit(const ScalarType alpha)
Replace each element of the vectors in *this with alpha.
std::vector< bool > ownership_
If true, then this object owns the vectors and must free them in dtor.
~MyMultiVec()
Destructor.
MyMultiVec * CloneCopy() const
Create a new MultiVec and copy contents of *this into it (deep copy).
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.
MyMultiVec * CloneViewNonConst(const std::vector< int > &index)
Returns a view of current std::vector (shallow copy)
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).
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.
void MvScale(const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
const MyMultiVec * CloneView(const std::vector< int > &index) const
Returns a view of current std::vector (shallow copy), const version.
NormType
The type of vector norm to compute.