Line data Source code
1 : /*
2 : * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef MATRIX_H_
34 : #define MATRIX_H_
35 :
36 : #include <cstddef>
37 : #include <ostream>
38 : #include <iomanip>
39 : #include "../../types.h"
40 : #include "math_vector.h"
41 : #include "common/assert.h"
42 :
43 : namespace ug
44 : {
45 :
46 : /**
47 : * \defgroup math_matrix Matrix
48 : * \ingroup ugbase_math
49 : * \{
50 : */
51 :
52 : template <std::size_t N, std::size_t M, typename T = number> class MathMatrix;
53 :
54 : /**
55 : * \class MathMatrix
56 : *
57 : * \brief A class for fixed size, dense matrices.
58 : *
59 : * A static memory NxM matrix
60 : */
61 : template <std::size_t N, std::size_t M, typename T>
62 : class MathMatrix
63 : {
64 : public:
65 : typedef T value_type;
66 : typedef std::size_t size_type;
67 : static const std::size_t RowSize = N;
68 : static const std::size_t ColSize = M;
69 :
70 : public:
71 137376 : MathMatrix() {}
72 0 : MathMatrix(const MathMatrix& v) {assign(v);}
73 :
74 : /**
75 : * \brief Assigns the elements of the given matrix to this one.
76 : *
77 : * \param v The matrix to be assigned.
78 : * \return A reference to this matrix.
79 : */
80 : MathMatrix& operator= (const MathMatrix& v)
81 : {
82 : assign(v);
83 : return *this;
84 : }
85 :
86 : /**
87 : * \brief Adds a matrix to 'this' one: \f$ A_{this} \leftarrow A_{this} + B\f$.
88 : *
89 : * \param B The matrix to be added.
90 : * \return A reference to this matrix.
91 : */
92 : MathMatrix& operator+= (const MathMatrix& B)
93 : {
94 0 : for(std::size_t i = 0; i < N; ++i){
95 0 : for(std::size_t j = 0; j < M; ++j){
96 0 : m_data[i][j] += B(i,j);
97 : }
98 : }
99 : return *this;
100 : }
101 :
102 : /**
103 : * \brief Subtracts a matrix from 'this' one: \f$ A_{this} \leftarrow A_{this} - B\f$.
104 : *
105 : * \param B The matrix to be subtracted.
106 : * \return A reference to this matrix.
107 : */
108 : MathMatrix& operator-= (const MathMatrix& B)
109 : {
110 0 : for(std::size_t i = 0; i < N; ++i){
111 0 : for(std::size_t j = 0; j < M; ++j){
112 0 : m_data[i][j] -= B(i,j);
113 : }
114 : }
115 : return *this;
116 : }
117 :
118 : /**
119 : * \brief Assigns the given value to all elements of the matrix.
120 : *
121 : * \param val The value to be assigned to the matrix.
122 : * \return A reference to this matrix.
123 : */
124 : MathMatrix& operator= (const value_type& val)
125 : {
126 0 : for(std::size_t i = 0; i < N; ++i){
127 0 : for(std::size_t j = 0; j < M; ++j){
128 0 : m_data[i][j] = val;
129 : }
130 : }
131 : return *this;
132 : }
133 :
134 : /**
135 : * \brief Adds the given value to all elements of the matrix.
136 : *
137 : * \param val The value to be added.
138 : * \return A reference to this matrix.
139 : */
140 : MathMatrix& operator+= (const value_type& val)
141 : {
142 : for(std::size_t i = 0; i < N; ++i){
143 : for(std::size_t j = 0; j < M; ++j){
144 : m_data[i][j] += val;
145 : }
146 : }
147 : return *this;
148 : }
149 :
150 : /**
151 : * \brief Subtracts the given value from all elements of the matrix.
152 : *
153 : * \param val The value to be subtracted.
154 : * \return A reference to this matrix.
155 : */
156 : MathMatrix& operator-= (const value_type& val)
157 : {
158 : for(std::size_t i = 0; i < N; ++i){
159 : for(std::size_t j = 0; j < M; ++j){
160 : m_data[i][j] -= val;
161 : }
162 : }
163 : return *this;
164 : }
165 :
166 : /**
167 : * \brief Divides all elements of the matrix by the given value.
168 : *
169 : * \param val The divisor.
170 : * \return A reference to this matrix.
171 : */
172 : MathMatrix& operator/= (const value_type& val)
173 : {
174 : for(std::size_t i = 0; i < N; ++i){
175 : for(std::size_t j = 0; j < M; ++j){
176 : m_data[i][j] /= val;
177 : }
178 : }
179 : return *this;
180 : }
181 :
182 : /**
183 : * \brief Multiplies all elements of the matrix with the given value.
184 : *
185 : * \param val The factor.
186 : * \return A reference to this matrix.
187 : */
188 : MathMatrix& operator*= (const value_type& val)
189 : {
190 : for(std::size_t i = 0; i < N; ++i){
191 : for(std::size_t j = 0; j < M; ++j){
192 : m_data[i][j] *= val;
193 : }
194 : }
195 : return *this;
196 : }
197 :
198 : /**
199 : * \brief Multiplies the matrix element-wise with another matrix and sums up the entries.
200 : *
201 : * \param v The Matrix.
202 : * \return A scalar value of the element-wise summed up products
203 : */
204 : value_type operator* (const MathMatrix& v) const
205 : {
206 : value_type res = 0.0;
207 0 : for(std::size_t i = 0; i < N; ++i){
208 0 : for(std::size_t j = 0; j < M; ++j){
209 0 : res += m_data[i][j] * v.m_data[i][j];
210 : }
211 : }
212 : return res;
213 : }
214 :
215 : inline std::size_t num_rows() const {return N;}
216 : inline std::size_t num_cols() const {return M;}
217 :
218 : inline value_type* operator[](std::size_t index){
219 : UG_ASSERT(index < N, "Invalid index");
220 0 : return m_data[index];
221 : }
222 : inline const value_type* operator[](std::size_t index) const{
223 : UG_ASSERT(index < N, "Invalid index");
224 0 : return m_data[index];
225 : }
226 :
227 : inline value_type& operator() (std::size_t row, std::size_t col){
228 : UG_ASSERT(row < N && col < M, "Accessing "<<N<<"x"<<M<<"Matrix at entry ("<<row<<","<<col<<")");
229 : return m_data[row][col];
230 : }
231 : inline const value_type& operator() (std::size_t row, std::size_t col) const{
232 : UG_ASSERT(row < N && col < M, "Accessing "<<N<<"x"<<M<<"Matrix at entry ("<<row<<","<<col<<")");
233 : return m_data[row][col];
234 : }
235 :
236 : inline value_type& entry(std::size_t row, std::size_t col) {return (*this)(row,col);}
237 : inline const value_type& entry(std::size_t row, std::size_t col) const {return (*this)(row,col);}
238 :
239 : inline void assign(const MathVector<N, value_type>& vec, const std::size_t row) {
240 : UG_ASSERT(vec.Size == N, "Wrong vector size");
241 : for(std::size_t j = 0; j < N; j++)
242 : m_data[row][j] = vec[j];
243 : }
244 :
245 : protected:
246 : value_type m_data[N][M];
247 :
248 : inline void assign(const MathMatrix& v)
249 : {
250 0 : for(std::size_t i = 0; i < N; ++i){
251 0 : for(std::size_t j = 0; j < M; ++j){
252 0 : m_data[i][j] = v.m_data[i][j] ;
253 : }
254 : }
255 : }
256 :
257 : };
258 :
259 : // this are explicit instantiations to avoid compiler errors on XL
260 : // {
261 : template <typename T> class MathMatrix<0,0,T>{
262 : public:
263 : typedef T value_type;
264 : typedef std::size_t size_type;
265 : };
266 : template <std::size_t N, typename T> class MathMatrix<N,0,T>{
267 : public:
268 : typedef T value_type;
269 : typedef std::size_t size_type;
270 : };
271 : template <std::size_t N, typename T> class MathMatrix<0,N,T>{
272 : public:
273 : typedef T value_type;
274 : typedef std::size_t size_type;
275 : };
276 : // }
277 :
278 : /// Print MathMatrix<N,M> to standard output
279 : template <std::size_t N, std::size_t M>
280 0 : std::ostream& operator<< (std::ostream& outStream, const ug::MathMatrix<N,M>& m)
281 : {
282 0 : for(std::size_t i = 0; i < N; ++i)
283 : {
284 0 : for(std::size_t j = 0; j < M; ++j)
285 : {
286 0 : outStream << "[" << i << "][" << j << "]: " << std::scientific << std::setprecision(8) << std::setw(15) << m.entry(i, j) << std::endl;
287 : }
288 : }
289 0 : return outStream;
290 : }
291 :
292 : std::ostream& operator<< (std::ostream& outStream, const ug::MathMatrix<2,2>& m);
293 : std::ostream& operator<< (std::ostream& outStream, const ug::MathMatrix<2,3>& m);
294 : std::ostream& operator<< (std::ostream& outStream, const ug::MathMatrix<3,2>& m);
295 : std::ostream& operator<< (std::ostream& outStream, const ug::MathMatrix<3,3>& m);
296 :
297 : // end group math_matrix
298 : /// \}
299 :
300 : } //end of namespace: lgmath
301 :
302 :
303 : #endif /* MathMatrix_H_ */
|
