Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
MyOperator.hpp
Go to the documentation of this file.
00001 #ifndef MY_OPERATOR_HPP
00002 #define MY_OPERATOR_HPP
00003 
00004 #include "Tpetra_Operator.hpp"
00005 #include "Tpetra_Vector.hpp"
00006 #include "Tpetra_VectorSpace.hpp"
00007 #include "Teuchos_BLAS.hpp"
00008 
00010 
00017 template <class OrdinalType, class ScalarType>
00018 class MyOperator : public Tpetra::Operator<OrdinalType,ScalarType> 
00019 {
00020 
00021 public:
00022 
00024   MyOperator(const Tpetra::VectorSpace<OrdinalType,ScalarType>& vs,
00025        const int nrows, const int *colptr,
00026        const int nnz, const int *rowin, const ScalarType *vals)
00027     : _vs(vs), _nr(nrows), _nnz(nnz), _cptr(nrows+1), _rind(nnz), _vals(nnz)
00028   {
00029     std::copy<const int*,IntIter>(colptr,colptr+nrows+1,_cptr.begin());
00030     std::copy<const int*,IntIter>(rowin,rowin+nnz,_rind.begin());
00031     std::copy<const ScalarType*,STIter>(vals,vals+nnz,_vals.begin());
00032   }
00033   
00035   ~MyOperator()
00036   {
00037   }
00038 
00041   
00043   Tpetra::VectorSpace<OrdinalType,ScalarType> const& getDomainDist() const { return _vs; };
00044   
00046   Tpetra::VectorSpace<OrdinalType,ScalarType> const& getRangeDist() const { return _vs; };
00047   
00049   void apply(Tpetra::Vector<OrdinalType,ScalarType> const& x, 
00050        Tpetra::Vector<OrdinalType, ScalarType> & y, 
00051        bool transpose=false) const 
00052   {
00053     // Get the indexes of the rows on this processor
00054     const int numMyElements = _vs.getNumMyEntries();
00055     const std::vector<int> &myGlobalElements = _vs.elementSpace().getMyGlobalElements();
00056     
00057     // Initialize output std::vector to zero.
00058     y.setAllToScalar( Teuchos::ScalarTraits<ScalarType>::zero() );
00059 
00060     assert (x.getNumGlobalEntries() == y.getNumGlobalEntries());
00061     assert (x.getNumGlobalEntries() == y.getNumGlobalEntries());
00062     
00063     // Apply operator
00064     int IA1, IA2, ri;
00065     int i,j;
00066     for (int myRow = 0; myRow < numMyElements; ++myRow ) {
00067 
00068       // For each row this processor owns, compute the value of A*x[myRow]
00069       const int rowIndex = myGlobalElements[myRow];
00070       for (j=0; j<_nr; j++) {   
00071   IA1 = _cptr[j]-1;
00072   IA2 = _cptr[j+1]-1;
00073   for (i=IA1; i<IA2; i++) {
00074     ri = _rind[i]-1;
00075     if (ri == rowIndex) 
00076       y[rowIndex] += _vals[i]*x[j];
00077   } // end for (i= ...)
00078       } // end for (j= ...)
00079     } // end for (myRow= ...)
00080     
00081   };
00082   
00084   
00085 private:
00086 
00087   typedef typename std::vector<ScalarType>::iterator STIter;
00088   typedef std::vector<int>::iterator        IntIter;
00089 
00091   Tpetra::VectorSpace<OrdinalType,ScalarType> _vs;
00092 
00094   int _nr, _nnz;
00096   std::vector<int> _cptr;
00098   std::vector<int> _rind;
00100   std::vector<ScalarType> _vals;
00101 };
00102 
00103 #endif //MY_OPERATOR_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines