44 #ifndef MY_OPERATOR_HPP 45 #define MY_OPERATOR_HPP 47 #include "Tpetra_Operator.hpp" 48 #include "Tpetra_Vector.hpp" 49 #include "Tpetra_VectorSpace.hpp" 50 #include "Teuchos_BLAS.hpp" 62 template <
class OrdinalType,
class ScalarType>
63 class MyOperator :
public Tpetra::Operator<OrdinalType,ScalarType>
69 MyOperator(
const Tpetra::VectorSpace<OrdinalType,ScalarType>& vs,
70 const int nrows,
const int *colptr,
71 const int nnz,
const int *rowin,
const ScalarType *vals)
74 std::copy<const int*,IntIter>(colptr,colptr+nrows+1,
_cptr.begin());
75 std::copy<const int*,IntIter>(rowin,rowin+nnz,
_rind.begin());
76 std::copy<const ScalarType*,STIter>(vals,vals+nnz,
_vals.begin());
88 Tpetra::VectorSpace<OrdinalType,ScalarType>
const&
getDomainDist()
const {
return _vs; };
91 Tpetra::VectorSpace<OrdinalType,ScalarType>
const&
getRangeDist()
const {
return _vs; };
94 void apply(Tpetra::Vector<OrdinalType,ScalarType>
const& x,
95 Tpetra::Vector<OrdinalType, ScalarType> & y,
96 bool transpose=
false)
const 99 const int numMyElements =
_vs.getNumMyEntries();
100 const std::vector<int> &myGlobalElements =
_vs.elementSpace().getMyGlobalElements();
103 y.setAllToScalar( Teuchos::ScalarTraits<ScalarType>::zero() );
105 assert (x.getNumGlobalEntries() == y.getNumGlobalEntries());
106 assert (x.getNumGlobalEntries() == y.getNumGlobalEntries());
111 for (
int myRow = 0; myRow < numMyElements; ++myRow ) {
114 const int rowIndex = myGlobalElements[myRow];
115 for (j=0; j<
_nr; j++) {
118 for (i=IA1; i<IA2; i++) {
121 y[rowIndex] +=
_vals[i]*x[j];
132 typedef typename std::vector<ScalarType>::iterator
STIter;
136 Tpetra::VectorSpace<OrdinalType,ScalarType>
_vs;
148 #endif //MY_OPERATOR_HPP int _nr
Number of rows and columns.
std::vector< ScalarType >::iterator STIter
std::vector< int >::iterator IntIter
MyOperator(const Tpetra::VectorSpace< OrdinalType, ScalarType > &vs, const int nrows, const int *colptr, const int nnz, const int *rowin, const ScalarType *vals)
Constructor.
Tpetra::VectorSpace< OrdinalType, ScalarType > const & getRangeDist() const
Returns the VectorSpace associated with the range of this linear operator.
std::vector< int > _rind
Row indices.
Tpetra::VectorSpace< OrdinalType, ScalarType > const & getDomainDist() const
Returns the VectorSpace associated with the domain of this linear operator.
Tpetra::VectorSpace< OrdinalType, ScalarType > _vs
Tpetra std::vector space.
std::vector< int > _cptr
Column pointers.
Simple example of a user's defined Tpetra::Operator class.
std::vector< ScalarType > _vals
Values.
void apply(Tpetra::Vector< OrdinalType, ScalarType > const &x, Tpetra::Vector< OrdinalType, ScalarType > &y, bool transpose=false) const
Computes the matrix-std::vector multiplication y = Ax.
~MyOperator()
Deconstructor.