ug4
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 
43 namespace ug
44 {
45 
52 //template <std::size_t N, std::size_t M, typename T = number> class MathSymmetricMatrix;
53 
61 template <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} }
245  inline value_type fnorm(){
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 
279 template <std::size_t N>
280 std::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 
292 std::ostream& operator<< (std::ostream& outStream, const ug::MathSymmetricMatrix<2>& m);
293 std::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
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
MathSymmetricMatrix(const MathSymmetricMatrix &v)
Definition: math_symmetric_matrix.h:72
MathSymmetricMatrix & operator-=(const MathSymmetricMatrix &B)
Subtracts a matrix from 'this' one: .
Definition: math_symmetric_matrix.h:107
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
const value_type & entry(std::size_t row, std::size_t col) const
Definition: math_symmetric_matrix.h:220
const value_type & operator[](std::size_t index) const
Definition: math_symmetric_matrix.h:228
value_type & operator[](std::size_t index)
Definition: math_symmetric_matrix.h:227
MathSymmetricMatrix & operator/=(const value_type &val)
Divides all elements of the matrix by the given value.
Definition: math_symmetric_matrix.h:167
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
MathSymmetricMatrix & operator=(const MathSymmetricMatrix &v)
Assigns the elements of the given matrix to this one.
Definition: math_symmetric_matrix.h:80
static const size_t m_size
Definition: math_symmetric_matrix.h:264
T value_type
Definition: math_symmetric_matrix.h:65
value_type & operator()(std::size_t row, std::size_t col)
Definition: math_symmetric_matrix.h:230
std::size_t num_rows() const
Definition: math_symmetric_matrix.h:210
std::size_t num_cols() const
Definition: math_symmetric_matrix.h:211
MathSymmetricMatrix & operator*=(const value_type &val)
Multiplies all elements of the matrix with the given value.
Definition: math_symmetric_matrix.h:182
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
#define UG_ASSERT(expr, msg)
Definition: assert.h:70
double number
Definition: types.h:124
the ug namespace