Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
MyOperator.hpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Stratimikos: Thyra-based strategies for linear solvers
6 // Copyright (2006) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 */
43 
44 #ifndef MY_OPERATOR_HPP
45 #define MY_OPERATOR_HPP
46 
47 #include "Tpetra_Operator.hpp"
48 #include "Tpetra_Vector.hpp"
49 #include "Tpetra_VectorSpace.hpp"
50 #include "Teuchos_BLAS.hpp"
51 
53 
62 template <class OrdinalType, class ScalarType>
63 class MyOperator : public Tpetra::Operator<OrdinalType,ScalarType>
64 {
65 
66 public:
67 
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)
72  : _vs(vs), _nr(nrows), _nnz(nnz), _cptr(nrows+1), _rind(nnz), _vals(nnz)
73  {
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());
77  }
78 
81  {
82  }
83 
86 
88  Tpetra::VectorSpace<OrdinalType,ScalarType> const& getDomainDist() const { return _vs; };
89 
91  Tpetra::VectorSpace<OrdinalType,ScalarType> const& getRangeDist() const { return _vs; };
92 
94  void apply(Tpetra::Vector<OrdinalType,ScalarType> const& x,
95  Tpetra::Vector<OrdinalType, ScalarType> & y,
96  bool transpose=false) const
97  {
98  // Get the indexes of the rows on this processor
99  const int numMyElements = _vs.getNumMyEntries();
100  const std::vector<int> &myGlobalElements = _vs.elementSpace().getMyGlobalElements();
101 
102  // Initialize output std::vector to zero.
103  y.setAllToScalar( Teuchos::ScalarTraits<ScalarType>::zero() );
104 
105  assert (x.getNumGlobalEntries() == y.getNumGlobalEntries());
106  assert (x.getNumGlobalEntries() == y.getNumGlobalEntries());
107 
108  // Apply operator
109  int IA1, IA2, ri;
110  int i,j;
111  for (int myRow = 0; myRow < numMyElements; ++myRow ) {
112 
113  // For each row this processor owns, compute the value of A*x[myRow]
114  const int rowIndex = myGlobalElements[myRow];
115  for (j=0; j<_nr; j++) {
116  IA1 = _cptr[j]-1;
117  IA2 = _cptr[j+1]-1;
118  for (i=IA1; i<IA2; i++) {
119  ri = _rind[i]-1;
120  if (ri == rowIndex)
121  y[rowIndex] += _vals[i]*x[j];
122  } // end for (i= ...)
123  } // end for (j= ...)
124  } // end for (myRow= ...)
125 
126  };
127 
129 
130 private:
131 
132  typedef typename std::vector<ScalarType>::iterator STIter;
133  typedef std::vector<int>::iterator IntIter;
134 
136  Tpetra::VectorSpace<OrdinalType,ScalarType> _vs;
137 
139  int _nr, _nnz;
141  std::vector<int> _cptr;
143  std::vector<int> _rind;
145  std::vector<ScalarType> _vals;
146 };
147 
148 #endif //MY_OPERATOR_HPP
int _nr
Number of rows and columns.
Definition: MyOperator.hpp:139
std::vector< ScalarType >::iterator STIter
Definition: MyOperator.hpp:126
std::vector< int >::iterator IntIter
Definition: MyOperator.hpp:133
MyOperator(const Tpetra::VectorSpace< OrdinalType, ScalarType > &vs, const int nrows, const int *colptr, const int nnz, const int *rowin, const ScalarType *vals)
Constructor.
Definition: MyOperator.hpp:69
Tpetra::VectorSpace< OrdinalType, ScalarType > const & getRangeDist() const
Returns the VectorSpace associated with the range of this linear operator.
Definition: MyOperator.hpp:91
std::vector< int > _rind
Row indices.
Definition: MyOperator.hpp:143
Tpetra::VectorSpace< OrdinalType, ScalarType > const & getDomainDist() const
Returns the VectorSpace associated with the domain of this linear operator.
Definition: MyOperator.hpp:88
Tpetra::VectorSpace< OrdinalType, ScalarType > _vs
Tpetra std::vector space.
Definition: MyOperator.hpp:136
std::vector< int > _cptr
Column pointers.
Definition: MyOperator.hpp:141
Simple example of a user&#39;s defined Tpetra::Operator class.
Definition: MyOperator.hpp:63
std::vector< ScalarType > _vals
Values.
Definition: MyOperator.hpp:145
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.
Definition: MyOperator.hpp:94
~MyOperator()
Deconstructor.
Definition: MyOperator.hpp:80