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 :
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 14390 : 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 swap (NumericVector<T> & v) override;
235 :
236 : virtual std::size_t max_allowed_id() const override;
237 :
238 : /**
239 : * References to the underlying Eigen data types.
240 : *
241 : * \note This is generally not required in user-level code.
242 : */
243 0 : DataType & vec () { return _vec; }
244 0 : const DataType & vec () const { return _vec; }
245 :
246 : private:
247 :
248 : /**
249 : * Actual Eigen::SparseVector<> we are wrapping.
250 : */
251 : DataType _vec;
252 :
253 : /**
254 : * Make other Eigen datatypes friends
255 : */
256 : friend class EigenSparseMatrix<T>;
257 : friend class EigenSparseLinearSolver<T>;
258 : };
259 :
260 :
261 :
262 : // ---------------------------------------------------------
263 : // EigenSparseVector inline methods
264 : template <typename T>
265 : inline
266 13566 : EigenSparseVector<T>::EigenSparseVector (const Parallel::Communicator & comm_in,
267 : const ParallelType ptype)
268 13566 : : NumericVector<T>(comm_in, ptype)
269 : {
270 13566 : this->_type = ptype;
271 0 : }
272 :
273 :
274 :
275 : template <typename T>
276 : inline
277 0 : EigenSparseVector<T>::EigenSparseVector (const Parallel::Communicator & comm_in,
278 : const numeric_index_type n,
279 : const ParallelType ptype)
280 0 : : NumericVector<T>(comm_in, ptype)
281 : {
282 0 : this->init(n, n, false, ptype);
283 0 : }
284 :
285 :
286 :
287 : template <typename T>
288 : inline
289 0 : EigenSparseVector<T>::EigenSparseVector (const Parallel::Communicator & comm_in,
290 : const numeric_index_type n,
291 : const numeric_index_type n_local,
292 : const ParallelType ptype)
293 0 : : NumericVector<T>(comm_in, ptype)
294 : {
295 0 : this->init(n, n_local, false, ptype);
296 0 : }
297 :
298 :
299 :
300 : template <typename T>
301 : inline
302 0 : EigenSparseVector<T>::EigenSparseVector (const Parallel::Communicator & comm_in,
303 : const numeric_index_type N,
304 : const numeric_index_type n_local,
305 : const std::vector<numeric_index_type> & ghost,
306 : const ParallelType ptype)
307 0 : : NumericVector<T>(comm_in, ptype)
308 : {
309 0 : this->init(N, n_local, ghost, false, ptype);
310 0 : }
311 :
312 :
313 :
314 : template <typename T>
315 : inline
316 8150 : void EigenSparseVector<T>::init (const numeric_index_type n,
317 : const numeric_index_type n_local,
318 : const bool fast,
319 : const ParallelType)
320 : {
321 : // Eigen vectors only for serial cases,
322 : // but can provide a "parallel" vector on one processor.
323 8150 : libmesh_error_msg_if(n != n_local, "Error: EigenSparseVectors can only be used in serial!");
324 :
325 8150 : this->_type = SERIAL;
326 :
327 : // Clear initialized vectors
328 8150 : if (this->initialized())
329 728 : this->clear();
330 :
331 8150 : _vec.resize(n);
332 :
333 8150 : this->_is_initialized = true;
334 8150 : this->_is_closed = true;
335 :
336 : // Optionally zero out all components
337 8150 : if (fast == false)
338 0 : this->zero ();
339 :
340 8150 : return;
341 : }
342 :
343 :
344 :
345 : template <typename T>
346 : inline
347 605 : void EigenSparseVector<T>::init (const numeric_index_type n,
348 : const bool fast,
349 : const ParallelType ptype)
350 : {
351 605 : this->init(n,n,fast,ptype);
352 605 : }
353 :
354 :
355 : template <typename T>
356 : inline
357 897 : void EigenSparseVector<T>::init (const numeric_index_type n,
358 : const numeric_index_type n_local,
359 : const std::vector<numeric_index_type> & libmesh_dbg_var(ghost),
360 : const bool fast,
361 : const ParallelType ptype)
362 : {
363 0 : libmesh_assert(ghost.empty());
364 897 : this->init(n,n_local,fast,ptype);
365 0 : }
366 :
367 :
368 :
369 : /* Default implementation for solver packages for which ghosted
370 : vectors are not yet implemented. */
371 : template <class T>
372 835 : void EigenSparseVector<T>::init (const NumericVector<T> & other,
373 : const bool fast)
374 : {
375 835 : this->init(other.size(),other.local_size(),fast,other.type());
376 835 : }
377 :
378 :
379 :
380 : template <typename T>
381 : inline
382 12778 : void EigenSparseVector<T>::close ()
383 : {
384 0 : libmesh_assert (this->initialized());
385 :
386 16387 : this->_is_closed = true;
387 12778 : }
388 :
389 :
390 :
391 : template <typename T>
392 : inline
393 1364 : void EigenSparseVector<T>::clear ()
394 : {
395 0 : _vec.resize(0);
396 :
397 1364 : this->_is_initialized = false;
398 1364 : this->_is_closed = false;
399 1364 : }
400 :
401 :
402 :
403 : template <typename T> inline
404 33764 : void EigenSparseVector<T>::zero ()
405 : {
406 0 : libmesh_assert (this->initialized());
407 0 : libmesh_assert (this->closed());
408 :
409 39940 : _vec.setZero();
410 6176 : }
411 :
412 :
413 :
414 : template <typename T>
415 : inline
416 107 : std::unique_ptr<NumericVector<T>> EigenSparseVector<T>::zero_clone () const
417 : {
418 107 : std::unique_ptr<NumericVector<T>> cloned_vector =
419 : std::make_unique<EigenSparseVector<T>>(this->comm());
420 107 : cloned_vector->init(*this);
421 107 : return cloned_vector;
422 0 : }
423 :
424 :
425 :
426 : template <typename T>
427 : inline
428 716 : std::unique_ptr<NumericVector<T>> EigenSparseVector<T>::clone () const
429 : {
430 716 : std::unique_ptr<NumericVector<T>> cloned_vector =
431 : std::make_unique<EigenSparseVector<T>>(this->comm());
432 716 : cloned_vector->init(*this, true);
433 716 : *cloned_vector = *this;
434 716 : return cloned_vector;
435 0 : }
436 :
437 :
438 :
439 : template <typename T>
440 : inline
441 2062 : numeric_index_type EigenSparseVector<T>::size () const
442 : {
443 0 : libmesh_assert (this->initialized());
444 :
445 528 : return static_cast<numeric_index_type>(_vec.size());
446 : }
447 :
448 :
449 :
450 : template <typename T>
451 : inline
452 1046 : numeric_index_type EigenSparseVector<T>::local_size () const
453 : {
454 0 : libmesh_assert (this->initialized());
455 :
456 1046 : return this->size();
457 : }
458 :
459 :
460 :
461 : template <typename T>
462 : inline
463 334455 : numeric_index_type EigenSparseVector<T>::first_local_index () const
464 : {
465 0 : libmesh_assert (this->initialized());
466 :
467 334455 : return 0;
468 : }
469 :
470 :
471 :
472 : template <typename T>
473 : inline
474 2123787 : numeric_index_type EigenSparseVector<T>::last_local_index () const
475 : {
476 0 : libmesh_assert (this->initialized());
477 :
478 2123787 : return this->size();
479 : }
480 :
481 :
482 :
483 : template <typename T>
484 : inline
485 3772826 : void EigenSparseVector<T>::set (const numeric_index_type i, const T value)
486 : {
487 0 : libmesh_assert (this->initialized());
488 0 : libmesh_assert_less (i, this->size());
489 :
490 3772826 : std::scoped_lock lock(this->_numeric_vector_mutex);
491 3772826 : _vec[static_cast<eigen_idx_type>(i)] = value;
492 :
493 3772826 : this->_is_closed = false;
494 3772826 : }
495 :
496 :
497 :
498 : template <typename T>
499 : inline
500 23004380 : void EigenSparseVector<T>::add (const numeric_index_type i, const T value)
501 : {
502 0 : libmesh_assert (this->initialized());
503 0 : libmesh_assert_less (i, this->size());
504 :
505 23004380 : std::scoped_lock lock(this->_numeric_vector_mutex);
506 23004380 : _vec[static_cast<eigen_idx_type>(i)] += value;
507 :
508 23004380 : this->_is_closed = false;
509 23004380 : }
510 :
511 :
512 :
513 : template <typename T>
514 : inline
515 154044059 : T EigenSparseVector<T>::operator() (const numeric_index_type i) const
516 : {
517 0 : libmesh_assert (this->initialized());
518 0 : libmesh_assert ( ((i >= this->first_local_index()) &&
519 : (i < this->last_local_index())) );
520 :
521 3210244 : return _vec[static_cast<eigen_idx_type>(i)];
522 : }
523 :
524 :
525 :
526 : template <typename T>
527 : inline
528 226 : void EigenSparseVector<T>::swap (NumericVector<T> & other)
529 : {
530 0 : EigenSparseVector<T> & v = cast_ref<EigenSparseVector<T> &>(other);
531 :
532 0 : _vec.swap(v._vec);
533 :
534 0 : std::swap (this->_is_closed, v._is_closed);
535 0 : std::swap (this->_is_initialized, v._is_initialized);
536 0 : std::swap (this->_type, v._type);
537 226 : }
538 :
539 :
540 :
541 : template <typename T>
542 : inline
543 332 : std::size_t EigenSparseVector<T>::max_allowed_id () const
544 : {
545 : // We use the Eigen::Matrix type which appears to be templated on
546 : // int for its sizes, see e.g. https://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html
547 332 : return std::numeric_limits<int>::max();
548 : }
549 :
550 : } // namespace libMesh
551 :
552 :
553 : #endif // #ifdef LIBMESH_HAVE_EIGEN
554 : #endif // LIBMESH_EIGEN_SPARSE_VECTOR_H
|
