libMesh
eigen_sparse_vector.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 // C++ includes
21 #include <algorithm> // for std::min
22 #include <limits>
23 
24 // Local Includes
25 #include "libmesh/dense_subvector.h"
26 #include "libmesh/dense_vector.h"
27 #include "libmesh/eigen_sparse_vector.h"
28 #include "libmesh/eigen_sparse_matrix.h"
29 #include "libmesh/int_range.h"
30 
31 #ifdef LIBMESH_HAVE_EIGEN
32 
33 namespace libMesh
34 {
35 
36 template <typename T>
38 {
39  libmesh_assert (this->closed());
40  libmesh_assert (this->initialized());
41 
42  return _vec.sum();
43 }
44 
45 
46 
47 template <typename T>
49 {
50  libmesh_assert (this->closed());
51  libmesh_assert (this->initialized());
52 
53  return _vec.lpNorm<1>();
54 }
55 
56 
57 
58 template <typename T>
60 {
61  libmesh_assert (this->closed());
62  libmesh_assert (this->initialized());
63 
64  return _vec.lpNorm<2>();
65 }
66 
67 
68 
69 template <typename T>
71 {
72  libmesh_assert (this->closed());
73  libmesh_assert (this->initialized());
74 
75  return _vec.lpNorm<Eigen::Infinity>();
76 }
77 
78 
79 
80 template <typename T>
82 {
83  libmesh_assert (this->closed());
84 
85  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
86 
87  _vec += v._vec;
88 
89  return *this;
90 }
91 
92 
93 
94 
95 template <typename T>
97 {
98  libmesh_assert (this->closed());
99 
100  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
101 
102  _vec -= v._vec;
103 
104  return *this;
105 }
106 
107 
108 
109 template <typename T>
111 {
112  libmesh_assert (this->closed());
113  libmesh_assert_equal_to(size(), v_in.size());
114 
115  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
116 
117  _vec = _vec.cwiseProduct(v._vec);
118 
119  return *this;
120 }
121 
122 
123 
124 template <typename T>
126 {
127  libmesh_assert (this->closed());
128  libmesh_assert_equal_to(size(), v_in.size());
129 
130  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
131 
132  _vec = _vec.cwiseQuotient(v._vec);
133 
134  return *this;
135 }
136 
137 
138 
139 
140 template <typename T>
142 {
143 #ifndef NDEBUG
144  const numeric_index_type n = this->size();
145 
146  for (numeric_index_type i=0; i<n; i++)
147  // Don't divide by zero!
148  libmesh_assert_not_equal_to ((*this)(i), T(0));
149 #endif
150 
151  _vec = _vec.cwiseInverse();
152 }
153 
154 
155 
156 template <typename T>
158 {
159  _vec = _vec.conjugate();
160 }
161 
162 
163 
164 template <typename T>
166 {
167  _vec += EigenSV::Constant(this->size(), v);
168 }
169 
170 
171 
172 
173 template <typename T>
175 {
176  libmesh_assert (this->initialized());
177 
178  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
179 
180  _vec += v._vec;
181 }
182 
183 
184 
185 template <typename T>
186 void EigenSparseVector<T>::add (const T a, const NumericVector<T> & v_in)
187 {
188  libmesh_assert (this->initialized());
189 
190  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
191 
192  _vec += v._vec*a;
193 }
194 
195 
196 
197 template <typename T>
199  const SparseMatrix<T> & mat_in)
200 {
201  // Make sure the data passed in are really in Eigen types
202  const EigenSparseVector<T> * e_vec = cast_ptr<const EigenSparseVector<T> *>(&vec_in);
203  const EigenSparseMatrix<T> * mat = cast_ptr<const EigenSparseMatrix<T> *>(&mat_in);
204 
205  libmesh_assert(e_vec);
206  libmesh_assert(mat);
207 
208  _vec += mat->_mat*e_vec->_vec;
209 }
210 
211 
212 
213 template <typename T>
215  const SparseMatrix<T> & mat_in)
216 {
217  // Make sure the data passed in are really in Eigen types
218  const EigenSparseVector<T> * e_vec = cast_ptr<const EigenSparseVector<T> *>(&vec_in);
219  const EigenSparseMatrix<T> * mat = cast_ptr<const EigenSparseMatrix<T> *>(&mat_in);
220 
221  libmesh_assert(e_vec);
222  libmesh_assert(mat);
223 
224  _vec += mat->_mat.transpose()*e_vec->_vec;
225 }
226 
227 
228 
229 template <typename T>
230 void EigenSparseVector<T>::scale (const T factor)
231 {
232  libmesh_assert (this->initialized());
233 
234  _vec *= factor;
235 }
236 
237 
238 
239 template <typename T>
241 {
242  libmesh_assert (this->initialized());
243 
244  const numeric_index_type n = this->size();
245 
246  for (numeric_index_type i=0; i!=n; ++i)
247  this->set(i,std::abs((*this)(i)));
248 }
249 
250 
251 
252 template <typename T>
254 {
255  libmesh_assert (this->initialized());
256 
257  // Make sure the NumericVector passed in is really a EigenSparseVector
258  const EigenSparseVector<T> * v = cast_ptr<const EigenSparseVector<T> *>(&v_in);
259  libmesh_assert(v);
260 
261  return _vec.dot(v->_vec);
262 }
263 
264 
265 
266 template <typename T>
269 {
270  libmesh_assert (this->initialized());
271  libmesh_assert (this->closed());
272 
273  _vec.fill(s);
274 
275  return *this;
276 }
277 
278 
279 
280 template <typename T>
283 {
284  // Make sure the NumericVector passed in is really a EigenSparseVector
285  const EigenSparseVector<T> * v =
286  cast_ptr<const EigenSparseVector<T> *>(&v_in);
287 
288  libmesh_assert(v);
289 
290  *this = *v;
291 
292  return *this;
293 }
294 
295 
296 
297 template <typename T>
300 {
301  libmesh_assert (this->initialized());
302  libmesh_assert (v.closed());
303  libmesh_assert_equal_to (this->size(), v.size());
304 
305  _vec = v._vec;
306 
307  this->_is_closed = true;
308 
309  return *this;
310 }
311 
312 
313 
314 template <typename T>
316 EigenSparseVector<T>::operator = (const std::vector<T> & v)
317 {
322  if (this->size() == v.size())
323  for (auto i : index_range(v))
324  this->set (i, v[i]);
325 
326  else
327  libmesh_error_msg("this->size() = " << this->size() << " must be equal to v.size() = " << v.size());
328 
329  return *this;
330 }
331 
332 
333 template <typename T>
335 {
336  // Make sure the NumericVector passed in is really a EigenSparseVector
337  EigenSparseVector<T> * v_local =
338  cast_ptr<EigenSparseVector<T> *>(&v_local_in);
339 
340  libmesh_assert(v_local);
341 
342  *v_local = *this;
343 }
344 
345 
346 
347 template <typename T>
349  const std::vector<numeric_index_type> & libmesh_dbg_var(send_list)) const
350 {
351  // Make sure the NumericVector passed in is really a EigenSparseVector
352  EigenSparseVector<T> * v_local =
353  cast_ptr<EigenSparseVector<T> *>(&v_local_in);
354 
355  libmesh_assert(v_local);
356  libmesh_assert_less_equal (send_list.size(), v_local->size());
357 
358  *v_local = *this;
359 }
360 
361 
362 
363 template <typename T>
364 void EigenSparseVector<T>::localize (std::vector<T> & v_local,
365  const std::vector<numeric_index_type> & indices) const
366 {
367  // EigenSparseVectors are serial, so we can just copy values
368  v_local.resize(indices.size());
369 
370  for (auto i : index_range(v_local))
371  v_local[i] = (*this)(indices[i]);
372 }
373 
374 
375 
376 template <typename T>
377 void EigenSparseVector<T>::localize (const numeric_index_type libmesh_dbg_var(first_local_idx),
378  const numeric_index_type libmesh_dbg_var(last_local_idx),
379  const std::vector<numeric_index_type> & libmesh_dbg_var(send_list))
380 {
381  libmesh_assert_equal_to (first_local_idx, 0);
382  libmesh_assert_equal_to (last_local_idx+1, this->size());
383 
384  libmesh_assert_less_equal (send_list.size(), this->size());
385 
386  this->_is_closed = true;
387 }
388 
389 
390 
391 template <typename T>
392 void EigenSparseVector<T>::localize (std::vector<T> & v_local) const
393 
394 {
395  v_local.resize(this->size());
396 
397  for (auto i : index_range(v_local))
398  v_local[i] = (*this)(i);
399 }
400 
401 
402 
403 template <typename T>
404 void EigenSparseVector<T>::localize_to_one (std::vector<T> & v_local,
405  const processor_id_type libmesh_dbg_var(pid)) const
406 {
407  libmesh_assert_equal_to (pid, 0);
408 
409  this->localize (v_local);
410 }
411 
412 
413 
414 template <typename T>
416  const NumericVector<T> & /*vec2*/)
417 {
418  libmesh_not_implemented();
419 }
420 
421 template <typename T>
423  const NumericVector<T> & /*vec2*/)
424 {
425  libmesh_not_implemented();
426 }
427 
428 
429 template <typename T>
431 {
432  libmesh_assert (this->initialized());
433  if (!this->size())
434  return -std::numeric_limits<Real>::max();
435 
436 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
437  Real the_max = libmesh_real((*this)(0));
438 
439  const numeric_index_type n = this->size();
440 
441  for (numeric_index_type i=1; i<n; i++)
442  the_max = std::max (the_max, libmesh_real((*this)(i)));
443 
444  return the_max;
445 #else
446  return libmesh_real(_vec.maxCoeff());
447 #endif
448 }
449 
450 
451 
452 template <typename T>
454 {
455  libmesh_assert (this->initialized());
456  if (!this->size())
457  return std::numeric_limits<Real>::max();
458 
459 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
460  Real the_min = libmesh_real((*this)(0));
461 
462  const numeric_index_type n = this->size();
463 
464  for (numeric_index_type i=1; i<n; i++)
465  the_min = std::min (the_min, libmesh_real((*this)(i)));
466 
467  return the_min;
468 #else
469  return libmesh_real(_vec.minCoeff());
470 #endif
471 }
472 
473 
474 //------------------------------------------------------------------
475 // Explicit instantiations
476 template class LIBMESH_EXPORT EigenSparseVector<Number>;
477 
478 } // namespace libMesh
479 
480 
481 #endif // #ifdef LIBMESH_HAVE_EIGEN
T libmesh_real(T a)
bool closed()
Checks that the library has been closed.
Definition: libmesh.C:273
The EigenSparseMatrix class wraps a sparse matrix object from the Eigen library.
virtual Real linfty_norm() const override
virtual void localize(std::vector< T > &v_local) const override
Creates a copy of the global vector in the local vector v_local.
virtual void conjugate() override
Negates the imaginary component of each entry in the vector.
virtual void pointwise_divide(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
virtual T sum() const override
virtual void add_vector(const NumericVector< T > &v, const SparseMatrix< T > &A) override
Computes , i.e.
DataType _vec
Actual Eigen::SparseVector<> we are wrapping.
virtual void reciprocal() override
Computes the component-wise reciprocal, .
virtual numeric_index_type size() const =0
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) override
Subtracts v from *this, .
virtual Real min() const override
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: vector_fe_ex5.C:43
The libMesh namespace provides an interface to certain functionality in the library.
uint8_t processor_id_type
Definition: id_types.h:104
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:57
Generic sparse matrix.
Definition: vector_fe_ex5.C:45
virtual Real l2_norm() const override
dof_id_type numeric_index_type
Definition: id_types.h:99
virtual Real max() const override
virtual NumericVector< T > & operator/=(const NumericVector< T > &v_in) override
Computes the component-wise division of this vector&#39;s entries by another&#39;s, .
virtual numeric_index_type size() const override
EigenSparseVector< T > & operator=(const EigenSparseVector< T > &v)
Copy assignment operator.
libmesh_assert(ctx)
virtual NumericVector< T > & operator*=(const NumericVector< T > &v_in) override
Computes the component-wise multiplication of this vector&#39;s entries by another&#39;s, ...
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A) override
Computes , i.e.
virtual void abs() override
Sets for each entry in the vector.
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const override
Creates a local copy of the global vector in v_local only on processor proc_id.
virtual bool closed() const
virtual void scale(const T factor) override
Scale each element of the vector by the given factor.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual Real l1_norm() const override
virtual void add(const numeric_index_type i, const T value) override
Adds value to the vector entry specified by i.
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:266
virtual T dot(const NumericVector< T > &v) const override
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
Adds v to *this, .
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:111
This class provides a nice interface to the Eigen C++-based data structures for serial vectors...