Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
64template <class ScalarType>
65class MyOperator : public Belos::Operator<ScalarType>
66{
67
68public:
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
155private:
159 ScalarType l_, d_, u_;
161 std::vector<ScalarType> diag_;
162};
163
164#endif //MY_OPERATOR_HPP
Belos header file which uses auto-configuration information to include necessary C++ headers.
Alternative run-time polymorphic interface for operators.
Interface for multivectors used by Belos' linear solvers.
virtual ptrdiff_t GetGlobalLength() const =0
The number of rows in the multivector.
virtual int GetNumberVecs() const =0
The number of vectors (i.e., columns) in the multivector.
Alternative run-time polymorphic interface for operators.
Simple example of a user's defined Belos::MultiVec class.
Definition: MyMultiVec.hpp:66
Simple example of a user's defined Belos::Operator class.
Definition: MyOperator.hpp:66
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
MyOperator(const int NumRows, std::vector< ScalarType > ldu)
Definition: MyOperator.hpp:83
~MyOperator()
Dtor.
Definition: MyOperator.hpp:101
ScalarType l_
Elements on subdiagonal, diagonal, and superdiagonal.
Definition: MyOperator.hpp:159
ScalarType u_
Definition: MyOperator.hpp:159
MyOperator(const int NumRows)
Definition: MyOperator.hpp:74
ScalarType d_
Definition: MyOperator.hpp:159
std::vector< ScalarType > diag_
Elements on diagonal (for variable-diagonal case).
Definition: MyOperator.hpp:161
int NumRows_
Number of rows and columns.
Definition: MyOperator.hpp:157
MyOperator(std::vector< ScalarType > diag)
Definition: MyOperator.hpp:92
ETrans
Whether to apply the (conjugate) transpose of an operator.
Definition: BelosTypes.hpp:81
@ NOTRANS
Definition: BelosTypes.hpp:81