ug4
Loading...
Searching...
No Matches
math_symmetric_matrix.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 * Author: Christian Wehner
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 SYM_MATRIX_H_
34#define SYM_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
43namespace ug
44{
45
52//template <std::size_t N, std::size_t M, typename T = number> class MathSymmetricMatrix;
53
61template <std::size_t N, typename T = number>
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 = N;
69
70 public:
73
81 {
82 assign(v);
83 return *this;
84 }
85
93 {
94 for(std::size_t i = 0; i < m_size; ++i)
95 {
96 m_data[i] += B[i];
97 }
98 return *this;
99 }
100
108 {
109 for(std::size_t i = 0; i < m_size; ++i)
110 {
111 m_data[i] -= B[i];
112 }
113 return *this;
114 }
115
123 {
124 for(std::size_t i = 0; i < m_size; ++i)
125 {
126 m_data[i] = val;
127 }
128 return *this;
129 }
130
138 {
139 for(std::size_t i = 0; i < m_size; ++i)
140 {
141 m_data[i] += val;
142 }
143 return *this;
144 }
145
153 {
154 for(std::size_t i = 0; i < m_size; ++i)
155 {
156 m_data[i] -= val;
157 }
158 return *this;
159 }
160
168 {
169 for(std::size_t i = 0; i < m_size; ++i)
170 {
171 m_data[i] /= val;
172 }
173 return *this;
174 }
175
183 {
184 for(std::size_t i = 0; i < m_size; ++i)
185 {
186 m_data[i] *= val;
187 }
188 return *this;
189 }
190
198 {
199 value_type res = 0.0;
200 for(std::size_t i = 0; i < N; ++i)
201 {
202 for(std::size_t j = 0; j < N; ++j)
203 {
204 res += this->entry(i,j) * v(i,j);
205 }
206 }
207 return res;
208 }
209
210 inline std::size_t num_rows() const {return N;}
211 inline std::size_t num_cols() const {return N;}
212
213 inline value_type& entry(std::size_t row, std::size_t col)
214 {
215 UG_ASSERT(row < N && col < N, "Accessing "<<N<<"x"<<N<<"Matrix at entry ("<<row<<","<<col<<")");
216 if (row<col) return m_data[row * N - (row - 1) * row / 2 + col - row];
217 else return m_data[col * N - (col - 1) * col / 2 + row - col];
218 }
219
220 inline const value_type& entry(std::size_t row, std::size_t col) const
221 {
222 UG_ASSERT(row < N && col < N, "Accessing "<<N<<"x"<<N<<"Matrix at entry ("<<row<<","<<col<<")");
223 if (row<col) return m_data[row * N - (row - 1) * row / 2 + col - row];
224 else return m_data[col * N - (col - 1) * col / 2 + row - col];
225 }
226
227 inline value_type& operator[](std::size_t index) {UG_ASSERT(index < m_size, "Invalid index"); return m_data[index];}
228 inline const value_type& operator[](std::size_t index) const {UG_ASSERT(index < m_size, "Invalid index"); return m_data[index];}
229
230 inline value_type& operator() (std::size_t row, std::size_t col)
231 {
232 UG_ASSERT(row < N && col < N, "Accessing "<<N<<"x"<<N<<"Matrix at entry ("<<row<<","<<col<<")");
233 if (row<col) return m_data[row * N - (row - 1) * row / 2 + col - row];
234 else return m_data[col * N - (col - 1) * col / 2 + row - col];
235 }
236
237 inline const value_type& operator() (std::size_t row, std::size_t col) const
238 {
239 UG_ASSERT(row < N && col < N, "Accessing "<<N<<"x"<<N<<"Matrix at entry ("<<row<<","<<col<<")");
240 if (row<col) return m_data[row * N - (row - 1) * row / 2 + col - row];
241 else return m_data[col * N - (col - 1) * col / 2 + row - col];
242 }
243
244 // frobenius norm |A|_{F} = \sqrt{ 2 \sum_{i,j} A_{i,j} }
246 value_type norm=m_data[0]*m_data[0];
247 for (size_t i=1;i<m_size-1;i++){
248 norm += 2*m_data[i]*m_data[i];
249 }
250 norm+=m_data[m_size-1]*m_data[m_size-1];
251 if (N>2)
252 for (size_t diag=1;diag<N-1;diag++)
253 norm -= m_data[diag * N - (diag - 1) * diag / 2]*m_data[diag * N - (diag - 1) * diag / 2];
254 return sqrt(2.0*norm);
255 }
256
257 // scale matrix by frobenius norm
259 number norm=fnorm();
260 for (size_t i=0;i<m_size;i++) m_data[i]*=norm;
261 }
262
263 protected:
264 static const size_t m_size = (size_t)((N*N+N)/2);
265
267
268 inline void assign(const MathSymmetricMatrix& v)
269 {
270 for(std::size_t i = 0; i < m_size; ++i)
271 {
272 m_data[i] = v[i];
273 }
274 }
275};
276
277
279template <std::size_t N>
280std::ostream& operator<< (std::ostream& outStream, const ug::MathSymmetricMatrix<N>& m)
281{
282 for(std::size_t i = 0; i < N; ++i)
283 {
284 for(std::size_t j = 0; j < N; ++j)
285 {
286 outStream << "[" << i << "][" << j << "]: " << std::scientific << std::setprecision(8) << std::setw(15) << m.entry(i, j) << std::endl;
287 }
288 }
289 return outStream;
290}
291
292std::ostream& operator<< (std::ostream& outStream, const ug::MathSymmetricMatrix<2>& m);
293std::ostream& operator<< (std::ostream& outStream, const ug::MathSymmetricMatrix<3>& m);
294
295
296} //end of namespace: lgmath
297
298
299#endif /* SYM_Matrix_H_ */
A class for fixed size, dense matrices.
Definition math_symmetric_matrix.h:63
static const std::size_t ColSize
Definition math_symmetric_matrix.h:68
MathSymmetricMatrix(const MathSymmetricMatrix &v)
Definition math_symmetric_matrix.h:72
MathSymmetricMatrix & operator/=(const value_type &val)
Divides all elements of the matrix by the given value.
Definition math_symmetric_matrix.h:167
value_type operator*(const MathSymmetricMatrix &v) const
Multiplies the matrix element-wise with another matrix and sums up the entries.
Definition math_symmetric_matrix.h:197
MathSymmetricMatrix & operator=(const MathSymmetricMatrix &v)
Assigns the elements of the given matrix to this one.
Definition math_symmetric_matrix.h:80
const value_type & operator[](std::size_t index) const
Definition math_symmetric_matrix.h:228
void scale_by_fnorm()
Definition math_symmetric_matrix.h:258
MathSymmetricMatrix()
Definition math_symmetric_matrix.h:71
std::size_t size_type
Definition math_symmetric_matrix.h:66
static const std::size_t RowSize
Definition math_symmetric_matrix.h:67
std::ostream & operator<<(std::ostream &outStream, const ug::MathMatrix< 2, 2 > &m)
Definition math_matrix.cpp:38
static const size_t m_size
Definition math_symmetric_matrix.h:264
value_type & operator[](std::size_t index)
Definition math_symmetric_matrix.h:227
T value_type
Definition math_symmetric_matrix.h:65
std::size_t num_rows() const
Definition math_symmetric_matrix.h:210
MathSymmetricMatrix & operator-=(const MathSymmetricMatrix &B)
Subtracts a matrix from 'this' one: .
Definition math_symmetric_matrix.h:107
MathSymmetricMatrix & operator*=(const value_type &val)
Multiplies all elements of the matrix with the given value.
Definition math_symmetric_matrix.h:182
value_type & operator()(std::size_t row, std::size_t col)
Definition math_symmetric_matrix.h:230
std::size_t num_cols() const
Definition math_symmetric_matrix.h:211
value_type m_data[m_size]
Definition math_symmetric_matrix.h:266
value_type fnorm()
Definition math_symmetric_matrix.h:245
void assign(const MathSymmetricMatrix &v)
Definition math_symmetric_matrix.h:268
value_type & entry(std::size_t row, std::size_t col)
Definition math_symmetric_matrix.h:213
MathSymmetricMatrix & operator+=(const MathSymmetricMatrix &B)
Adds a matrix to 'this' one: .
Definition math_symmetric_matrix.h:92
const value_type & entry(std::size_t row, std::size_t col) const
Definition math_symmetric_matrix.h:220
#define UG_ASSERT(expr, msg)
Definition assert.h:70
double number
Definition types.h:124
the ug namespace