Belos Package Browser (Single Doxygen Collection) Development
MyMultiVec.hpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 //
00005 //                 Belos: Block Linear Solvers Package
00006 //                  Copyright 2004 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #ifndef MY_MULTIVECTOR_HPP
00045 #define MY_MULTIVECTOR_HPP
00046 
00047 #include "BelosConfigDefs.hpp"
00048 #include "BelosMultiVec.hpp"
00049 
00051 
00061 template <class ScalarType>
00062 class MyMultiVec : public Belos::MultiVec<ScalarType>
00063 {
00064 public:
00065 
00067   MyMultiVec(const int Length, const int NumberVecs) :
00068     Length_(Length),
00069     NumberVecs_(NumberVecs)
00070   {
00071     Check();
00072 
00073     data_.resize(NumberVecs);
00074     ownership_.resize(NumberVecs);
00075     
00076     // Allocates memory to store the vectors.
00077     for (int v = 0 ; v < NumberVecs_ ; ++v)
00078     {
00079       data_[v] = new ScalarType[Length];
00080       ownership_[v] = true;
00081     }
00082 
00083     // Initializes all elements to zero.
00084     MvInit(0.0);
00085   }
00086   
00088   MyMultiVec(const int Length, const std::vector<ScalarType*>& rhs) :
00089     Length_(Length),
00090     NumberVecs_(rhs.size())
00091   {
00092     Check();
00093 
00094     data_.resize(NumberVecs_);
00095     ownership_.resize(NumberVecs_);
00096     
00097     // Copies pointers from input array, set ownership so that
00098     // this memory is not free'd in the destructor
00099     for (int v = 0 ; v < NumberVecs_ ; ++v)
00100     {
00101       data_[v] = rhs[v];
00102       ownership_[v] = false;
00103     }
00104   }
00105   
00107   MyMultiVec(const MyMultiVec& rhs) :
00108     Length_(rhs.GetVecLength()),
00109     NumberVecs_(rhs.NumberVecs_)
00110   {
00111     Check();
00112 
00113     data_.resize(NumberVecs_);
00114     ownership_.resize(NumberVecs_);
00115     
00116     for (int v = 0 ; v < NumberVecs_ ; ++v)
00117     {
00118       data_[v] = new ScalarType[Length_];
00119       ownership_[v] = true;
00120     }
00121     
00122     for (int v = 0 ; v < NumberVecs_ ; ++v)
00123     {
00124       for (int i = 0 ; i < Length_ ; ++i)
00125         (*this)(i, v) = rhs(i, v);
00126     }
00127   }
00128   
00130   ~MyMultiVec()
00131   {
00132     for (int v = 0 ; v < NumberVecs_ ; ++v)
00133       if (ownership_[v]) 
00134         delete[] data_[v];
00135   }
00136   
00138   MyMultiVec* Clone(const int NumberVecs) const
00139   {
00140     // FIXME
00141     MyMultiVec* tmp = new MyMultiVec(Length_, NumberVecs);
00142     
00143     //   for (int v = 0 ; v < NumberVecs ; ++v)
00144     //         for (int i = 0 ; i < Length_ ; ++i)
00145     //           (*tmp)(i, v) = (*this)(i, v);
00146     
00147     return(tmp);
00148   }
00149   
00150   // Returns a clone of the corrent multi-std::vector.
00151   MyMultiVec* CloneCopy() const
00152   {
00153     return(new MyMultiVec(*this));
00154   }
00155   
00157   MyMultiVec* CloneCopy(const std::vector< int > &index) const
00158   {
00159     int size = index.size();
00160     MyMultiVec* tmp = new MyMultiVec(Length_, size);
00161     
00162     for (unsigned int v = 0 ; v < index.size() ; ++v) {
00163       for (int i = 0 ; i < Length_ ; ++i) {
00164         (*tmp)(i, v) = (*this)(i, index[v]);
00165       }
00166     }
00167     
00168     return(tmp);
00169   }
00170   
00172   MyMultiVec* CloneViewNonConst(const std::vector< int > &index) 
00173   {
00174     int size = index.size();
00175     std::vector<ScalarType*> values(size);
00176     
00177     for (unsigned int v = 0 ; v < index.size() ; ++v)
00178       values[v] = data_[index[v]];
00179     
00180     return(new MyMultiVec(Length_, values));
00181   }
00182   
00184   const MyMultiVec* CloneView(const std::vector< int > &index) const
00185   {
00186     int size = index.size();
00187     std::vector<ScalarType*> values(size);
00188     
00189     for (unsigned int v = 0 ; v < index.size() ; ++v)
00190       values[v] = data_[index[v]];
00191     
00192     return(new MyMultiVec(Length_, values));
00193   }
00194   
00195   int GetVecLength () const
00196   {
00197     return(Length_);
00198   }
00199   
00200   int GetNumberVecs () const
00201   {
00202     return(NumberVecs_);
00203   }
00204   
00205   // Update *this with alpha * A * B + beta * (*this). 
00206   void MvTimesMatAddMv (const ScalarType alpha, const Belos::MultiVec<ScalarType> &A, 
00207                         const Teuchos::SerialDenseMatrix<int, ScalarType> &B, 
00208                         const ScalarType beta)
00209   {
00210     
00211     assert (Length_ == A.GetVecLength());
00212     assert (B.numRows() == A.GetNumberVecs());
00213     assert (B.numCols() <= NumberVecs_);
00214 
00215     MyMultiVec* MyA;
00216     MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A)); 
00217     assert(MyA!=NULL);
00218 
00219     if ((*this)[0] == (*MyA)[0]) {
00220       // If this == A, then need additional storage ...
00221       // This situation is a bit peculiar but it may be required by
00222       // certain algorithms.
00223       
00224       std::vector<ScalarType> tmp(NumberVecs_);
00225 
00226       for (int i = 0 ; i < Length_ ; ++i) {
00227         for (int v = 0; v < A.GetNumberVecs() ; ++v) {
00228           tmp[v] = (*MyA)(i, v);
00229         }
00230 
00231         for (int v = 0 ; v < B.numCols() ; ++v) {
00232           (*this)(i, v) *= beta; 
00233           ScalarType res = Teuchos::ScalarTraits<ScalarType>::zero();
00234 
00235           for (int j = 0 ; j < A.GetNumberVecs() ; ++j) {
00236             res +=  tmp[j] * B(j, v);
00237           }
00238 
00239           (*this)(i, v) += alpha * res;
00240         }
00241       }
00242     }
00243     else {
00244       for (int i = 0 ; i < Length_ ; ++i) {
00245         for (int v = 0 ; v < B.numCols() ; ++v) {
00246           (*this)(i, v) *= beta; 
00247           ScalarType res = 0.0;
00248           for (int j = 0 ; j < A.GetNumberVecs() ; ++j) {
00249             res +=  (*MyA)(i, j) * B(j, v);
00250           }
00251 
00252           (*this)(i, v) += alpha * res;
00253         }
00254       }
00255     }
00256   }
00257   
00258   // Replace *this with alpha * A + beta * B. 
00259   void MvAddMv (const ScalarType alpha, const Belos::MultiVec<ScalarType>& A, 
00260                 const ScalarType beta,  const Belos::MultiVec<ScalarType>& B)
00261   {
00262     MyMultiVec* MyA;
00263     MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A)); 
00264     assert (MyA != 0);
00265     
00266     MyMultiVec* MyB;
00267     MyB = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(B)); 
00268     assert (MyB != 0);
00269     
00270     assert (NumberVecs_ == A.GetNumberVecs());
00271     assert (NumberVecs_ == B.GetNumberVecs());
00272     
00273     assert (Length_ == A.GetVecLength());
00274     assert (Length_ == B.GetVecLength());
00275     
00276     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00277       for (int i = 0 ; i < Length_ ; ++i) {
00278         (*this)(i, v) = alpha * (*MyA)(i, v) + beta * (*MyB)(i, v);
00279       }
00280     }
00281   }
00282   
00283   // Replace each element of the vectors in *this with (*this)*alpha.
00284   void MvScale (const ScalarType alpha)
00285   {
00286     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00287       for (int i = 0 ; i < Length_ ; ++i) {
00288         (*this)(i, v) *= alpha;
00289       }
00290     }
00291   }
00292   
00293   // Replace each element of the vectors in *this with (*this)*alpha[i].
00294   void MvScale (const std::vector<ScalarType>& alpha)
00295   {
00296     assert((int)alpha.size() >= NumberVecs_);
00297     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00298       for (int i = 0 ; i < Length_ ; ++i) {
00299         (*this)(i, v) *= alpha[v];
00300       }
00301     }
00302   }
00303   
00304   // Compute a dense matrix B through the matrix-matrix multiply alpha * A^H * (*this). 
00305   void MvTransMv (const ScalarType alpha, const Belos::MultiVec<ScalarType>& A, 
00306                   Teuchos::SerialDenseMatrix< int, ScalarType >& B) const
00307   {
00308     MyMultiVec* MyA;
00309     MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A)); 
00310     assert (MyA != 0);
00311     
00312     assert (A.GetVecLength() == Length_);
00313     assert (NumberVecs_ <= B.numCols());
00314     assert (A.GetNumberVecs() <= B.numRows());
00315     
00316       for (int v = 0 ; v < A.GetNumberVecs() ; ++v) {
00317         for (int w = 0 ; w < NumberVecs_ ; ++w) {
00318           ScalarType value = 0.0;
00319           for (int i = 0 ; i < Length_ ; ++i) {
00320             value += Teuchos::ScalarTraits<ScalarType>::conjugate((*MyA)(i, v)) * (*this)(i, w);
00321           }
00322           B(v, w) = alpha * value;
00323         }
00324       }
00325   }
00326   
00327   
00328   // 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. 
00329   void MvDot (const Belos::MultiVec<ScalarType>& A, std::vector<ScalarType>& b) const
00330   {
00331     MyMultiVec* MyA;
00332     MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A)); 
00333     assert (MyA != 0);
00334     
00335     assert (NumberVecs_ <= (int)b.size());
00336     assert (NumberVecs_ == A.GetNumberVecs());
00337     assert (Length_ == A.GetVecLength());
00338     
00339       for (int v = 0 ; v < NumberVecs_ ; ++v) {
00340         ScalarType value = 0.0;
00341         for (int i = 0 ; i < Length_ ; ++i) {
00342           value += (*this)(i, v) * Teuchos::ScalarTraits<ScalarType>::conjugate((*MyA)(i, v));
00343         }
00344         b[v] = value;
00345       }
00346   }  
00347   
00348   
00349   void MvNorm ( std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec,
00350     Belos::NormType type = Belos::TwoNorm ) const 
00351   {
00352     assert (NumberVecs_ <= (int)normvec.size());
00353 
00354     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00355     
00356     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00357       MagnitudeType value = Teuchos::ScalarTraits<MagnitudeType>::zero();
00358       for (int i = 0 ; i < Length_ ; ++i) {
00359         MagnitudeType val = Teuchos::ScalarTraits<ScalarType>::magnitude((*this)(i, v));
00360         value += val * val;
00361       }
00362       normvec[v] = Teuchos::ScalarTraits<MagnitudeType>::squareroot(value);
00363     }
00364   }
00365   
00366   // Copy the vectors in A to a set of vectors in *this. The numvecs vectors in 
00367   // A are copied to a subset of vectors in *this indicated by the indices given 
00368   // in index.
00369   // FIXME: not so clear what the size of A and index.size() are...
00370   void SetBlock (const Belos::MultiVec<ScalarType>& A, 
00371                  const std::vector<int> &index)
00372   {
00373     MyMultiVec* MyA;
00374     MyA = dynamic_cast<MyMultiVec*>(&const_cast<Belos::MultiVec<ScalarType> &>(A)); 
00375     assert (MyA != 0);
00376     
00377     assert (A.GetNumberVecs() >= (int)index.size());
00378     assert (A.GetVecLength() == Length_);
00379     
00380     for (unsigned int v = 0 ; v < index.size() ; ++v) {
00381       for (int i = 0 ; i < Length_ ; ++i) {
00382         (*this)(i, index[v])  = (*MyA)(i, v);
00383       }
00384     }
00385   }
00386   
00387   // Fill the vectors in *this with random numbers.
00388   void MvRandom ()
00389   {
00390     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00391       for (int i = 0 ; i < Length_ ; ++i) {
00392         (*this)(i, v) = Teuchos::ScalarTraits<ScalarType>::random();
00393       }
00394     }
00395   }
00396   
00397   // Replace each element of the vectors in *this with alpha.
00398   void MvInit (const ScalarType alpha)
00399   {
00400     for (int v = 0 ; v < NumberVecs_ ; ++v) {
00401       for (int i = 0 ; i < Length_ ; ++i) {
00402         (*this)(i, v) = alpha;
00403       }
00404     }
00405   }
00406   
00407   void MvPrint (std::ostream &os) const
00408   {
00409     os << "Object MyMultiVec" << std::endl;
00410     os << "Number of rows = " << Length_ << std::endl;
00411     os << "Number of vecs = " << NumberVecs_ << std::endl;
00412     
00413     for (int i = 0 ; i < Length_ ; ++i)
00414       {
00415         for (int v = 0 ; v < NumberVecs_ ; ++v)
00416           os << (*this)(i, v) << " ";
00417         os << std::endl;
00418       }
00419   }
00420   
00421   inline ScalarType& operator()(const int i, const int j)
00422   {
00423     if (j < 0 || j >= NumberVecs_) throw(-1);
00424     if (i < 0 || i >= Length_) throw(-2);
00425     // 
00426     return(data_[j][i]);
00427   }
00428   
00429   inline const ScalarType& operator()(const int i, const int j) const
00430   {
00431     if (j < 0 || j >= NumberVecs_) throw(-1);
00432     if (i < 0 || i >= Length_) throw(-2);
00433     // 
00434     return(data_[j][i]);
00435   }
00436   
00437   ScalarType* operator[](int v)
00438   {
00439     if (v < 0 || v >= NumberVecs_) throw(-1);
00440     return(data_[v]);
00441   }
00442   
00443   ScalarType* operator[](int v) const
00444   {
00445     return(data_[v]);
00446   }
00447   
00448 private:
00449   void Check()
00450   {
00451     if (Length_ <= 0)
00452       throw("Length must be positive");
00453 
00454     if (NumberVecs_ <= 0)
00455       throw("Number of vectors must be positive");
00456   }
00457 
00459   const int Length_;
00461   const int NumberVecs_;
00463   std::vector<ScalarType*> data_;
00465   std::vector<bool> ownership_;
00466 
00467 };
00468 
00469 
00470 #endif // MY_MULTIVECTOR_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines