MyBetterOperator.hpp

Go to the documentation of this file.
00001 #ifndef MY_BETTER_OPERATOR_HPP
00002 #define MY_BETTER_OPERATOR_HPP
00003 
00004 #include "BelosConfigDefs.hpp"
00005 #include "BelosOperator.hpp"
00006 #include "MyMultiVec.hpp"
00007 #include "Teuchos_BLAS.hpp"
00008 
00010 
00022 template <class ScalarType>
00023 class MyBetterOperator : public Belos::Operator<ScalarType> 
00024 {
00025 
00026 public:
00027 
00029   MyBetterOperator(const int nrows, const int *colptr,
00030                    const int nnz, const int *rowin, const ScalarType *vals)
00031   : _nr(nrows), _nnz(nnz), _cptr(nrows+1), _rind(nnz), _vals(nnz)
00032   {
00033     std::copy<const int*,IntIter>(colptr,colptr+nrows+1,_cptr.begin());
00034     std::copy<const int*,IntIter>(rowin,rowin+nnz,_rind.begin());
00035     std::copy<const ScalarType*,STIter>(vals,vals+nnz,_vals.begin());
00036   }
00037 
00039   ~MyBetterOperator()
00040   { }
00041 
00043   void Apply(const Belos::MultiVec<ScalarType>& X, 
00044              Belos::MultiVec<ScalarType>& Y,
00045              Belos::ETrans trans = Belos::NOTRANS) const
00046   {
00047     const MyMultiVec<ScalarType>* MyX;
00048     MyX = dynamic_cast<const MyMultiVec<ScalarType>*>(&X); 
00049     assert (MyX != 0);
00050     
00051     MyMultiVec<ScalarType>* MyY;
00052     MyY = dynamic_cast<MyMultiVec<ScalarType>*>(&Y); 
00053     assert (MyY != 0);
00054 
00055     // Initialize output std::vector to zero.
00056     MyY->MvInit( Teuchos::ScalarTraits<ScalarType>::zero() );
00057     
00058     assert (X.GetNumberVecs() == Y.GetNumberVecs());
00059     assert (X.GetVecLength() == Y.GetVecLength());
00060     
00061     int nv = X.GetNumberVecs();
00062 
00063     // Apply operator
00064     int IA1, IA2, ri;
00065     ScalarType aval;
00066     int i,j,v;
00067     for (j=0; j<_nr; j++) {
00068       IA1 = _cptr[j]-1;
00069       IA2 = _cptr[j+1]-1;
00070       for (i=IA1; i<IA2; i++) {
00071   ri = _rind[i]-1;
00072         aval = _vals[i];
00073         for (v=0; v<nv; v++) {
00074           (*MyY)[v][ri] += aval*(*MyX)[v][j];
00075         }
00076       }
00077     }
00078   }
00079 
00080   void Print( std::ostream& os ) {
00081     for (int j=0; j<_nr; j++) {
00082       int IA1 = _cptr[j]-1;
00083       int IA2 = _cptr[j+1]-1;
00084       for (int i=IA1; i<IA2; i++) {
00085         os << "("<<_rind[i]-1<<","<<j<<")\t"<<_vals[i]<< std::endl;
00086       }
00087     } 
00088   }
00089 
00090   private:
00091   typedef typename std::vector<ScalarType>::iterator STIter;
00092   typedef std::vector<int>::iterator        IntIter;
00094   int _nr, _nnz;
00096   std::vector<int> _cptr;
00098   std::vector<int> _rind;
00100   std::vector<ScalarType> _vals;
00101 };
00102 
00103 #endif //MY_BETTER_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