libMesh
eigen_sparse_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 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 namespace libMesh
35 {
36 
37 // Forward declarations
38 template <typename T> class EigenSparseMatrix;
39 template <typename T> class EigenSparseLinearSolver;
40 template <typename T> class SparseMatrix;
41 
50 template <typename T>
51 class EigenSparseVector final : public NumericVector<T>
52 {
53 public:
54 
58  explicit
59  EigenSparseVector (const Parallel::Communicator & comm_in,
60  const ParallelType = AUTOMATIC);
61 
65  explicit
66  EigenSparseVector (const Parallel::Communicator & comm_in,
67  const numeric_index_type n,
68  const ParallelType = AUTOMATIC);
69 
74  EigenSparseVector (const Parallel::Communicator & comm_in,
75  const numeric_index_type n,
76  const numeric_index_type n_local,
77  const ParallelType = AUTOMATIC);
78 
84  EigenSparseVector (const Parallel::Communicator & comm_in,
85  const numeric_index_type N,
86  const numeric_index_type n_local,
87  const std::vector<numeric_index_type> & ghost,
88  const ParallelType = AUTOMATIC);
89 
95  EigenSparseVector<T> & operator= (const EigenSparseVector<T> & v);
96 
101  EigenSparseVector (EigenSparseVector &&) = default;
102  EigenSparseVector (const EigenSparseVector &) = default;
104  virtual ~EigenSparseVector () = default;
105 
109  typedef EigenSV DataType;
110 
111  virtual void close () override;
112 
113  virtual void clear () override;
114 
115  virtual void zero () override;
116 
117  virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
118 
119  virtual std::unique_ptr<NumericVector<T>> clone () const override;
120 
121  virtual void init (const numeric_index_type N,
122  const numeric_index_type n_local,
123  const bool fast=false,
124  const ParallelType ptype=AUTOMATIC) override;
125 
126  virtual void init (const numeric_index_type N,
127  const bool fast=false,
128  const ParallelType ptype=AUTOMATIC) override;
129 
130  virtual void init (const numeric_index_type N,
131  const numeric_index_type n_local,
132  const std::vector<numeric_index_type> & ghost,
133  const bool fast = false,
134  const ParallelType = AUTOMATIC) override;
135 
136  virtual void init (const NumericVector<T> & other,
137  const bool fast = false) override;
138 
139  virtual NumericVector<T> & operator= (const T s) override;
140 
141  virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
142 
143  virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
144 
145  virtual Real min () const override;
146 
147  virtual Real max () const override;
148 
149  virtual T sum () const override;
150 
151  virtual Real l1_norm () const override;
152 
153  virtual Real l2_norm () const override;
154 
155  virtual Real linfty_norm () const override;
156 
157  virtual numeric_index_type size () const override;
158 
159  virtual numeric_index_type local_size() const override;
160 
161  virtual numeric_index_type first_local_index() const override;
162 
163  virtual numeric_index_type last_local_index() const override;
164 
165  virtual T operator() (const numeric_index_type i) const override;
166 
167  virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
168 
169  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
170 
171  virtual NumericVector<T> & operator /= (const NumericVector<T> & v_in) override;
172 
173  virtual void reciprocal() override;
174 
175  virtual void conjugate() override;
176 
177  virtual void set (const numeric_index_type i, const T value) override;
178 
179  virtual void add (const numeric_index_type i, const T value) override;
180 
181  virtual void add (const T s) override;
182 
183  virtual void add (const NumericVector<T> & v) override;
184 
185  virtual void add (const T a, const NumericVector<T> & v) override;
186 
192 
193  virtual void add_vector (const NumericVector<T> & v,
194  const SparseMatrix<T> & A) override;
195 
196  virtual void add_vector_transpose (const NumericVector<T> & v,
197  const SparseMatrix<T> & A) override;
198 
199  virtual void scale (const T factor) override;
200 
201  virtual void abs() override;
202 
203  virtual T dot(const NumericVector<T> & v) const override;
204 
205  virtual void localize (std::vector<T> & v_local) const override;
206 
207  virtual void localize (NumericVector<T> & v_local) const override;
208 
209  virtual void localize (NumericVector<T> & v_local,
210  const std::vector<numeric_index_type> & send_list) const override;
211 
212  virtual void localize (std::vector<T> & v_local,
213  const std::vector<numeric_index_type> & indices) const override;
214 
215  virtual void localize (const numeric_index_type first_local_idx,
216  const numeric_index_type last_local_idx,
217  const std::vector<numeric_index_type> & send_list) override;
218 
219  virtual void localize_to_one (std::vector<T> & v_local,
220  const processor_id_type proc_id=0) const override;
221 
222  virtual void pointwise_mult (const NumericVector<T> & vec1,
223  const NumericVector<T> & vec2) override;
224 
225  virtual void swap (NumericVector<T> & v) override;
226 
232  DataType & vec () { return _vec; }
233  const DataType & vec () const { return _vec; }
234 
235 private:
236 
241 
245  friend class EigenSparseMatrix<T>;
246  friend class EigenSparseLinearSolver<T>;
247 };
248 
249 
250 
251 // ---------------------------------------------------------
252 // EigenSparseVector inline methods
253 template <typename T>
254 inline
256  const ParallelType ptype)
257  : NumericVector<T>(comm_in, ptype)
258 {
259  this->_type = ptype;
260 }
261 
262 
263 
264 template <typename T>
265 inline
267  const numeric_index_type n,
268  const ParallelType ptype)
269  : NumericVector<T>(comm_in, ptype)
270 {
271  this->init(n, n, false, ptype);
272 }
273 
274 
275 
276 template <typename T>
277 inline
279  const numeric_index_type n,
280  const numeric_index_type n_local,
281  const ParallelType ptype)
282  : NumericVector<T>(comm_in, ptype)
283 {
284  this->init(n, n_local, false, ptype);
285 }
286 
287 
288 
289 template <typename T>
290 inline
292  const numeric_index_type N,
293  const numeric_index_type n_local,
294  const std::vector<numeric_index_type> & ghost,
295  const ParallelType ptype)
296  : NumericVector<T>(comm_in, ptype)
297 {
298  this->init(N, n_local, ghost, false, ptype);
299 }
300 
301 
302 
303 template <typename T>
304 inline
306  const numeric_index_type n_local,
307  const bool fast,
308  const ParallelType)
309 {
310  // Eigen vectors only for serial cases,
311  // but can provide a "parallel" vector on one processor.
312  if (n != n_local)
313  libmesh_error_msg("Error: EigenSparseVectors can only be used in serial!");
314 
315  this->_type = SERIAL;
316 
317  // Clear initialized vectors
318  if (this->initialized())
319  this->clear();
320 
321  _vec.resize(n);
322 
323  this->_is_initialized = true;
324 #ifndef NDEBUG
325  this->_is_closed = true;
326 #endif
327 
328  // Optionally zero out all components
329  if (fast == false)
330  this->zero ();
331 
332  return;
333 }
334 
335 
336 
337 template <typename T>
338 inline
340  const bool fast,
341  const ParallelType ptype)
342 {
343  this->init(n,n,fast,ptype);
344 }
345 
346 
347 template <typename T>
348 inline
350  const numeric_index_type n_local,
351  const std::vector<numeric_index_type> & libmesh_dbg_var(ghost),
352  const bool fast,
353  const ParallelType ptype)
354 {
355  libmesh_assert(ghost.empty());
356  this->init(n,n_local,fast,ptype);
357 }
358 
359 
360 
361 /* Default implementation for solver packages for which ghosted
362  vectors are not yet implemented. */
363 template <class T>
365  const bool fast)
366 {
367  this->init(other.size(),other.local_size(),fast,other.type());
368 }
369 
370 
371 
372 template <typename T>
373 inline
375 {
376  libmesh_assert (this->initialized());
377 
378 #ifndef NDEBUG
379  this->_is_closed = true;
380 #endif
381 }
382 
383 
384 
385 template <typename T>
386 inline
388 {
389  _vec.resize(0);
390 
391  this->_is_initialized = false;
392 #ifndef NDEBUG
393  this->_is_closed = false;
394 #endif
395 }
396 
397 
398 
399 template <typename T> inline
401 {
402  libmesh_assert (this->initialized());
403  libmesh_assert (this->closed());
404 
405  _vec.setZero();
406 }
407 
408 
409 
410 template <typename T>
411 inline
412 std::unique_ptr<NumericVector<T>> EigenSparseVector<T>::zero_clone () const
413 {
414  NumericVector<T> * cloned_vector = new EigenSparseVector<T>(this->comm());
415  cloned_vector->init(*this);
416  return std::unique_ptr<NumericVector<T>>(cloned_vector);
417 }
418 
419 
420 
421 template <typename T>
422 inline
423 std::unique_ptr<NumericVector<T>> EigenSparseVector<T>::clone () const
424 {
425  NumericVector<T> * cloned_vector = new EigenSparseVector<T>(this->comm());
426  cloned_vector->init(*this, true);
427  *cloned_vector = *this;
428  return std::unique_ptr<NumericVector<T>>(cloned_vector);
429 }
430 
431 
432 
433 template <typename T>
434 inline
436 {
437  libmesh_assert (this->initialized());
438 
439  return static_cast<numeric_index_type>(_vec.size());
440 }
441 
442 
443 
444 template <typename T>
445 inline
447 {
448  libmesh_assert (this->initialized());
449 
450  return this->size();
451 }
452 
453 
454 
455 template <typename T>
456 inline
458 {
459  libmesh_assert (this->initialized());
460 
461  return 0;
462 }
463 
464 
465 
466 template <typename T>
467 inline
469 {
470  libmesh_assert (this->initialized());
471 
472  return this->size();
473 }
474 
475 
476 
477 template <typename T>
478 inline
479 void EigenSparseVector<T>::set (const numeric_index_type i, const T value)
480 {
481  libmesh_assert (this->initialized());
482  libmesh_assert_less (i, this->size());
483 
484  _vec[static_cast<eigen_idx_type>(i)] = value;
485 
486 #ifndef NDEBUG
487  this->_is_closed = false;
488 #endif
489 }
490 
491 
492 
493 template <typename T>
494 inline
495 void EigenSparseVector<T>::add (const numeric_index_type i, const T value)
496 {
497  libmesh_assert (this->initialized());
498  libmesh_assert_less (i, this->size());
499 
500  _vec[static_cast<eigen_idx_type>(i)] += value;
501 
502 #ifndef NDEBUG
503  this->_is_closed = false;
504 #endif
505 }
506 
507 
508 
509 template <typename T>
510 inline
512 {
513  libmesh_assert (this->initialized());
514  libmesh_assert ( ((i >= this->first_local_index()) &&
515  (i < this->last_local_index())) );
516 
517  return _vec[static_cast<eigen_idx_type>(i)];
518 }
519 
520 
521 
522 template <typename T>
523 inline
525 {
526  EigenSparseVector<T> & v = cast_ref<EigenSparseVector<T> &>(other);
527 
528  _vec.swap(v._vec);
529 
530  std::swap (this->_is_closed, v._is_closed);
531  std::swap (this->_is_initialized, v._is_initialized);
532  std::swap (this->_type, v._type);
533 }
534 
535 
536 } // namespace libMesh
537 
538 
539 #endif // #ifdef LIBMESH_HAVE_EIGEN
540 #endif // LIBMESH_EIGEN_SPARSE_VECTOR_H
bool closed()
Checks that the library has been closed.
The EigenSparseMatrix class wraps a sparse matrix object from the Eigen library.
virtual NumericVector< T > & operator/=(const NumericVector< T > &v_in) override
Computes the pointwise division of this vector&#39;s entries by another&#39;s, .
Encapsulates the MPI_Comm object.
Definition: communicator.h:92
virtual Real min() const override
virtual void abs() override
Sets for each entry in the vector.
DataType _vec
Actual Eigen::SparseVector<> we are wrapping.
virtual void add_vector(const NumericVector< T > &v, const SparseMatrix< T > &A) override
Computes , i.e.
virtual numeric_index_type size() const =0
This class provides an interface to Eigen iterative solvers that is compatible with the libMesh Linea...
virtual std::unique_ptr< NumericVector< T > > clone() const override
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: dof_map.h:76
virtual numeric_index_type last_local_index() const override
virtual void conjugate() override
Negates the imaginary component of each entry in the vector.
int32_t eigen_idx_type
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:239
virtual void reciprocal() override
Computes the pointwise reciprocal, .
virtual void close() override
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
Change the dimension of the vector to n.
virtual void swap(NumericVector< T > &v) override
Swaps the contents of this with v.
virtual std::unique_ptr< NumericVector< T > > zero_clone() const override
uint8_t processor_id_type
Definition: id_types.h:104
virtual Real l1_norm() const override
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.
Generic sparse matrix.
Definition: dof_map.h:75
EigenSV DataType
Convenient typedefs.
dof_id_type numeric_index_type
Definition: id_types.h:99
virtual void clear() override
Restores the NumericVector<T> to a pristine state.
virtual void localize(std::vector< T > &v_local) const override
Creates a copy of the global vector in the local vector v_local.
Eigen::Matrix< Number, Eigen::Dynamic, 1 > EigenSV
virtual Real linfty_norm() const override
virtual numeric_index_type size() const override
const DataType & vec() const
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
virtual numeric_index_type local_size() const override
virtual Real max() const override
virtual void init(const numeric_index_type N, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC) override
Change the dimension of the vector to n.
DataType & vec()
References to the underlying Eigen data types.
virtual T dot(const NumericVector< T > &v) const override
ParallelType _type
Type of vector.
virtual Real l2_norm() const override
virtual void set(const numeric_index_type i, const T value) override
Sets v(i) = value.
EigenSparseVector< T > & operator=(const EigenSparseVector< T > &v)
Copy assignment operator.
virtual numeric_index_type first_local_index() const override
ParallelType type() 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 void add(const numeric_index_type i, const T value) override
Adds value to each entry of the vector.
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
virtual numeric_index_type local_size() const =0
virtual void zero() override
Set all entries to zero.
bool initialized()
Checks that library initialization has been done.
EigenSparseVector(const Parallel::Communicator &comm_in, const ParallelType=AUTOMATIC)
Dummy-Constructor.
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
Adds v to *this, .
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) override
Subtracts v from *this, .
virtual T sum() const override
virtual T operator()(const numeric_index_type i) const override
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A) override
Computes , i.e.
virtual ~EigenSparseVector()=default
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
This class provides a nice interface to the Eigen C++-based data structures for serial vectors...
ParallelType
Defines an enum for parallel data structure types.