MyMultiVec.hpp

Go to the documentation of this file.
00001 #ifndef MY_MULTIVECTOR_HPP
00002 #define MY_MULTIVECTOR_HPP
00003 
00004 #include "BelosConfigDefs.hpp"
00005 #include "BelosMultiVec.hpp"
00006 
00008 
00018 template <class ScalarType>
00019 class MyMultiVec : public Belos::MultiVec<ScalarType>
00020 {
00021 public:
00022 
00024   MyMultiVec(const int Length, const int NumberVecs) :
00025     Length_(Length),
00026     NumberVecs_(NumberVecs)
00027   {
00028     Check();
00029 
00030     data_.resize(NumberVecs);
00031     ownership_.resize(NumberVecs);
00032     
00033     // Allocates memory to store the vectors.
00034     for (int v = 0 ; v < NumberVecs_ ; ++v)
00035     {
00036       data_[v] = new ScalarType[Length];
00037       ownership_[v] = true;
00038     }
00039 
00040     // Initializes all elements to zero.
00041     MvInit(0.0);
00042   }
00043   
00045   MyMultiVec(const int Length, const std::vector<ScalarType*>& rhs) :
00046     Length_(Length),
00047     NumberVecs_(rhs.size())
00048   {
00049     Check();
00050 
00051     data_.resize(NumberVecs_);
00052     ownership_.resize(NumberVecs_);
00053     
00054     // Copies pointers from input array, set ownership so that
00055     // this memory is not free'd in the destructor
00056     for (int v = 0 ; v < NumberVecs_ ; ++v)
00057     {
00058       data_[v] = rhs[v];
00059       ownership_[v] = false;
00060     }
00061   }
00062   
00064   MyMultiVec(const MyMultiVec& rhs) :
00065     Length_(rhs.GetVecLength()),
00066     NumberVecs_(rhs.NumberVecs_)
00067   {
00068     Check();
00069 
00070     data_.resize(NumberVecs_);
00071     ownership_.resize(NumberVecs_);
00072     
00073     for (int v = 0 ; v < NumberVecs_ ; ++v)
00074     {
00075       data_[v] = new ScalarType[Length_];
00076       ownership_[v] = true;
00077     }
00078     
00079     for (int v = 0 ; v < NumberVecs_ ; ++v)
00080     {
00081       for (int i = 0 ; i < Length_ ; ++i)
00082         (*this)(i, v) = rhs(i, v);
00083     }
00084   }
00085   
00087   ~MyMultiVec()
00088   {
00089     for (int v = 0 ; v < NumberVecs_ ; ++v)
00090       if (ownership_[v]) 
00091         delete[] data_[v];
00092   }
00093   
00095   MyMultiVec* Clone(const int NumberVecs) const
00096   {
00097     // FIXME
00098     MyMultiVec* tmp = new MyMultiVec(Length_, NumberVecs);
00099     
00100     //   for (int v = 0 ; v < NumberVecs ; ++v)
00101     //         for (int i = 0 ; i < Length_ ; ++i)
00102     //           (*tmp)(i, v) = (*this)(i, v);
00103     
00104     return(tmp);
00105   }
00106   
00107   // Returns a clone of the corrent multi-std::vector.
00108   MyMultiVec* CloneCopy() const
00109   {
00110     return(new MyMultiVec(*this));
00111   }
00112   
00114   MyMultiVec* CloneCopy(const std::vector< int > &index) const
00115   {
00116     int size = index.size();
00117     MyMultiVec* tmp = new MyMultiVec(Length_, size);
00118     
00119     for (unsigned int v = 0 ; v < index.size() ; ++v) {
00120       for (int i = 0 ; i < Length_ ; ++i) {
00121         (*tmp)(i, v) = (*this)(i, index[v]);
00122       }
00123     }
00124     
00125     return(tmp);
00126   }
00127   
00129   MyMultiVec* CloneViewNonConst(const std::vector< int > &index) 
00130   {
00131     int size = index.size();
00132     std::vector<ScalarType*> values(size);
00133     
00134     for (unsigned int v = 0 ; v < index.size() ; ++v)
00135       values[v] = data_[index[v]];
00136     
00137     return(new MyMultiVec(Length_, values));
00138   }
00139   
00141   const MyMultiVec* CloneView(const std::vector< int > &index) const
00142   {
00143     int size = index.size();
00144     std::vector<ScalarType*> values(size);
00145     
00146     for (unsigned int v = 0 ; v < index.size() ; ++v)
00147       values[v] = data_[index[v]];
00148     
00149     return(new MyMultiVec(Length_, values));
00150   }
00151   
00152   int GetVecLength () const
00153   {
00154     return(Length_);
00155   }
00156   
00157   int GetNumberVecs () const
00158   {
00159     return(NumberVecs_);
00160   }
00161   
00162   // Update *this with alpha * A * B + beta * (*this). 
00163   void MvTimesMatAddMv (const ScalarType alpha, const Belos::MultiVec<ScalarType> &A, 
00164                         const Teuchos::SerialDenseMatrix<int, ScalarType> &B, 
00165                         const ScalarType beta)
00166   {
00167     
00168     assert (Length_ == A.GetVecLength());
00169     assert (B.numRows() == A.GetNumberVecs());
00170     assert (B.numCols() <= NumberVecs_);
00171 
00172     MyMultiVec* MyA;
00173     MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A)); 
00174     assert(MyA!=NULL);
00175 
00176     if ((*this)[0] == (*MyA)[0]) {
00177       // If this == A, then need additional storage ...
00178       // This situation is a bit peculiar but it may be required by
00179       // certain algorithms.
00180       
00181       std::vector<ScalarType> tmp(NumberVecs_);
00182 
00183       for (int i = 0 ; i < Length_ ; ++i) {
00184         for (int v = 0; v < A.GetNumberVecs() ; ++v) {
00185           tmp[v] = (*MyA)(i, v);
00186         }
00187 
00188         for (int v = 0 ; v < B.numCols() ; ++v) {
00189           (*this)(i, v) *= beta; 
00190           ScalarType res = Teuchos::ScalarTraits<ScalarType>::zero();
00191 
00192           for (int j = 0 ; j < A.GetNumberVecs() ; ++j) {
00193             res +=  tmp[j] * B(j, v);
00194           }
00195 
00196           (*this)(i, v) += alpha * res;
00197         }
00198       }
00199     }
00200     else {
00201       for (int i = 0 ; i < Length_ ; ++i) {
00202         for (int v = 0 ; v < B.numCols() ; ++v) {
00203           (*this)(i, v) *= beta; 
00204           ScalarType res = 0.0;
00205           for (int j = 0 ; j < A.GetNumberVecs() ; ++j) {
00206             res +=  (*MyA)(i, j) * B(j, v);
00207           }
00208 
00209           (*this)(i, v) += alpha * res;
00210         }
00211       }
00212     }
00213   }
00214   
00215   // Replace *this with alpha * A + beta * B. 
00216   void MvAddMv (const ScalarType alpha, const Belos::MultiVec<ScalarType>& A, 
00217                 const ScalarType beta,  const Belos::MultiVec<ScalarType>& B)
00218   {
00219     MyMultiVec* MyA;
00220     MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A)); 
00221     assert (MyA != 0);
00222     
00223     MyMultiVec* MyB;
00224     MyB = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(B)); 
00225     assert (MyB != 0);
00226     
00227     assert (NumberVecs_ == A.GetNumberVecs());
00228     assert (NumberVecs_ == B.GetNumberVecs());
00229     
00230     assert (Length_ == A.GetVecLength());
00231     assert (Length_ == B.GetVecLength());
00232     
00233     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00234       for (int i = 0 ; i < Length_ ; ++i) {
00235         (*this)(i, v) = alpha * (*MyA)(i, v) + beta * (*MyB)(i, v);
00236       }
00237     }
00238   }
00239   
00240   // Replace each element of the vectors in *this with (*this)*alpha.
00241   void MvScale (const ScalarType alpha)
00242   {
00243     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00244       for (int i = 0 ; i < Length_ ; ++i) {
00245         (*this)(i, v) *= alpha;
00246       }
00247     }
00248   }
00249   
00250   // Replace each element of the vectors in *this with (*this)*alpha[i].
00251   void MvScale (const std::vector<ScalarType>& alpha)
00252   {
00253     assert((int)alpha.size() >= NumberVecs_);
00254     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00255       for (int i = 0 ; i < Length_ ; ++i) {
00256         (*this)(i, v) *= alpha[v];
00257       }
00258     }
00259   }
00260   
00261   // Compute a dense matrix B through the matrix-matrix multiply alpha * A^H * (*this). 
00262   void MvTransMv (const ScalarType alpha, const Belos::MultiVec<ScalarType>& A, 
00263                   Teuchos::SerialDenseMatrix< int, ScalarType >& B) const
00264   {
00265     MyMultiVec* MyA;
00266     MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A)); 
00267     assert (MyA != 0);
00268     
00269     assert (A.GetVecLength() == Length_);
00270     assert (NumberVecs_ <= B.numCols());
00271     assert (A.GetNumberVecs() <= B.numRows());
00272     
00273       for (int v = 0 ; v < A.GetNumberVecs() ; ++v) {
00274         for (int w = 0 ; w < NumberVecs_ ; ++w) {
00275           ScalarType value = 0.0;
00276           for (int i = 0 ; i < Length_ ; ++i) {
00277             value += Teuchos::ScalarTraits<ScalarType>::conjugate((*MyA)(i, v)) * (*this)(i, w);
00278           }
00279           B(v, w) = alpha * value;
00280         }
00281       }
00282   }
00283   
00284   
00285   // 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. 
00286   void MvDot (const Belos::MultiVec<ScalarType>& A, std::vector<ScalarType>& b) const
00287   {
00288     MyMultiVec* MyA;
00289     MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A)); 
00290     assert (MyA != 0);
00291     
00292     assert (NumberVecs_ <= (int)b.size());
00293     assert (NumberVecs_ == A.GetNumberVecs());
00294     assert (Length_ == A.GetVecLength());
00295     
00296       for (int v = 0 ; v < NumberVecs_ ; ++v) {
00297         ScalarType value = 0.0;
00298         for (int i = 0 ; i < Length_ ; ++i) {
00299           value += (*this)(i, v) * Teuchos::ScalarTraits<ScalarType>::conjugate((*MyA)(i, v));
00300         }
00301         b[v] = value;
00302       }
00303   }  
00304   
00305   
00306   void MvNorm ( std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec,
00307     Belos::NormType type = Belos::TwoNorm ) const 
00308   {
00309     assert (NumberVecs_ <= (int)normvec.size());
00310 
00311     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00312     
00313     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00314       MagnitudeType value = Teuchos::ScalarTraits<MagnitudeType>::zero();
00315       for (int i = 0 ; i < Length_ ; ++i) {
00316         MagnitudeType val = Teuchos::ScalarTraits<ScalarType>::magnitude((*this)(i, v));
00317         value += val * val;
00318       }
00319       normvec[v] = Teuchos::ScalarTraits<MagnitudeType>::squareroot(value);
00320     }
00321   }
00322   
00323   // Copy the vectors in A to a set of vectors in *this. The numvecs vectors in 
00324   // A are copied to a subset of vectors in *this indicated by the indices given 
00325   // in index.
00326   // FIXME: not so clear what the size of A and index.size() are...
00327   void SetBlock (const Belos::MultiVec<ScalarType>& A, 
00328                  const std::vector<int> &index)
00329   {
00330     MyMultiVec* MyA;
00331     MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A)); 
00332     assert (MyA != 0);
00333     
00334     assert (A.GetNumberVecs() >= (int)index.size());
00335     assert (A.GetVecLength() == Length_);
00336     
00337     for (unsigned int v = 0 ; v < index.size() ; ++v) {
00338       for (int i = 0 ; i < Length_ ; ++i) {
00339         (*this)(i, index[v])  = (*MyA)(i, v);
00340       }
00341     }
00342   }
00343   
00344   // Fill the vectors in *this with random numbers.
00345   void MvRandom ()
00346   {
00347     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00348       for (int i = 0 ; i < Length_ ; ++i) {
00349         (*this)(i, v) = Teuchos::ScalarTraits<ScalarType>::random();
00350       }
00351     }
00352   }
00353   
00354   // Replace each element of the vectors in *this with alpha.
00355   void MvInit (const ScalarType alpha)
00356   {
00357     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00358       for (int i = 0 ; i < Length_ ; ++i) {
00359         (*this)(i, v) = alpha;
00360       }
00361     }
00362   }
00363   
00364   void MvPrint (std::ostream &os) const
00365   {
00366     os << "Object MyMultiVec" << std::endl;
00367     os << "Number of rows = " << Length_ << std::endl;
00368     os << "Number of vecs = " << NumberVecs_ << std::endl;
00369     
00370     for (int i = 0 ; i < Length_ ; ++i)
00371       {
00372         for (int v = 0 ; v < NumberVecs_ ; ++v)
00373           os << (*this)(i, v) << " ";
00374         os << std::endl;
00375       }
00376   }
00377   
00378   inline ScalarType& operator()(const int i, const int j)
00379   {
00380     if (j < 0 || j >= NumberVecs_) throw(-1);
00381     if (i < 0 || i >= Length_) throw(-2);
00382     // 
00383     return(data_[j][i]);
00384   }
00385   
00386   inline const ScalarType& operator()(const int i, const int j) const
00387   {
00388     if (j < 0 || j >= NumberVecs_) throw(-1);
00389     if (i < 0 || i >= Length_) throw(-2);
00390     // 
00391     return(data_[j][i]);
00392   }
00393   
00394   ScalarType* operator[](int v)
00395   {
00396     if (v < 0 || v >= NumberVecs_) throw(-1);
00397     return(data_[v]);
00398   }
00399   
00400   ScalarType* operator[](int v) const
00401   {
00402     return(data_[v]);
00403   }
00404   
00405 private:
00406   void Check()
00407   {
00408     if (Length_ <= 0)
00409       throw("Length must be positive");
00410 
00411     if (NumberVecs_ <= 0)
00412       throw("Number of vectors must be positive");
00413   }
00414 
00416   const int Length_;
00418   const int NumberVecs_;
00420   std::vector<ScalarType*> data_;
00422   std::vector<bool> ownership_;
00423 
00424 };
00425 
00426 
00427 #endif // MY_MULTIVECTOR_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:15 2011 for Belos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3