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 : #include "libmesh/diagonal_matrix.h"
19 : #include "libmesh/numeric_vector.h"
20 : #include "libmesh/dense_matrix.h"
21 : #include "libmesh/dof_map.h"
22 : #include "libmesh/libmesh_common.h"
23 :
24 :
25 : // C++ Includes
26 : #include <memory>
27 :
28 :
29 : namespace libMesh
30 : {
31 :
32 : template <typename T>
33 639 : DiagonalMatrix<T>::DiagonalMatrix(const Parallel::Communicator & comm_in) : SparseMatrix<T>(comm_in)
34 : {
35 1242 : _diagonal = NumericVector<T>::build(comm_in);
36 639 : }
37 :
38 : template <typename T>
39 : DiagonalMatrix<T> &
40 4 : DiagonalMatrix<T>::operator=(const NumericVector<T> & vec)
41 : {
42 71 : *_diagonal = vec;
43 4 : return *this;
44 : }
45 :
46 : template <typename T>
47 : DiagonalMatrix<T> &
48 75 : DiagonalMatrix<T>::operator=(NumericVector<T> && vec)
49 : {
50 : // Don't get confused by the &&: vec is an lvalue reference; the && just
51 : // indicates that we are receiving an object that is safe to move from. Note
52 : // that we are not going to use std::move here because we do not have
53 : // (virtual) move assignment operations defined for NumericVector sub-classes
54 142 : _diagonal->swap(vec);
55 75 : return *this;
56 : }
57 :
58 : template <typename T>
59 : void
60 284 : DiagonalMatrix<T>::init(const numeric_index_type m,
61 : const numeric_index_type /*n*/,
62 : const numeric_index_type m_l,
63 : const numeric_index_type /*n_l*/,
64 : const numeric_index_type /*nnz*/,
65 : const numeric_index_type /*noz*/,
66 : const numeric_index_type /*blocksize*/)
67 : {
68 284 : _diagonal->init(m, m_l);
69 284 : }
70 :
71 : template <typename T>
72 : void
73 0 : DiagonalMatrix<T>::init(const ParallelType type)
74 : {
75 0 : libmesh_assert(this->_dof_map);
76 :
77 0 : _diagonal->init(this->_dof_map->n_dofs(),
78 0 : this->_dof_map->n_local_dofs(),
79 : /*fast=*/false,
80 : type);
81 0 : }
82 :
83 : template <typename T>
84 : void
85 213 : DiagonalMatrix<T>::init(const NumericVector<T> & other, const bool fast)
86 : {
87 213 : _diagonal->init(other, fast);
88 213 : }
89 :
90 : template <typename T>
91 : void
92 213 : DiagonalMatrix<T>::init(const DiagonalMatrix<T> & other, const bool fast)
93 : {
94 213 : init(other.diagonal(), fast);
95 213 : }
96 :
97 : template <typename T>
98 : void
99 0 : DiagonalMatrix<T>::clear()
100 : {
101 0 : _diagonal->clear();
102 0 : }
103 :
104 : template <typename T>
105 : void
106 213 : DiagonalMatrix<T>::zero()
107 : {
108 213 : _diagonal->zero();
109 213 : }
110 :
111 : template <typename T>
112 71 : std::unique_ptr<SparseMatrix<T>> DiagonalMatrix<T>::zero_clone () const
113 : {
114 : // Make empty copy with matching comm
115 73 : auto mat_copy = std::make_unique<DiagonalMatrix<T>>(this->comm());
116 :
117 : // Initialize copy with our same nonzero structure, and explicitly
118 : // zero values using fast == false.
119 71 : mat_copy->init(*this, /*fast=*/false);
120 :
121 73 : return mat_copy;
122 67 : }
123 :
124 :
125 :
126 : template <typename T>
127 71 : std::unique_ptr<SparseMatrix<T>> DiagonalMatrix<T>::clone () const
128 : {
129 : // Make empty copy with matching comm
130 73 : auto mat_copy = std::make_unique<DiagonalMatrix<T>>(this->comm());
131 :
132 : // Make copy of our diagonal
133 73 : auto diag_copy = _diagonal->clone();
134 :
135 : // Swap diag_copy with diagonal in mat_copy
136 4 : *mat_copy = std::move(*diag_copy);
137 :
138 73 : return mat_copy;
139 67 : }
140 :
141 : template <typename T>
142 : void
143 426 : DiagonalMatrix<T>::close()
144 : {
145 426 : _diagonal->close();
146 426 : }
147 :
148 : template <typename T>
149 : numeric_index_type
150 657 : DiagonalMatrix<T>::m() const
151 : {
152 639 : return _diagonal->size();
153 : }
154 :
155 : template <typename T>
156 : numeric_index_type
157 639 : DiagonalMatrix<T>::n() const
158 : {
159 639 : return _diagonal->size();
160 : }
161 :
162 : template <typename T>
163 : numeric_index_type
164 1491 : DiagonalMatrix<T>::row_start() const
165 : {
166 1491 : return _diagonal->first_local_index();
167 : }
168 :
169 : template <typename T>
170 : numeric_index_type
171 1491 : DiagonalMatrix<T>::row_stop() const
172 : {
173 1491 : return _diagonal->last_local_index();
174 : }
175 :
176 : template <typename T>
177 : numeric_index_type
178 0 : DiagonalMatrix<T>::col_start() const
179 : {
180 0 : return _diagonal->first_local_index();
181 : }
182 :
183 : template <typename T>
184 : numeric_index_type
185 0 : DiagonalMatrix<T>::col_stop() const
186 : {
187 0 : return _diagonal->last_local_index();
188 : }
189 :
190 : template <typename T>
191 : void
192 524 : DiagonalMatrix<T>::set(const numeric_index_type i, const numeric_index_type j, const T value)
193 : {
194 524 : if (i == j)
195 453 : _diagonal->set(i, value);
196 524 : }
197 :
198 : template <typename T>
199 : void
200 595 : DiagonalMatrix<T>::add(const numeric_index_type i, const numeric_index_type j, const T value)
201 : {
202 595 : if (i == j)
203 524 : _diagonal->add(i, value);
204 595 : }
205 :
206 : template <typename T>
207 : void
208 71 : DiagonalMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
209 : const std::vector<numeric_index_type> & rows,
210 : const std::vector<numeric_index_type> & cols)
211 : {
212 4 : auto m = dm.m();
213 4 : auto n = dm.n();
214 :
215 524 : for (decltype(m) i = 0; i < m; ++i)
216 4364 : for (decltype(n) j = 0; j < n; ++j)
217 : {
218 3911 : auto global_i = rows[i];
219 3911 : auto global_j = cols[j];
220 3911 : if (global_i == global_j)
221 458 : _diagonal->add(global_i, dm(i, j));
222 : }
223 71 : }
224 :
225 : template <typename T>
226 : void
227 71 : DiagonalMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
228 : const std::vector<numeric_index_type> & dof_indices)
229 : {
230 71 : _diagonal->add_vector(dm.diagonal(), dof_indices);
231 71 : }
232 :
233 : template <typename T>
234 : void
235 71 : DiagonalMatrix<T>::add(const T a, const SparseMatrix<T> & X)
236 : {
237 73 : auto x_diagonal = _diagonal->zero_clone();
238 73 : X.get_diagonal(*x_diagonal);
239 73 : _diagonal->add(a, *x_diagonal);
240 71 : }
241 :
242 : template <typename T>
243 : T
244 22904 : DiagonalMatrix<T>::operator()(const numeric_index_type i, const numeric_index_type j) const
245 : {
246 22904 : if (i == j)
247 5460 : return (*_diagonal)(i);
248 : else
249 68 : return 0;
250 : }
251 :
252 : template <typename T>
253 : Real
254 284 : DiagonalMatrix<T>::l1_norm() const
255 : {
256 284 : return _diagonal->l1_norm();
257 : }
258 :
259 : template <typename T>
260 : Real
261 71 : DiagonalMatrix<T>::linfty_norm() const
262 : {
263 71 : return _diagonal->linfty_norm();
264 : }
265 :
266 : template <typename T>
267 : bool
268 284 : DiagonalMatrix<T>::closed() const
269 : {
270 284 : return _diagonal->closed();
271 : }
272 :
273 : template <typename T>
274 : void
275 0 : DiagonalMatrix<T>::print_personal(std::ostream & os) const
276 : {
277 0 : _diagonal->print(os);
278 0 : }
279 :
280 :
281 :
282 : template <typename T>
283 : void
284 142 : DiagonalMatrix<T>::get_diagonal(NumericVector<T> & dest) const
285 : {
286 146 : dest = *_diagonal;
287 142 : }
288 :
289 :
290 :
291 : template <typename T>
292 : void
293 71 : DiagonalMatrix<T>::get_transpose(SparseMatrix<T> & dest) const
294 : {
295 71 : auto diagonal_dest = dynamic_cast<DiagonalMatrix<T> *>(&dest);
296 71 : if (diagonal_dest)
297 4 : *diagonal_dest = *_diagonal;
298 : else
299 0 : libmesh_error_msg("DenseMatrix<T>::get_transpose currently only accepts another DenseMatrix<T> "
300 : "as its argument");
301 71 : }
302 :
303 :
304 :
305 : template <typename T>
306 : void
307 0 : DiagonalMatrix<T>::get_row(numeric_index_type i,
308 : std::vector<numeric_index_type> & indices,
309 : std::vector<T> & values) const
310 : {
311 0 : indices.clear();
312 0 : indices.push_back(i);
313 0 : values.clear();
314 0 : values.push_back((*_diagonal)(i));
315 0 : }
316 :
317 :
318 :
319 : template <typename T>
320 : void
321 71 : DiagonalMatrix<T>::zero_rows(std::vector<numeric_index_type> & rows, T val/*=0*/)
322 : {
323 524 : for (auto row : rows)
324 453 : _diagonal->set(row, val);
325 71 : }
326 :
327 : template <typename T>
328 : const NumericVector<T> &
329 83 : DiagonalMatrix<T>::diagonal() const
330 : {
331 83 : return *_diagonal;
332 : }
333 :
334 : template <typename T>
335 : void
336 0 : DiagonalMatrix<T>::restore_original_nonzero_pattern()
337 : {
338 0 : _diagonal->zero();
339 0 : }
340 :
341 : template class LIBMESH_EXPORT DiagonalMatrix<Number>;
342 : }
|
