44 #ifndef MY_MULTIVECTOR_HPP 45 #define MY_MULTIVECTOR_HPP 61 template <
class ScalarType>
67 MyMultiVec (
const ptrdiff_t Length,
const int NumberVecs) :
73 data_.resize(NumberVecs);
78 data_[v] =
new ScalarType[Length];
87 MyMultiVec(
const ptrdiff_t Length,
const std::vector<ScalarType*>& rhs) :
122 for (
int i = 0 ; i <
Length_ ; ++i)
123 (*
this)(i, v) = rhs(i, v);
157 int size = index.size();
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]);
172 int size = index.size();
173 std::vector<ScalarType*> values(size);
175 for (
size_t v = 0 ; v < index.size() ; ++v)
176 values[v] =
data_[index[v]];
184 int size = index.size();
185 std::vector<ScalarType*> values (size);
187 for (
size_t v = 0; v < index.size (); ++v) {
188 values[v] =
data_[index[v]];
207 const ScalarType beta)
209 assert (
Length_ ==
A.GetGlobalLength());
210 assert (
B.numRows() ==
A.GetNumberVecs());
216 if ((*
this)[0] == (*MyA)[0]) {
223 for (
int i = 0 ; i <
Length_ ; ++i) {
224 for (
int v = 0; v <
A.GetNumberVecs() ; ++v) {
225 tmp[v] = (*MyA)(i, v);
228 for (
int v = 0 ; v <
B.numCols() ; ++v) {
229 (*this)(i, v) *= beta;
232 for (
int j = 0 ; j <
A.GetNumberVecs() ; ++j) {
233 res += tmp[j] *
B(j, v);
236 (*this)(i, v) += alpha * res;
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);
249 (*this)(i, v) += alpha * res;
268 assert (
Length_ ==
A.GetGlobalLength());
269 assert (
Length_ ==
B.GetGlobalLength());
272 for (
int i = 0 ; i <
Length_ ; ++i) {
273 (*this)(i, v) = alpha * (*MyA)(i, v) + beta * (*MyB)(i, v);
282 for (
int i = 0 ; i <
Length_ ; ++i) {
283 (*this)(i, v) *= alpha;
289 void MvScale (
const std::vector<ScalarType>& alpha)
293 for (
int i = 0 ; i <
Length_ ; ++i) {
294 (*this)(i, v) *= alpha[v];
307 assert (
A.GetGlobalLength() ==
Length_);
309 assert (
A.GetNumberVecs() <=
B.numRows());
311 for (
int v = 0 ; v <
A.GetNumberVecs() ; ++v) {
313 ScalarType value = 0.0;
314 for (
int i = 0 ; i <
Length_ ; ++i) {
317 B(v, w) = alpha * value;
332 assert (
Length_ ==
A.GetGlobalLength());
335 ScalarType value = 0.0;
336 for (
int i = 0 ; i <
Length_ ; ++i) {
353 for (
int i = 0 ; i <
Length_ ; ++i) {
366 const std::vector<int> &index)
372 assert (
A.GetNumberVecs() >= (int)index.size());
373 assert (
A.GetGlobalLength() ==
Length_);
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);
386 for (
int i = 0 ; i <
Length_ ; ++i) {
396 for (
int i = 0 ; i <
Length_ ; ++i) {
397 (*this)(i, v) = alpha;
404 os <<
"Object MyMultiVec" << std::endl;
405 os <<
"Number of rows = " <<
Length_ << std::endl;
406 os <<
"Number of vecs = " <<
NumberVecs_ << std::endl;
408 for (
int i = 0 ; i <
Length_ ; ++i)
411 os << (*
this)(i, v) <<
" ";
419 if (i < 0 || i >=
Length_)
throw(-2);
424 inline const ScalarType&
operator()(
const int i,
const int j)
const 427 if (i < 0 || i >=
Length_)
throw(-2);
447 throw(
"Length must be positive");
450 throw(
"Number of vectors must be positive");
465 #endif // MY_MULTIVECTOR_HPP
MyMultiVec(const MyMultiVec &rhs)
Copy constructor, performs a deep copy.
std::vector< ScalarType * > data_
Pointers to the storage of the vectors.
MyMultiVec * CloneCopy() const
Create a new MultiVec and copy contents of *this into it (deep copy).
const ScalarType & operator()(const int i, const int j) const
ScalarType & operator()(const int i, const int j)
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 MvPrint(std::ostream &os) const
Print *this multivector to the os output stream.
void MvScale(const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
MyMultiVec * CloneViewNonConst(const std::vector< int > &index)
Returns a view of current std::vector (shallow copy)
const MyMultiVec * CloneView(const std::vector< int > &index) const
Returns a view of current std::vector (shallow copy), const version.
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).
Simple example of a user's defined Belos::MultiVec class.
int GetNumberVecs() const
The number of vectors (i.e., columns) in the multivector.
MyMultiVec(const ptrdiff_t Length, const std::vector< ScalarType *> &rhs)
Constructor with already allocated memory.
MyMultiVec * CloneCopy(const std::vector< int > &index) const
Returns a clone copy of specified vectors.
ScalarType * operator[](int v)
static magnitudeType magnitude(T a)
NormType
The type of vector norm to compute.
void MvRandom()
Fill all the vectors in *this with random numbers.
ptrdiff_t GetGlobalLength() const
The number of rows in the multivector.
const int NumberVecs_
Number of multi-vectors.
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).
Interface for multivectors used by Belos' 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.
void MvScale(const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
const ptrdiff_t Length_
Length of the vectors.
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.
void SetBlock(const Belos::MultiVec< ScalarType > &A, const std::vector< int > &index)
Copy the vectors in A to a set of vectors in *this.
std::vector< bool > ownership_
If true, then this object owns the vectors and must free them in dtor.
MyMultiVec(const ptrdiff_t Length, const int NumberVecs)
Constructor for a NumberVecs vectors of length Length.
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.
ScalarType * operator[](int v) const
void MvInit(const ScalarType alpha)
Replace each element of the vectors in *this with alpha.
Interface for multivectors used by Belos' linear solvers.