OGS
EigenMatrix.h
Go to the documentation of this file.
1 
11 #pragma once
12 
13 #include <Eigen/Sparse>
14 #include <iosfwd>
15 #include <string>
16 
19 
20 namespace MathLib
21 {
27 class EigenMatrix final
28 {
29 public:
30  using RawMatrixType = Eigen::SparseMatrix<double, Eigen::RowMajor>;
31  using IndexType = RawMatrixType::Index;
32 
33  // TODO The matrix constructor should take num_rows and num_cols as
34  // arguments
35  // that is left for a later refactoring.
36 
40  explicit EigenMatrix(IndexType const n,
41  IndexType const n_nonzero_columns = 0)
42  : mat_(n, n)
43  {
44  if (n_nonzero_columns > 0)
45  {
46  mat_.reserve(Eigen::Matrix<IndexType, Eigen::Dynamic, 1>::Constant(
47  n, n_nonzero_columns));
48  }
49  }
50 
52  IndexType getNumberOfRows() const { return mat_.rows(); }
53 
55  IndexType getNumberOfColumns() const { return mat_.cols(); }
56 
58  static constexpr IndexType getRangeBegin() { return 0; }
59 
61  IndexType getRangeEnd() const { return getNumberOfRows(); }
62 
64  void setZero()
65  {
66  auto const N = mat_.nonZeros();
67  for (auto i = decltype(N){0}; i < N; i++)
68  {
69  mat_.valuePtr()[i] = 0;
70  }
71  // don't use mat_.setZero(). it makes a matrix uncompressed
72  }
73 
76  int setValue(IndexType row, IndexType col, double val)
77  {
78  assert(row < (IndexType)getNumberOfRows() &&
79  col < (IndexType)getNumberOfColumns());
80  mat_.coeffRef(row, col) = val;
81  return 0;
82  }
83 
86  int add(IndexType row, IndexType col, double val)
87  {
88  mat_.coeffRef(row, col) += val;
89  return 0;
90  }
91 
94  template <class T_DENSE_MATRIX>
95  void add(std::vector<IndexType> const& row_pos,
96  const T_DENSE_MATRIX& sub_matrix,
97  double fkt = 1.0)
98  {
99  this->add(row_pos, row_pos, sub_matrix, fkt);
100  }
101 
104  template <class T_DENSE_MATRIX>
105  void add(RowColumnIndices<IndexType> const& indices,
106  const T_DENSE_MATRIX& sub_matrix,
107  double fkt = 1.0)
108  {
109  this->add(indices.rows, indices.columns, sub_matrix, fkt);
110  }
111 
124  template <class T_DENSE_MATRIX>
125  void add(std::vector<IndexType> const& row_pos,
126  std::vector<IndexType> const& col_pos,
127  const T_DENSE_MATRIX& sub_matrix, double fkt = 1.0);
128 
130  double get(IndexType row, IndexType col) const
131  {
132  return mat_.coeff(row, col);
133  }
134 
136  static constexpr bool isAssembled() { return true; }
137 
139  void write(const std::string& filename) const;
140 
142  void write(std::ostream& os) const;
143 
145  const RawMatrixType& getRawMatrix() const { return mat_; }
146 
147 protected:
149 };
150 
151 template <class T_DENSE_MATRIX>
152 void EigenMatrix::add(std::vector<IndexType> const& row_pos,
153  std::vector<IndexType> const& col_pos,
154  const T_DENSE_MATRIX& sub_matrix, double fkt)
155 {
156  auto const n_rows = row_pos.size();
157  auto const n_cols = col_pos.size();
158  for (auto i = decltype(n_rows){0}; i < n_rows; i++)
159  {
160  auto const row = row_pos[i];
161  for (auto j = decltype(n_cols){0}; j < n_cols; j++)
162  {
163  auto const col = col_pos[j];
164  add(row, col, fkt * sub_matrix(i, j));
165  }
166  }
167 };
168 
170 template <typename SPARSITY_PATTERN>
171 struct SetMatrixSparsity<EigenMatrix, SPARSITY_PATTERN>
172 {
175  void operator()(EigenMatrix& matrix,
176  SPARSITY_PATTERN const& sparsity_pattern) const
177  {
178  static_assert(EigenMatrix::RawMatrixType::IsRowMajor,
179  "Set matrix sparsity relies on the EigenMatrix to be in "
180  "row-major storage order.");
181 
182  assert(matrix.getNumberOfRows() ==
183  static_cast<EigenMatrix::IndexType>(sparsity_pattern.size()));
184 
185  matrix.getRawMatrix().reserve(sparsity_pattern);
186  }
187 };
188 
189 } // end namespace MathLib
RawMatrixType::Index IndexType
Definition: EigenMatrix.h:31
IndexType getNumberOfColumns() const
return the number of columns
Definition: EigenMatrix.h:55
void add(std::vector< IndexType > const &row_pos, const T_DENSE_MATRIX &sub_matrix, double fkt=1.0)
Definition: EigenMatrix.h:95
static constexpr IndexType getRangeBegin()
return a start index of the active data range
Definition: EigenMatrix.h:58
void write(const std::string &filename) const
printout this matrix for debugging
Definition: EigenMatrix.cpp:19
void add(RowColumnIndices< IndexType > const &indices, const T_DENSE_MATRIX &sub_matrix, double fkt=1.0)
Definition: EigenMatrix.h:105
void setZero()
reset data entries to zero.
Definition: EigenMatrix.h:64
int add(IndexType row, IndexType col, double val)
Definition: EigenMatrix.h:86
RawMatrixType mat_
Definition: EigenMatrix.h:148
const RawMatrixType & getRawMatrix() const
Definition: EigenMatrix.h:145
IndexType getNumberOfRows() const
return the number of rows
Definition: EigenMatrix.h:52
int setValue(IndexType row, IndexType col, double val)
Definition: EigenMatrix.h:76
Eigen::SparseMatrix< double, Eigen::RowMajor > RawMatrixType
Definition: EigenMatrix.h:30
RawMatrixType & getRawMatrix()
Definition: EigenMatrix.h:144
double get(IndexType row, IndexType col) const
get value. This function returns zero if the element doesn't exist.
Definition: EigenMatrix.h:130
EigenMatrix(IndexType const n, IndexType const n_nonzero_columns=0)
Definition: EigenMatrix.h:40
IndexType getRangeEnd() const
return an end index of the active data range
Definition: EigenMatrix.h:61
static constexpr bool isAssembled()
return always true, i.e. the matrix is always ready for use
Definition: EigenMatrix.h:136
LineIndex const & columns
void operator()(EigenMatrix &matrix, SPARSITY_PATTERN const &sparsity_pattern) const
Definition: EigenMatrix.h:175