libMesh
eigen_sparse_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 
22 // Local includes
23 #include "libmesh/libmesh_config.h"
24 
25 #ifdef LIBMESH_HAVE_EIGEN
26 
27 #include "libmesh/eigen_sparse_vector.h"
28 #include "libmesh/eigen_sparse_matrix.h"
29 #include "libmesh/dense_matrix.h"
30 #include "libmesh/dof_map.h"
31 #include "libmesh/sparsity_pattern.h"
32 
33 namespace libMesh
34 {
35 
36 
37 //-----------------------------------------------------------------------
38 // EigenSparseMatrix members
39 template <typename T>
41  const numeric_index_type n_in,
44  const numeric_index_type nnz,
45  const numeric_index_type,
46  const numeric_index_type)
47 {
48  // noz ignored... only used for multiple processors!
49  libmesh_assert_equal_to (m_in, m_l);
50  libmesh_assert_equal_to (n_in, n_l);
51  libmesh_assert_equal_to (m_in, n_in);
52  libmesh_assert_greater (nnz, 0);
53 
54  _mat.resize(m_in, n_in);
55  _mat.reserve(Eigen::Matrix<numeric_index_type, Eigen::Dynamic, 1>::Constant(m_in,nnz));
56 
57  this->_is_initialized = true;
58 }
59 
60 
61 
62 template <typename T>
64 {
65  // Ignore calls on initialized objects
66  if (this->initialized())
67  return;
68 
69  // We need the DofMap for this!
70  libmesh_assert(this->_dof_map);
71 
72  // Clear initialized matrices
73  if (this->initialized())
74  this->clear();
75 
76  const numeric_index_type n_rows = this->_dof_map->n_dofs();
77  const numeric_index_type n_cols = n_rows;
78 
79 #ifndef NDEBUG
80  // The following variables are only used for assertions,
81  // so avoid declaring them when asserts are inactive.
82  const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(0);
83  const numeric_index_type m_l = n_l;
84 #endif
85 
86  // Laspack Matrices only work for uniprocessor cases
87  libmesh_assert_equal_to (n_rows, n_cols);
88  libmesh_assert_equal_to (m_l, n_rows);
89  libmesh_assert_equal_to (n_l, n_cols);
90 
91  const std::vector<numeric_index_type> & n_nz = this->_dof_map->get_n_nz();
92 
93 #ifndef NDEBUG
94  // The following variables are only used for assertions,
95  // so avoid declaring them when asserts are inactive.
96  const std::vector<numeric_index_type> & n_oz = this->_dof_map->get_n_oz();
97 #endif
98 
99  // Make sure the sparsity pattern isn't empty
100  libmesh_assert_equal_to (n_nz.size(), n_l);
101  libmesh_assert_equal_to (n_oz.size(), n_l);
102 
103  if (n_rows==0)
104  {
105  _mat.resize(0,0);
106  return;
107  }
108 
109  _mat.resize(n_rows,n_cols);
110  _mat.reserve(n_nz);
111 
112  this->_is_initialized = true;
113 
114  libmesh_assert_equal_to (n_rows, this->m());
115  libmesh_assert_equal_to (n_cols, this->n());
116 }
117 
118 
119 
120 template <typename T>
121 void EigenSparseMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
122  const std::vector<numeric_index_type> & rows,
123  const std::vector<numeric_index_type> & cols)
124 
125 {
126  libmesh_assert (this->initialized());
127  unsigned int n_rows = cast_int<unsigned int>(rows.size());
128  unsigned int n_cols = cast_int<unsigned int>(cols.size());
129  libmesh_assert_equal_to (dm.m(), n_rows);
130  libmesh_assert_equal_to (dm.n(), n_cols);
131 
132 
133  for (unsigned int i=0; i<n_rows; i++)
134  for (unsigned int j=0; j<n_cols; j++)
135  this->add(rows[i],cols[j],dm(i,j));
136 }
137 
138 
139 
140 template <typename T>
141 void EigenSparseMatrix<T>::get_diagonal (NumericVector<T> & dest_in) const
142 {
143  EigenSparseVector<T> & dest = cast_ref<EigenSparseVector<T> &>(dest_in);
144 
145  dest._vec = _mat.diagonal();
146 }
147 
148 
149 
150 template <typename T>
151 void EigenSparseMatrix<T>::get_transpose (SparseMatrix<T> & dest_in) const
152 {
153  EigenSparseMatrix<T> & dest = cast_ref<EigenSparseMatrix<T> &>(dest_in);
154 
155  dest._mat = _mat.transpose();
156 }
157 
158 
159 
160 template <typename T>
161 EigenSparseMatrix<T>::EigenSparseMatrix (const Parallel::Communicator & comm_in) :
162  SparseMatrix<T>(comm_in),
163  _closed (false)
164 {
165 }
166 
167 
168 
169 template <typename T>
170 EigenSparseMatrix<T>::~EigenSparseMatrix ()
171 {
172  this->clear ();
173 }
174 
175 
176 
177 template <typename T>
178 void EigenSparseMatrix<T>::clear ()
179 {
180  _mat.resize(0,0);
181 
182  _closed = false;
183  this->_is_initialized = false;
184 }
185 
186 
187 
188 template <typename T>
190 {
191  _mat.setZero();
192 }
193 
194 
195 
196 template <typename T>
197 numeric_index_type EigenSparseMatrix<T>::m () const
198 {
199  libmesh_assert (this->initialized());
200 
201  return _mat.rows();
202 }
203 
204 
205 
206 template <typename T>
207 numeric_index_type EigenSparseMatrix<T>::n () const
208 {
209  libmesh_assert (this->initialized());
210 
211  return _mat.cols();
212 }
213 
214 
215 
216 template <typename T>
217 numeric_index_type EigenSparseMatrix<T>::row_start () const
218 {
219  return 0;
220 }
221 
222 
223 
224 template <typename T>
225 numeric_index_type EigenSparseMatrix<T>::row_stop () const
226 {
227  return this->m();
228 }
229 
230 
231 
232 template <typename T>
233 void EigenSparseMatrix<T>::set (const numeric_index_type i,
234  const numeric_index_type j,
235  const T value)
236 {
237  libmesh_assert (this->initialized());
238  libmesh_assert_less (i, this->m());
239  libmesh_assert_less (j, this->n());
240 
241  _mat.coeffRef(i,j) = value;
242 }
243 
244 
245 
246 template <typename T>
247 void EigenSparseMatrix<T>::add (const numeric_index_type i,
248  const numeric_index_type j,
249  const T value)
250 {
251  libmesh_assert (this->initialized());
252  libmesh_assert_less (i, this->m());
253  libmesh_assert_less (j, this->n());
254 
255  _mat.coeffRef(i,j) += value;
256 }
257 
258 
259 
260 template <typename T>
261 void EigenSparseMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
262  const std::vector<numeric_index_type> & dof_indices)
263 {
264  this->add_matrix (dm, dof_indices, dof_indices);
265 }
266 
267 
268 
269 template <typename T>
270 void EigenSparseMatrix<T>::add (const T a_in, SparseMatrix<T> & X_in)
271 {
272  libmesh_assert (this->initialized());
273  libmesh_assert_equal_to (this->m(), X_in.m());
274  libmesh_assert_equal_to (this->n(), X_in.n());
275 
276  EigenSparseMatrix<T> & X = cast_ref<EigenSparseMatrix<T> &> (X_in);
277 
278  _mat += X._mat*a_in;
279 }
280 
281 
282 
283 
284 template <typename T>
285 T EigenSparseMatrix<T>::operator () (const numeric_index_type i,
286  const numeric_index_type j) const
287 {
288  libmesh_assert (this->initialized());
289  libmesh_assert_less (i, this->m());
290  libmesh_assert_less (j, this->n());
291 
292  return _mat.coeff(i,j);
293 }
294 
295 
296 
297 template <typename T>
298 Real EigenSparseMatrix<T>::l1_norm () const
299 {
300  // There does not seem to be a straightforward way to iterate over
301  // the columns of an EigenSparseMatrix. So we use some extra
302  // storage and keep track of the column sums while going over the
303  // row entries...
304  std::vector<Real> abs_col_sums(this->n());
305 
306  // For a row-major Eigen SparseMatrix like we're using, the
307  // InnerIterator iterates over the non-zero entries of rows.
308  for (unsigned row=0; row<this->m(); ++row)
309  {
310  EigenSM::InnerIterator it(_mat, row);
311  for (; it; ++it)
312  abs_col_sums[it.col()] += std::abs(it.value());
313  }
314 
315  return *(std::max_element(abs_col_sums.begin(), abs_col_sums.end()));
316 }
317 
318 
319 
320 template <typename T>
321 Real EigenSparseMatrix<T>::linfty_norm () const
322 {
323  Real max_abs_row_sum = 0.;
324 
325  // For a row-major Eigen SparseMatrix like we're using, the
326  // InnerIterator iterates over the non-zero entries of rows.
327  for (unsigned row=0; row<this->m(); ++row)
328  {
329  Real current_abs_row_sum = 0.;
330  EigenSM::InnerIterator it(_mat, row);
331  for (; it; ++it)
332  current_abs_row_sum += std::abs(it.value());
333 
334  max_abs_row_sum = std::max(max_abs_row_sum, current_abs_row_sum);
335  }
336 
337  return max_abs_row_sum;
338 }
339 
340 
341 //------------------------------------------------------------------
342 // Explicit instantiations
343 template class EigenSparseMatrix<Number>;
344 
345 } // namespace libMesh
346 
347 
348 #endif // #ifdef LIBMESH_HAVE_EIGEN
double abs(double a)
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:178
long double max(long double a, double b)
libmesh_assert(j)
dof_id_type numeric_index_type
Definition: id_types.h:92
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:255
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
PetscErrorCode Vec Mat libmesh_dbg_var(j)
bool _is_initialized
Flag indicating if the data structures have been initialized.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const bool value
Definition: xdr_io.C:108
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:274
bool initialized() const
virtual void clear() libmesh_override
Release all memory and clear data structures.