Belos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MyOperator.hpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ************************************************************************
4 //
5 // Belos: Block Linear Solvers Package
6 // Copyright 2004 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 */
43 
44 #ifndef MY_OPERATOR_HPP
45 #define MY_OPERATOR_HPP
46 
47 #include "BelosConfigDefs.hpp"
48 #include "BelosOperator.hpp"
49 #include "MyMultiVec.hpp"
50 
52 
64 template <class ScalarType>
65 class MyOperator : public Belos::Operator<ScalarType>
66 {
67 
68 public:
69 
70  /* Constructs a square matrix with \c NumRows rows and columns.
71  * The matrix is tridiagonal, and the computational stencil is
72  * [-1, 2, -1]
73  */
74  MyOperator(const int NumRows) :
75  NumRows_(NumRows)
76  {
77  l_ = -1.0;
78  d_ = 2.0;
79  u_ = -1.0;
80  }
81 
82  // Constructor for tridiagonal matrix.
83  MyOperator(const int NumRows, std::vector<ScalarType> ldu) :
84  NumRows_(NumRows)
85  {
86  l_ = ldu[0];
87  d_ = ldu[1];
88  u_ = ldu[2];
89  }
90 
91  // Constructor for a diagonal matrix with variable entries.
92  MyOperator(std::vector<ScalarType> diag) :
93  NumRows_(diag.size())
94  {
95  diag_.resize(diag.size());
96  for(unsigned int i=0; i<diag_.size(); ++i)
97  diag_[i] = diag[i];
98  }
99 
102  {}
103 
107  Belos::ETrans trans = Belos::NOTRANS) const
108  {
109  const MyMultiVec<ScalarType>* MyX;
110  MyX = dynamic_cast<const MyMultiVec<ScalarType>*>(&X);
111  assert (MyX != 0);
112 
114  MyY = dynamic_cast<MyMultiVec<ScalarType>*>(&Y);
115  assert (MyY != 0);
116 
117  assert (X.GetNumberVecs() == Y.GetNumberVecs());
118  assert (X.GetGlobalLength() == Y.GetGlobalLength());
119 
120  if (diag_.size() == 0)
121  {
122  // This is a tridiagonal matrix
123  for (int v = 0 ; v < X.GetNumberVecs() ; ++v)
124  {
125  for (ptrdiff_t i = 0 ; i < X.GetGlobalLength() ; ++i)
126  {
127  if (i == 0)
128  {
129  (*MyY)[v][i] = (d_ * (*MyX)[v][i] + u_ * (*MyX)[v][i + 1]);
130  }
131  else if (i == X.GetGlobalLength() - 1)
132  {
133  (*MyY)[v][i] = (d_ * (*MyX)[v][i] + l_ * (*MyX)[v][i-1]);
134  }
135  else
136  {
137  (*MyY)[v][i] = (d_ * (*MyX)[v][i] + l_ * (*MyX)[v][i-1] + u_ * (*MyX)[v][i+1]);
138  }
139  }
140  }
141  }
142  else
143  {
144  // This is a diagonal matrix
145  for (int v = 0 ; v < X.GetNumberVecs() ; ++v)
146  {
147  for (ptrdiff_t i = 0 ; i < X.GetGlobalLength() ; ++i)
148  {
149  (*MyY)[v][i] = diag_[i] * (*MyX)[v][i];
150  }
151  }
152  }
153  }
154 
155 private:
157  int NumRows_;
159  ScalarType l_, d_, u_;
161  std::vector<ScalarType> diag_;
162 };
163 
164 #endif //MY_OPERATOR_HPP
ScalarType u_
Definition: MyOperator.hpp:159
MyOperator(const int NumRows, std::vector< ScalarType > ldu)
Definition: MyOperator.hpp:83
virtual int GetNumberVecs() const =0
The number of vectors (i.e., columns) in the multivector.
void Apply(const Belos::MultiVec< ScalarType > &X, Belos::MultiVec< ScalarType > &Y, Belos::ETrans trans=Belos::NOTRANS) const
Applies the tridiagonal or diagonal matrix to a multivector.
Definition: MyOperator.hpp:105
virtual ptrdiff_t GetGlobalLength() const =0
The number of rows in the multivector.
ETrans
Whether to apply the (conjugate) transpose of an operator.
Definition: BelosTypes.hpp:81
MyOperator(std::vector< ScalarType > diag)
Definition: MyOperator.hpp:92
Simple example of a user&#39;s defined Belos::MultiVec class.
Definition: MyMultiVec.hpp:65
std::vector< ScalarType > diag_
Elements on diagonal (for variable-diagonal case).
Definition: MyOperator.hpp:161
Alternative run-time polymorphic interface for operators.
~MyOperator()
Dtor.
Definition: MyOperator.hpp:101
ScalarType d_
Definition: MyOperator.hpp:159
int NumRows_
Number of rows and columns.
Definition: MyOperator.hpp:157
ScalarType l_
Elements on subdiagonal, diagonal, and superdiagonal.
Definition: MyOperator.hpp:159
Alternative run-time polymorphic interface for operators.
Simple example of a user&#39;s defined Belos::Operator class.
Definition: MyOperator.hpp:65
Interface for multivectors used by Belos&#39; linear solvers.
Belos header file which uses auto-configuration information to include necessary C++ headers...
MyOperator(const int NumRows)
Definition: MyOperator.hpp:74