Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 : // Local includes
21 : #include "libmesh/libmesh_config.h"
22 :
23 : #ifdef LIBMESH_HAVE_EIGEN
24 :
25 : #include "libmesh/eigen_sparse_vector.h"
26 : #include "libmesh/eigen_sparse_matrix.h"
27 : #include "libmesh/dense_matrix.h"
28 : #include "libmesh/dof_map.h"
29 : #include "libmesh/sparsity_pattern.h"
30 :
31 : // C++ Includes
32 : #include <memory>
33 :
34 :
35 : namespace libMesh
36 : {
37 :
38 :
39 : //-----------------------------------------------------------------------
40 : // EigenSparseMatrix members
41 : template <typename T>
42 927 : void EigenSparseMatrix<T>::init (const numeric_index_type m_in,
43 : const numeric_index_type n_in,
44 : const numeric_index_type libmesh_dbg_var(m_l),
45 : const numeric_index_type libmesh_dbg_var(n_l),
46 : const numeric_index_type nnz,
47 : const numeric_index_type,
48 : const numeric_index_type)
49 : {
50 : // noz ignored... only used for multiple processors!
51 26 : libmesh_assert_equal_to (m_in, m_l);
52 26 : libmesh_assert_equal_to (n_in, n_l);
53 26 : libmesh_assert_greater (nnz, 0);
54 :
55 927 : _mat.resize(m_in, n_in);
56 927 : _mat.reserve(Eigen::Matrix<numeric_index_type, Eigen::Dynamic, 1>::Constant(m_in,nnz));
57 :
58 927 : this->_is_initialized = true;
59 927 : }
60 :
61 :
62 :
63 : template <typename T>
64 550 : void EigenSparseMatrix<T>::init (const ParallelType)
65 : {
66 : // Ignore calls on initialized objects
67 550 : if (this->initialized())
68 0 : return;
69 :
70 : // We need the DofMap for this!
71 0 : libmesh_assert(this->_dof_map);
72 :
73 : // Clear initialized matrices
74 0 : if (this->initialized())
75 0 : this->clear();
76 :
77 550 : const numeric_index_type n_rows = this->_dof_map->n_dofs();
78 0 : const numeric_index_type n_cols = n_rows;
79 :
80 : #ifndef NDEBUG
81 : // The following variables are only used for assertions,
82 : // so avoid declaring them when asserts are inactive.
83 0 : const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(0);
84 0 : const numeric_index_type m_l = n_l;
85 : #endif
86 :
87 : // Eigen Matrices only work for uniprocessor cases
88 0 : libmesh_assert_equal_to (m_l, n_rows);
89 0 : libmesh_assert_equal_to (n_l, n_cols);
90 :
91 550 : const std::vector<numeric_index_type> & n_nz = this->_sp->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 0 : const std::vector<numeric_index_type> & n_oz = this->_sp->get_n_oz();
97 : #endif
98 :
99 : // Make sure the sparsity pattern isn't empty
100 0 : libmesh_assert_equal_to (n_nz.size(), n_l);
101 0 : libmesh_assert_equal_to (n_oz.size(), n_l);
102 :
103 550 : if (n_rows==0)
104 : {
105 0 : _mat.resize(0,0);
106 0 : return;
107 : }
108 :
109 550 : _mat.resize(n_rows,n_cols);
110 0 : _mat.reserve(n_nz);
111 :
112 550 : this->_is_initialized = true;
113 :
114 0 : libmesh_assert_equal_to (n_rows, this->m());
115 0 : libmesh_assert_equal_to (n_cols, this->n());
116 : }
117 :
118 :
119 :
120 : template <typename T>
121 1044562 : 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 6 : libmesh_assert (this->initialized());
127 12 : unsigned int n_rows = cast_int<unsigned int>(rows.size());
128 12 : unsigned int n_cols = cast_int<unsigned int>(cols.size());
129 6 : libmesh_assert_equal_to (dm.m(), n_rows);
130 6 : libmesh_assert_equal_to (dm.n(), n_cols);
131 :
132 :
133 7820206 : for (unsigned int i=0; i<n_rows; i++)
134 133404971 : for (unsigned int j=0; j<n_cols; j++)
135 126629519 : this->add(rows[i],cols[j],dm(i,j));
136 1044562 : }
137 :
138 :
139 :
140 : template <typename T>
141 0 : void EigenSparseMatrix<T>::get_diagonal (NumericVector<T> & dest_in) const
142 : {
143 0 : EigenSparseVector<T> & dest = cast_ref<EigenSparseVector<T> &>(dest_in);
144 :
145 0 : dest._vec = _mat.diagonal();
146 0 : }
147 :
148 :
149 :
150 : template <typename T>
151 76 : void EigenSparseMatrix<T>::get_transpose (SparseMatrix<T> & dest_in) const
152 : {
153 0 : EigenSparseMatrix<T> & dest = cast_ref<EigenSparseMatrix<T> &>(dest_in);
154 :
155 76 : dest._mat = _mat.transpose();
156 76 : }
157 :
158 :
159 :
160 : template <typename T>
161 962 : EigenSparseMatrix<T>::EigenSparseMatrix (const Parallel::Communicator & comm_in) :
162 : SparseMatrix<T>(comm_in),
163 962 : _closed (false)
164 : {
165 962 : }
166 :
167 :
168 :
169 : template <typename T>
170 806 : void EigenSparseMatrix<T>::clear ()
171 : {
172 806 : _mat.resize(0,0);
173 :
174 806 : _closed = false;
175 806 : this->_is_initialized = false;
176 806 : }
177 :
178 :
179 :
180 : template <typename T>
181 1946 : void EigenSparseMatrix<T>::zero ()
182 : {
183 : // This doesn't just zero, it clears the entire non-zero structure!
184 1946 : _mat.setZero();
185 :
186 1946 : if (this->_sp)
187 : {
188 : // Re-reserve our non-zero structure
189 0 : const std::vector<numeric_index_type> & n_nz = this->_sp->get_n_nz();
190 0 : _mat.reserve(n_nz);
191 : }
192 1946 : }
193 :
194 :
195 :
196 : template <typename T>
197 71 : std::unique_ptr<SparseMatrix<T>> EigenSparseMatrix<T>::zero_clone () const
198 : {
199 : // TODO: If there is a more efficient way to make a zeroed-out copy
200 : // of an EigenSM, we should call that instead.
201 73 : auto ret = std::make_unique<EigenSparseMatrix<T>>(*this);
202 71 : ret->zero();
203 :
204 73 : return ret;
205 : }
206 :
207 :
208 :
209 : template <typename T>
210 213 : std::unique_ptr<SparseMatrix<T>> EigenSparseMatrix<T>::clone () const
211 : {
212 213 : return std::make_unique<EigenSparseMatrix<T>>(*this);
213 : }
214 :
215 :
216 :
217 : template <typename T>
218 4179 : numeric_index_type EigenSparseMatrix<T>::m () const
219 : {
220 3818 : libmesh_assert (this->initialized());
221 :
222 4257 : return cast_int<numeric_index_type>(_mat.rows());
223 : }
224 :
225 :
226 :
227 : template <typename T>
228 4218 : numeric_index_type EigenSparseMatrix<T>::n () const
229 : {
230 3795 : libmesh_assert (this->initialized());
231 :
232 4273 : return cast_int<numeric_index_type>(_mat.cols());
233 : }
234 :
235 :
236 :
237 : template <typename T>
238 5700 : numeric_index_type EigenSparseMatrix<T>::row_start () const
239 : {
240 5700 : return 0;
241 : }
242 :
243 :
244 :
245 : template <typename T>
246 943 : numeric_index_type EigenSparseMatrix<T>::row_stop () const
247 : {
248 943 : return this->m();
249 : }
250 :
251 :
252 :
253 : template <typename T>
254 779 : numeric_index_type EigenSparseMatrix<T>::col_start () const
255 : {
256 779 : return 0;
257 : }
258 :
259 :
260 :
261 : template <typename T>
262 284 : numeric_index_type EigenSparseMatrix<T>::col_stop () const
263 : {
264 284 : return this->n();
265 : }
266 :
267 :
268 :
269 : template <typename T>
270 126867 : void EigenSparseMatrix<T>::set (const numeric_index_type i,
271 : const numeric_index_type j,
272 : const T value)
273 : {
274 3572 : libmesh_assert (this->initialized());
275 3572 : libmesh_assert_less (i, this->m());
276 3572 : libmesh_assert_less (j, this->n());
277 :
278 126867 : _mat.coeffRef(i,j) = value;
279 126867 : }
280 :
281 :
282 :
283 : template <typename T>
284 17017 : void EigenSparseMatrix<T>::add (const numeric_index_type i,
285 : const numeric_index_type j,
286 : const T value)
287 : {
288 96 : libmesh_assert (this->initialized());
289 96 : libmesh_assert_less (i, this->m());
290 96 : libmesh_assert_less (j, this->n());
291 :
292 126646152 : _mat.coeffRef(i,j) += value;
293 17017 : }
294 :
295 :
296 :
297 : template <typename T>
298 1003727 : void EigenSparseMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
299 : const std::vector<numeric_index_type> & dof_indices)
300 : {
301 1003727 : this->add_matrix (dm, dof_indices, dof_indices);
302 1003727 : }
303 :
304 :
305 :
306 : template <typename T>
307 305 : void EigenSparseMatrix<T>::add (const T a_in, const SparseMatrix<T> & X_in)
308 : {
309 4 : libmesh_assert (this->initialized());
310 4 : libmesh_assert_equal_to (this->m(), X_in.m());
311 4 : libmesh_assert_equal_to (this->n(), X_in.n());
312 :
313 : const EigenSparseMatrix<T> & X =
314 4 : cast_ref<const EigenSparseMatrix<T> &> (X_in);
315 :
316 309 : _mat += X._mat*a_in;
317 305 : }
318 :
319 :
320 :
321 :
322 : template <typename T>
323 1520 : T EigenSparseMatrix<T>::operator () (const numeric_index_type i,
324 : const numeric_index_type j) const
325 : {
326 64 : libmesh_assert (this->initialized());
327 64 : libmesh_assert_less (i, this->m());
328 64 : libmesh_assert_less (j, this->n());
329 :
330 1520 : return _mat.coeff(i,j);
331 : }
332 :
333 :
334 :
335 : template <typename T>
336 1065 : Real EigenSparseMatrix<T>::l1_norm () const
337 : {
338 : // There does not seem to be a straightforward way to iterate over
339 : // the columns of an EigenSparseMatrix. So we use some extra
340 : // storage and keep track of the column sums while going over the
341 : // row entries...
342 1065 : std::vector<Real> abs_col_sums(this->n());
343 :
344 : // For a row-major Eigen SparseMatrix like we're using, the
345 : // InnerIterator iterates over the non-zero entries of rows.
346 19413 : for (auto row : make_range(this->m()))
347 : {
348 17802 : EigenSM::InnerIterator it(_mat, row);
349 151940 : for (; it; ++it)
350 141150 : abs_col_sums[it.col()] += std::abs(it.value());
351 : }
352 :
353 1125 : return *(std::max_element(abs_col_sums.begin(), abs_col_sums.end()));
354 : }
355 :
356 :
357 :
358 : template <typename T>
359 426 : Real EigenSparseMatrix<T>::linfty_norm () const
360 : {
361 426 : Real max_abs_row_sum = 0.;
362 :
363 : // For a row-major Eigen SparseMatrix like we're using, the
364 : // InnerIterator iterates over the non-zero entries of rows.
365 16188 : for (auto row : make_range(this->m()))
366 : {
367 15762 : Real current_abs_row_sum = 0.;
368 15318 : EigenSM::InnerIterator it(_mat, row);
369 140296 : for (; it; ++it)
370 128042 : current_abs_row_sum += std::abs(it.value());
371 :
372 15762 : max_abs_row_sum = std::max(max_abs_row_sum, current_abs_row_sum);
373 : }
374 :
375 426 : return max_abs_row_sum;
376 : }
377 :
378 :
379 :
380 : template <typename T>
381 1208 : void EigenSparseMatrix<T>::get_row(numeric_index_type i,
382 : std::vector<numeric_index_type> & indices,
383 : std::vector<T> & values) const
384 : {
385 32 : indices.clear();
386 32 : values.clear();
387 :
388 : // InnerIterator is over rows in RowMajor ordering
389 : static_assert(EigenSM::IsRowMajor);
390 :
391 5826 : for (EigenSM::InnerIterator it(_mat, i); it; ++it)
392 : {
393 4618 : indices.push_back(it.col());
394 4618 : values.push_back(it.value());
395 : }
396 1208 : }
397 :
398 :
399 :
400 : //------------------------------------------------------------------
401 : // Explicit instantiations
402 : template class LIBMESH_EXPORT EigenSparseMatrix<Number>;
403 :
404 : } // namespace libMesh
405 :
406 :
407 : #endif // #ifdef LIBMESH_HAVE_EIGEN
|
