Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 :
21 : #ifndef LIBMESH_EIGEN_SPARSE_VECTOR_H
22 : #define LIBMESH_EIGEN_SPARSE_VECTOR_H
23 :
24 :
25 :
26 : #include "libmesh/libmesh_common.h"
27 :
28 : #ifdef LIBMESH_HAVE_EIGEN
29 :
30 : // Local includes
31 : #include "libmesh/eigen_core_support.h"
32 : #include "libmesh/numeric_vector.h"
33 :
34 : // C++ includes
35 : #include <limits>
36 : #include <mutex>
37 :
38 : namespace libMesh
39 : {
40 :
41 : // Forward declarations
42 : template <typename T> class EigenSparseMatrix;
43 : template <typename T> class EigenSparseLinearSolver;
44 : template <typename T> class SparseMatrix;
45 :
46 : /**
47 : * This class provides a nice interface to the Eigen C++-based data
48 : * structures for serial vectors. All overridden virtual functions are
49 : * documented in numeric_vector.h.
50 : *
51 : * \author Benjamin S. Kirk
52 : * \date 2002
53 : */
54 : template <typename T>
55 : class EigenSparseVector final : public NumericVector<T>
56 : {
57 : public:
58 :
59 : /**
60 : * Dummy-Constructor. Dimension=0
61 : */
62 : explicit
63 : EigenSparseVector (const Parallel::Communicator & comm_in,
64 : const ParallelType = AUTOMATIC);
65 :
66 : /**
67 : * Constructor. Set dimension to \p n and initialize all elements with zero.
68 : */
69 : explicit
70 : EigenSparseVector (const Parallel::Communicator & comm_in,
71 : const numeric_index_type n,
72 : const ParallelType = AUTOMATIC);
73 :
74 : /**
75 : * Constructor. Set local dimension to \p n_local, the global dimension
76 : * to \p n, and initialize all elements with zero.
77 : */
78 : EigenSparseVector (const Parallel::Communicator & comm_in,
79 : const numeric_index_type n,
80 : const numeric_index_type n_local,
81 : const ParallelType = AUTOMATIC);
82 :
83 : /**
84 : * Constructor. Set local dimension to \p n_local, the global
85 : * dimension to \p n, but additionally reserve memory for the
86 : * indices specified by the \p ghost argument.
87 : */
88 : EigenSparseVector (const Parallel::Communicator & comm_in,
89 : const numeric_index_type N,
90 : const numeric_index_type n_local,
91 : const std::vector<numeric_index_type> & ghost,
92 : const ParallelType = AUTOMATIC);
93 :
94 : /**
95 : * Copy assignment operator. Does some state checks before copying
96 : * the underlying Eigen data type.
97 : * \returns A reference to *this as the derived type.
98 : */
99 : EigenSparseVector<T> & operator= (const EigenSparseVector<T> & v);
100 :
101 : /**
102 : * The 5 special functions can be defaulted for this class, as it
103 : * does not manage any memory itself.
104 : */
105 : EigenSparseVector (EigenSparseVector &&) = default;
106 : EigenSparseVector (const EigenSparseVector &) = default;
107 : EigenSparseVector & operator= (EigenSparseVector &&) = default;
108 14248 : virtual ~EigenSparseVector () = default;
109 :
110 : /**
111 : * Convenient typedefs
112 : */
113 : typedef EigenSV DataType;
114 :
115 : virtual void close () override;
116 :
117 : virtual void clear () override;
118 :
119 : virtual void zero () override;
120 :
121 : virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
122 :
123 : virtual std::unique_ptr<NumericVector<T>> clone () const override;
124 :
125 : virtual void init (const numeric_index_type N,
126 : const numeric_index_type n_local,
127 : const bool fast=false,
128 : const ParallelType ptype=AUTOMATIC) override;
129 :
130 : virtual void init (const numeric_index_type N,
131 : const bool fast=false,
132 : const ParallelType ptype=AUTOMATIC) override;
133 :
134 : virtual void init (const numeric_index_type N,
135 : const numeric_index_type n_local,
136 : const std::vector<numeric_index_type> & ghost,
137 : const bool fast = false,
138 : const ParallelType = AUTOMATIC) override;
139 :
140 : virtual void init (const NumericVector<T> & other,
141 : const bool fast = false) override;
142 :
143 : virtual NumericVector<T> & operator= (const T s) override;
144 :
145 : virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
146 :
147 : virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
148 :
149 : virtual Real min () const override;
150 :
151 : virtual Real max () const override;
152 :
153 : virtual T sum () const override;
154 :
155 : virtual Real l1_norm () const override;
156 :
157 : virtual Real l2_norm () const override;
158 :
159 : virtual Real linfty_norm () const override;
160 :
161 : virtual numeric_index_type size () const override;
162 :
163 : virtual numeric_index_type local_size() const override;
164 :
165 : virtual numeric_index_type first_local_index() const override;
166 :
167 : virtual numeric_index_type last_local_index() const override;
168 :
169 : virtual T operator() (const numeric_index_type i) const override;
170 :
171 : virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
172 :
173 : virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
174 :
175 : virtual NumericVector<T> & operator *= (const NumericVector<T> & v_in) override;
176 :
177 : virtual NumericVector<T> & operator /= (const NumericVector<T> & v_in) override;
178 :
179 : virtual void reciprocal() override;
180 :
181 : virtual void conjugate() override;
182 :
183 : virtual void set (const numeric_index_type i, const T value) override;
184 :
185 : virtual void add (const numeric_index_type i, const T value) override;
186 :
187 : virtual void add (const T s) override;
188 :
189 : virtual void add (const NumericVector<T> & v) override;
190 :
191 : virtual void add (const T a, const NumericVector<T> & v) override;
192 :
193 : /**
194 : * We override one NumericVector<T>::add_vector() method but don't
195 : * want to hide the other defaults.
196 : */
197 : using NumericVector<T>::add_vector;
198 :
199 : virtual void add_vector (const NumericVector<T> & v,
200 : const SparseMatrix<T> & A) override;
201 :
202 : virtual void add_vector_transpose (const NumericVector<T> & v,
203 : const SparseMatrix<T> & A) override;
204 :
205 : virtual void scale (const T factor) override;
206 :
207 : virtual void abs() override;
208 :
209 : virtual T dot(const NumericVector<T> & v) const override;
210 :
211 : virtual void localize (std::vector<T> & v_local) const override;
212 :
213 : virtual void localize (NumericVector<T> & v_local) const override;
214 :
215 : virtual void localize (NumericVector<T> & v_local,
216 : const std::vector<numeric_index_type> & send_list) const override;
217 :
218 : virtual void localize (std::vector<T> & v_local,
219 : const std::vector<numeric_index_type> & indices) const override;
220 :
221 : virtual void localize (const numeric_index_type first_local_idx,
222 : const numeric_index_type last_local_idx,
223 : const std::vector<numeric_index_type> & send_list) override;
224 :
225 : virtual void localize_to_one (std::vector<T> & v_local,
226 : const processor_id_type proc_id=0) const override;
227 :
228 : virtual void pointwise_mult (const NumericVector<T> & vec1,
229 : const NumericVector<T> & vec2) override;
230 :
231 : virtual void pointwise_divide (const NumericVector<T> & vec1,
232 : const NumericVector<T> & vec2) override;
233 :
234 : virtual void create_subvector(NumericVector<T> & subvector,
235 : const std::vector<numeric_index_type> & rows,
236 : bool supplying_global_rows = true) const override;
237 :
238 : virtual void swap (NumericVector<T> & v) override;
239 :
240 : virtual std::size_t max_allowed_id() const override;
241 :
242 : virtual std::unique_ptr<NumericVector<T>>
243 : get_subvector(const std::vector<numeric_index_type> & rows) override;
244 :
245 : virtual void restore_subvector(std::unique_ptr<NumericVector<T>> subvector,
246 : const std::vector<numeric_index_type> & rows) override;
247 :
248 : /**
249 : * References to the underlying Eigen data types.
250 : *
251 : * \note This is generally not required in user-level code.
252 : */
253 158 : DataType & vec () { return _vec; }
254 426 : const DataType & vec () const { return _vec; }
255 :
256 : private:
257 :
258 : /**
259 : * Actual Eigen::SparseVector<> we are wrapping.
260 : */
261 : DataType _vec;
262 :
263 : /**
264 : * Make other Eigen datatypes friends
265 : */
266 : friend class EigenSparseMatrix<T>;
267 : friend class EigenSparseLinearSolver<T>;
268 : };
269 :
270 :
271 :
272 : // ---------------------------------------------------------
273 : // EigenSparseVector inline methods
274 : template <typename T>
275 : inline
276 13780 : EigenSparseVector<T>::EigenSparseVector (const Parallel::Communicator & comm_in,
277 : const ParallelType ptype)
278 13780 : : NumericVector<T>(comm_in, ptype)
279 : {
280 13780 : this->_type = ptype;
281 0 : }
282 :
283 :
284 :
285 : template <typename T>
286 : inline
287 142 : EigenSparseVector<T>::EigenSparseVector (const Parallel::Communicator & comm_in,
288 : const numeric_index_type n,
289 : const ParallelType ptype)
290 142 : : NumericVector<T>(comm_in, ptype)
291 : {
292 142 : this->init(n, n, false, ptype);
293 142 : }
294 :
295 :
296 :
297 : template <typename T>
298 : inline
299 0 : EigenSparseVector<T>::EigenSparseVector (const Parallel::Communicator & comm_in,
300 : const numeric_index_type n,
301 : const numeric_index_type n_local,
302 : const ParallelType ptype)
303 0 : : NumericVector<T>(comm_in, ptype)
304 : {
305 0 : this->init(n, n_local, false, ptype);
306 0 : }
307 :
308 :
309 :
310 : template <typename T>
311 : inline
312 0 : EigenSparseVector<T>::EigenSparseVector (const Parallel::Communicator & comm_in,
313 : const numeric_index_type N,
314 : const numeric_index_type n_local,
315 : const std::vector<numeric_index_type> & ghost,
316 : const ParallelType ptype)
317 0 : : NumericVector<T>(comm_in, ptype)
318 : {
319 0 : this->init(N, n_local, ghost, false, ptype);
320 0 : }
321 :
322 :
323 :
324 : template <typename T>
325 : inline
326 8090 : void EigenSparseVector<T>::init (const numeric_index_type n,
327 : const numeric_index_type n_local,
328 : const bool fast,
329 : const ParallelType)
330 : {
331 : // Eigen vectors only for serial cases,
332 : // but can provide a "parallel" vector on one processor.
333 8090 : libmesh_error_msg_if(n != n_local, "Error: EigenSparseVectors can only be used in serial!");
334 :
335 8090 : this->_type = SERIAL;
336 :
337 : // Clear initialized vectors
338 8090 : if (this->initialized())
339 728 : this->clear();
340 :
341 8090 : _vec.resize(n);
342 :
343 8090 : this->_is_initialized = true;
344 8090 : this->_is_closed = true;
345 :
346 : // Optionally zero out all components
347 8090 : if (fast == false)
348 0 : this->zero ();
349 :
350 8090 : return;
351 : }
352 :
353 :
354 :
355 : template <typename T>
356 : inline
357 605 : void EigenSparseVector<T>::init (const numeric_index_type n,
358 : const bool fast,
359 : const ParallelType ptype)
360 : {
361 605 : this->init(n,n,fast,ptype);
362 0 : }
363 :
364 :
365 : template <typename T>
366 : inline
367 886 : void EigenSparseVector<T>::init (const numeric_index_type n,
368 : const numeric_index_type n_local,
369 : const std::vector<numeric_index_type> & libmesh_dbg_var(ghost),
370 : const bool fast,
371 : const ParallelType ptype)
372 : {
373 0 : libmesh_assert(ghost.empty());
374 886 : this->init(n,n_local,fast,ptype);
375 0 : }
376 :
377 :
378 :
379 : /* Default implementation for solver packages for which ghosted
380 : vectors are not yet implemented. */
381 : template <class T>
382 830 : void EigenSparseVector<T>::init (const NumericVector<T> & other,
383 : const bool fast)
384 : {
385 830 : this->init(other.size(),other.local_size(),fast,other.type());
386 830 : }
387 :
388 :
389 :
390 : template <typename T>
391 : inline
392 12715 : void EigenSparseVector<T>::close ()
393 : {
394 0 : libmesh_assert (this->initialized());
395 :
396 16399 : this->_is_closed = true;
397 12715 : }
398 :
399 :
400 :
401 : template <typename T>
402 : inline
403 1364 : void EigenSparseVector<T>::clear ()
404 : {
405 0 : _vec.resize(0);
406 :
407 1364 : this->_is_initialized = false;
408 1364 : this->_is_closed = false;
409 1364 : }
410 :
411 :
412 :
413 : template <typename T> inline
414 33743 : void EigenSparseVector<T>::zero ()
415 : {
416 0 : libmesh_assert (this->initialized());
417 0 : libmesh_assert (this->closed());
418 :
419 39859 : _vec.setZero();
420 6116 : }
421 :
422 :
423 :
424 : template <typename T>
425 : inline
426 102 : std::unique_ptr<NumericVector<T>> EigenSparseVector<T>::zero_clone () const
427 : {
428 102 : std::unique_ptr<NumericVector<T>> cloned_vector =
429 : std::make_unique<EigenSparseVector<T>>(this->comm());
430 102 : cloned_vector->init(*this);
431 102 : return cloned_vector;
432 0 : }
433 :
434 :
435 :
436 : template <typename T>
437 : inline
438 716 : std::unique_ptr<NumericVector<T>> EigenSparseVector<T>::clone () const
439 : {
440 716 : std::unique_ptr<NumericVector<T>> cloned_vector =
441 : std::make_unique<EigenSparseVector<T>>(this->comm());
442 716 : cloned_vector->init(*this, true);
443 716 : *cloned_vector = *this;
444 716 : return cloned_vector;
445 0 : }
446 :
447 :
448 :
449 : template <typename T>
450 : inline
451 2057 : numeric_index_type EigenSparseVector<T>::size () const
452 : {
453 0 : libmesh_assert (this->initialized());
454 :
455 528 : return static_cast<numeric_index_type>(_vec.size());
456 : }
457 :
458 :
459 :
460 : template <typename T>
461 : inline
462 1041 : numeric_index_type EigenSparseVector<T>::local_size () const
463 : {
464 0 : libmesh_assert (this->initialized());
465 :
466 1041 : return this->size();
467 : }
468 :
469 :
470 :
471 : template <typename T>
472 : inline
473 334455 : numeric_index_type EigenSparseVector<T>::first_local_index () const
474 : {
475 0 : libmesh_assert (this->initialized());
476 :
477 334455 : return 0;
478 : }
479 :
480 :
481 :
482 : template <typename T>
483 : inline
484 2123787 : numeric_index_type EigenSparseVector<T>::last_local_index () const
485 : {
486 0 : libmesh_assert (this->initialized());
487 :
488 2123787 : return this->size();
489 : }
490 :
491 :
492 :
493 : template <typename T>
494 : inline
495 3772826 : void EigenSparseVector<T>::set (const numeric_index_type i, const T value)
496 : {
497 0 : libmesh_assert (this->initialized());
498 0 : libmesh_assert_less (i, this->size());
499 :
500 3772826 : std::scoped_lock lock(this->_numeric_vector_mutex);
501 3772826 : _vec[static_cast<eigen_idx_type>(i)] = value;
502 :
503 3772826 : this->_is_closed = false;
504 3772826 : }
505 :
506 :
507 :
508 : template <typename T>
509 : inline
510 21841284 : void EigenSparseVector<T>::add (const numeric_index_type i, const T value)
511 : {
512 0 : libmesh_assert (this->initialized());
513 0 : libmesh_assert_less (i, this->size());
514 :
515 21841284 : std::scoped_lock lock(this->_numeric_vector_mutex);
516 21841284 : _vec[static_cast<eigen_idx_type>(i)] += value;
517 :
518 21841284 : this->_is_closed = false;
519 21841284 : }
520 :
521 :
522 :
523 : template <typename T>
524 : inline
525 150097576 : T EigenSparseVector<T>::operator() (const numeric_index_type i) const
526 : {
527 0 : libmesh_assert (this->initialized());
528 0 : libmesh_assert ( ((i >= this->first_local_index()) &&
529 : (i < this->last_local_index())) );
530 :
531 3039163 : return _vec[static_cast<eigen_idx_type>(i)];
532 : }
533 :
534 :
535 :
536 : template <typename T>
537 : inline
538 226 : void EigenSparseVector<T>::swap (NumericVector<T> & other)
539 : {
540 0 : EigenSparseVector<T> & v = cast_ref<EigenSparseVector<T> &>(other);
541 :
542 0 : _vec.swap(v._vec);
543 :
544 0 : std::swap (this->_is_closed, v._is_closed);
545 0 : std::swap (this->_is_initialized, v._is_initialized);
546 0 : std::swap (this->_type, v._type);
547 226 : }
548 :
549 :
550 :
551 : template <typename T>
552 : inline
553 321 : std::size_t EigenSparseVector<T>::max_allowed_id () const
554 : {
555 : // We use the Eigen::Matrix type which appears to be templated on
556 : // int for its sizes, see e.g. https://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html
557 321 : return std::numeric_limits<int>::max();
558 : }
559 :
560 : } // namespace libMesh
561 :
562 :
563 : #endif // #ifdef LIBMESH_HAVE_EIGEN
564 : #endif // LIBMESH_EIGEN_SPARSE_VECTOR_H
|
