MyOperator.hpp

Go to the documentation of this file.
00001 #ifndef MY_OPERATOR_HPP
00002 #define MY_OPERATOR_HPP
00003 
00004 #include "BelosConfigDefs.hpp"
00005 #include "BelosOperator.hpp"
00006 #include "MyMultiVec.hpp"
00007 
00009 
00021 template <class ScalarType>
00022 class MyOperator : public Belos::Operator<ScalarType>
00023 {
00024 
00025 public:
00026   
00027   /* Constructs a square matrix with \c NumRows rows and columns.
00028    * The matrix is tridiagonal, and the computational stencil is 
00029    * [-1, 2, -1]
00030    */
00031   MyOperator(const int NumRows) :
00032     NumRows_(NumRows)
00033   {
00034     l_ = -1.0;
00035     d_ =  2.0;
00036     u_ = -1.0;
00037   }
00038 
00039   // Constructor for tridiagonal matrix.
00040   MyOperator(const int NumRows, std::vector<ScalarType> ldu) :
00041     NumRows_(NumRows)
00042   {
00043     l_ = ldu[0];
00044     d_ = ldu[1];
00045     u_ = ldu[2];    
00046   }
00047 
00048   // Constructor for a diagonal matrix with variable entries.
00049   MyOperator(std::vector<ScalarType> diag) :
00050     NumRows_(diag.size())
00051   {
00052     diag_.resize(diag.size());
00053     for(unsigned int i=0; i<diag_.size(); ++i)
00054       diag_[i] = diag[i];
00055   }
00056 
00058   ~MyOperator()
00059   {}
00060   
00062   void Apply(const Belos::MultiVec<ScalarType>& X, 
00063              Belos::MultiVec<ScalarType>& Y,
00064              Belos::ETrans trans = Belos::NOTRANS) const
00065   {
00066     const MyMultiVec<ScalarType>* MyX;
00067     MyX = dynamic_cast<const MyMultiVec<ScalarType>*>(&X); 
00068     assert (MyX != 0);
00069     
00070     MyMultiVec<ScalarType>* MyY;
00071     MyY = dynamic_cast<MyMultiVec<ScalarType>*>(&Y); 
00072     assert (MyY != 0);
00073     
00074     assert (X.GetNumberVecs() == Y.GetNumberVecs());
00075     assert (X.GetVecLength() == Y.GetVecLength());
00076    
00077     if (diag_.size() == 0)
00078     {
00079       // This is a tridiagonal matrix
00080       for (int v = 0 ; v < X.GetNumberVecs() ; ++v)
00081       {
00082         for (int i = 0 ; i < X.GetVecLength() ; ++i)
00083         {
00084           if (i == 0)
00085           {
00086             (*MyY)[v][i] = (d_ * (*MyX)[v][i] + u_ * (*MyX)[v][i + 1]);
00087           }
00088           else if (i == X.GetVecLength() - 1)
00089           {
00090             (*MyY)[v][i] = (d_ * (*MyX)[v][i] + l_ * (*MyX)[v][i-1]);
00091           }
00092           else
00093           {
00094             (*MyY)[v][i] = (d_ * (*MyX)[v][i] + l_ * (*MyX)[v][i-1] + u_ * (*MyX)[v][i+1]);
00095           }
00096         }
00097       }
00098     } 
00099     else
00100     {
00101       // This is a diagonal matrix
00102       for (int v = 0 ; v < X.GetNumberVecs() ; ++v)
00103       {
00104         for (int i = 0 ; i < X.GetVecLength() ; ++i)
00105         {
00106           (*MyY)[v][i] = diag_[i] * (*MyX)[v][i];
00107         }
00108       }      
00109     }
00110   }
00111   
00112 private:
00114   int NumRows_;
00116   ScalarType l_, d_, u_;
00118   std::vector<ScalarType> diag_;
00119 };
00120 
00121 #endif //MY_OPERATOR_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