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 1073 : 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 30 : libmesh_assert_equal_to (m_in, m_l);
52 30 : libmesh_assert_equal_to (n_in, n_l);
53 30 : libmesh_assert_greater (nnz, 0);
54 :
55 1073 : _mat.resize(m_in, n_in);
56 1073 : _mat.reserve(Eigen::Matrix<numeric_index_type, Eigen::Dynamic, 1>::Constant(m_in,nnz));
57 :
58 1073 : this->_is_initialized = true;
59 1073 : }
60 :
61 :
62 :
63 : template <typename T>
64 563 : void EigenSparseMatrix<T>::init (const ParallelType)
65 : {
66 : // Ignore calls on initialized objects
67 563 : 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 563 : 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 563 : 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 563 : if (n_rows==0)
104 : {
105 0 : _mat.resize(0,0);
106 0 : return;
107 : }
108 :
109 563 : _mat.resize(n_rows,n_cols);
110 0 : _mat.reserve(n_nz);
111 :
112 563 : 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 1058851 : 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 8 : libmesh_assert (this->initialized());
127 16 : unsigned int n_rows = cast_int<unsigned int>(rows.size());
128 16 : unsigned int n_cols = cast_int<unsigned int>(cols.size());
129 8 : libmesh_assert_equal_to (dm.m(), n_rows);
130 8 : libmesh_assert_equal_to (dm.n(), n_cols);
131 :
132 :
133 8073654 : for (unsigned int i=0; i<n_rows; i++)
134 133793805 : for (unsigned int j=0; j<n_cols; j++)
135 126779258 : this->add(rows[i],cols[j],dm(i,j));
136 1058851 : }
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 :
147 0 : dest.close();
148 0 : }
149 :
150 :
151 :
152 : template <typename T>
153 147 : void EigenSparseMatrix<T>::get_transpose (SparseMatrix<T> & dest_in) const
154 : {
155 2 : EigenSparseMatrix<T> & dest = cast_ref<EigenSparseMatrix<T> &>(dest_in);
156 :
157 147 : dest._mat = _mat.transpose();
158 :
159 147 : dest._is_initialized = true;
160 147 : dest._closed = true;
161 147 : }
162 :
163 :
164 :
165 : template <typename T>
166 1190 : EigenSparseMatrix<T>::EigenSparseMatrix (const Parallel::Communicator & comm_in) :
167 : SparseMatrix<T>(comm_in),
168 1190 : _closed (false)
169 : {
170 1190 : }
171 :
172 :
173 :
174 : template <typename T>
175 821 : void EigenSparseMatrix<T>::clear ()
176 : {
177 821 : _mat.resize(0,0);
178 :
179 821 : _closed = false;
180 821 : this->_is_initialized = false;
181 821 : }
182 :
183 :
184 :
185 : template <typename T>
186 2496 : void EigenSparseMatrix<T>::zero ()
187 : {
188 : // This doesn't just zero, it clears the entire non-zero structure!
189 2496 : _mat.setZero();
190 :
191 2496 : if (this->_sp)
192 : {
193 : // Re-reserve our non-zero structure
194 0 : const std::vector<numeric_index_type> & n_nz = this->_sp->get_n_nz();
195 0 : _mat.reserve(n_nz);
196 : }
197 2496 : }
198 :
199 :
200 :
201 : template <typename T>
202 71 : std::unique_ptr<SparseMatrix<T>> EigenSparseMatrix<T>::zero_clone () const
203 : {
204 : // TODO: If there is a more efficient way to make a zeroed-out copy
205 : // of an EigenSM, we should call that instead.
206 73 : auto ret = std::make_unique<EigenSparseMatrix<T>>(*this);
207 71 : ret->zero();
208 :
209 73 : return ret;
210 : }
211 :
212 :
213 :
214 : template <typename T>
215 213 : std::unique_ptr<SparseMatrix<T>> EigenSparseMatrix<T>::clone () const
216 : {
217 213 : return std::make_unique<EigenSparseMatrix<T>>(*this);
218 : }
219 :
220 :
221 :
222 : template <typename T>
223 4227 : numeric_index_type EigenSparseMatrix<T>::m () const
224 : {
225 3858 : libmesh_assert (this->initialized());
226 :
227 4313 : return cast_int<numeric_index_type>(_mat.rows());
228 : }
229 :
230 :
231 :
232 : template <typename T>
233 4258 : numeric_index_type EigenSparseMatrix<T>::n () const
234 : {
235 3831 : libmesh_assert (this->initialized());
236 :
237 4317 : return cast_int<numeric_index_type>(_mat.cols());
238 : }
239 :
240 :
241 :
242 : template <typename T>
243 5771 : numeric_index_type EigenSparseMatrix<T>::row_start () const
244 : {
245 5771 : return 0;
246 : }
247 :
248 :
249 :
250 : template <typename T>
251 943 : numeric_index_type EigenSparseMatrix<T>::row_stop () const
252 : {
253 943 : return this->m();
254 : }
255 :
256 :
257 :
258 : template <typename T>
259 850 : numeric_index_type EigenSparseMatrix<T>::col_start () const
260 : {
261 850 : return 0;
262 : }
263 :
264 :
265 :
266 : template <typename T>
267 284 : numeric_index_type EigenSparseMatrix<T>::col_stop () const
268 : {
269 284 : return this->n();
270 : }
271 :
272 :
273 :
274 : template <typename T>
275 127371 : void EigenSparseMatrix<T>::set (const numeric_index_type i,
276 : const numeric_index_type j,
277 : const T value)
278 : {
279 3572 : libmesh_assert (this->initialized());
280 3572 : libmesh_assert_less (i, this->m());
281 3572 : libmesh_assert_less (j, this->n());
282 :
283 127371 : _mat.coeffRef(i,j) = value;
284 127371 : }
285 :
286 :
287 :
288 : template <typename T>
289 17081 : void EigenSparseMatrix<T>::add (const numeric_index_type i,
290 : const numeric_index_type j,
291 : const T value)
292 : {
293 128 : libmesh_assert (this->initialized());
294 128 : libmesh_assert_less (i, this->m());
295 128 : libmesh_assert_less (j, this->n());
296 :
297 126795827 : _mat.coeffRef(i,j) += value;
298 17081 : }
299 :
300 :
301 :
302 : template <typename T>
303 1017945 : void EigenSparseMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
304 : const std::vector<numeric_index_type> & dof_indices)
305 : {
306 1017945 : this->add_matrix (dm, dof_indices, dof_indices);
307 1017945 : }
308 :
309 :
310 :
311 : template <typename T>
312 305 : void EigenSparseMatrix<T>::add (const T a_in, const SparseMatrix<T> & X_in)
313 : {
314 4 : libmesh_assert (this->initialized());
315 4 : libmesh_assert_equal_to (this->m(), X_in.m());
316 4 : libmesh_assert_equal_to (this->n(), X_in.n());
317 :
318 : const EigenSparseMatrix<T> & X =
319 4 : cast_ref<const EigenSparseMatrix<T> &> (X_in);
320 :
321 309 : _mat += X._mat*a_in;
322 305 : }
323 :
324 :
325 :
326 :
327 : template <typename T>
328 1520 : T EigenSparseMatrix<T>::operator () (const numeric_index_type i,
329 : const numeric_index_type j) const
330 : {
331 64 : libmesh_assert (this->initialized());
332 64 : libmesh_assert_less (i, this->m());
333 64 : libmesh_assert_less (j, this->n());
334 :
335 1520 : return _mat.coeff(i,j);
336 : }
337 :
338 :
339 :
340 : template <typename T>
341 1211 : Real EigenSparseMatrix<T>::l1_norm () const
342 : {
343 : // There does not seem to be a straightforward way to iterate over
344 : // the columns of an EigenSparseMatrix. So we use some extra
345 : // storage and keep track of the column sums while going over the
346 : // row entries...
347 1211 : std::vector<Real> abs_col_sums(this->n());
348 :
349 : // For a row-major Eigen SparseMatrix like we're using, the
350 : // InnerIterator iterates over the non-zero entries of rows.
351 20239 : for (auto row : make_range(this->m()))
352 : {
353 18462 : EigenSM::InnerIterator it(_mat, row);
354 155392 : for (; it; ++it)
355 144054 : abs_col_sums[it.col()] += std::abs(it.value());
356 : }
357 :
358 1279 : return *(std::max_element(abs_col_sums.begin(), abs_col_sums.end()));
359 : }
360 :
361 :
362 :
363 : template <typename T>
364 572 : Real EigenSparseMatrix<T>::linfty_norm () const
365 : {
366 572 : Real max_abs_row_sum = 0.;
367 :
368 : // For a row-major Eigen SparseMatrix like we're using, the
369 : // InnerIterator iterates over the non-zero entries of rows.
370 17010 : for (auto row : make_range(this->m()))
371 : {
372 16438 : Real current_abs_row_sum = 0.;
373 15978 : EigenSM::InnerIterator it(_mat, row);
374 143748 : for (; it; ++it)
375 130882 : current_abs_row_sum += std::abs(it.value());
376 :
377 16438 : max_abs_row_sum = std::max(max_abs_row_sum, current_abs_row_sum);
378 : }
379 :
380 572 : return max_abs_row_sum;
381 : }
382 :
383 :
384 :
385 : template <typename T>
386 1208 : void EigenSparseMatrix<T>::get_row(numeric_index_type i,
387 : std::vector<numeric_index_type> & indices,
388 : std::vector<T> & values) const
389 : {
390 32 : indices.clear();
391 32 : values.clear();
392 :
393 : // InnerIterator is over rows in RowMajor ordering
394 : static_assert(EigenSM::IsRowMajor);
395 :
396 5826 : for (EigenSM::InnerIterator it(_mat, i); it; ++it)
397 : {
398 4618 : indices.push_back(it.col());
399 4618 : values.push_back(it.value());
400 : }
401 1208 : }
402 :
403 :
404 :
405 : //------------------------------------------------------------------
406 : // Explicit instantiations
407 : template class LIBMESH_EXPORT EigenSparseMatrix<Number>;
408 :
409 : } // namespace libMesh
410 :
411 :
412 : #endif // #ifdef LIBMESH_HAVE_EIGEN
|
