libMesh
eigen_sparse_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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>
52 {
53 public:
54 
58  explicit
60  const ParallelType = AUTOMATIC);
61 
65  explicit
67  const numeric_index_type n,
68  const ParallelType = AUTOMATIC);
69 
75  const numeric_index_type n,
76  const numeric_index_type n_local,
77  const ParallelType = AUTOMATIC);
78 
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 
99  typedef EigenSV DataType;
100 
101  virtual void close () libmesh_override;
102 
103  virtual void clear () libmesh_override;
104 
105  virtual void zero () libmesh_override;
106 
107  virtual UniquePtr<NumericVector<T>> zero_clone () const libmesh_override;
108 
109  virtual UniquePtr<NumericVector<T>> clone () const libmesh_override;
110 
111  virtual void init (const numeric_index_type N,
112  const numeric_index_type n_local,
113  const bool fast=false,
114  const ParallelType ptype=AUTOMATIC) libmesh_override;
115 
116  virtual void init (const numeric_index_type N,
117  const bool fast=false,
118  const ParallelType ptype=AUTOMATIC) libmesh_override;
119 
120  virtual void init (const numeric_index_type N,
121  const numeric_index_type n_local,
122  const std::vector<numeric_index_type> & ghost,
123  const bool fast = false,
124  const ParallelType = AUTOMATIC) libmesh_override;
125 
126  virtual void init (const NumericVector<T> & other,
127  const bool fast = false) libmesh_override;
128 
129  virtual NumericVector<T> & operator= (const T s) libmesh_override;
130 
131  virtual NumericVector<T> & operator= (const NumericVector<T> & v) libmesh_override;
132 
138  EigenSparseVector<T> & operator= (const EigenSparseVector<T> & v);
139 
140  virtual NumericVector<T> & operator= (const std::vector<T> & v) libmesh_override;
141 
142  virtual Real min () const libmesh_override;
143 
144  virtual Real max () const libmesh_override;
145 
146  virtual T sum () const libmesh_override;
147 
148  virtual Real l1_norm () const libmesh_override;
149 
150  virtual Real l2_norm () const libmesh_override;
151 
152  virtual Real linfty_norm () const libmesh_override;
153 
154  virtual numeric_index_type size () const libmesh_override;
155 
156  virtual numeric_index_type local_size() const libmesh_override;
157 
158  virtual numeric_index_type first_local_index() const libmesh_override;
159 
160  virtual numeric_index_type last_local_index() const libmesh_override;
161 
162  virtual T operator() (const numeric_index_type i) const libmesh_override;
163 
164  virtual NumericVector<T> & operator += (const NumericVector<T> & v) libmesh_override;
165 
166  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) libmesh_override;
167 
168  virtual NumericVector<T> & operator /= (NumericVector<T> & v_in) libmesh_override;
169 
170  virtual void reciprocal() libmesh_override;
171 
172  virtual void conjugate() libmesh_override;
173 
174  virtual void set (const numeric_index_type i, const T value) libmesh_override;
175 
176  virtual void add (const numeric_index_type i, const T value) libmesh_override;
177 
178  virtual void add (const T s) libmesh_override;
179 
180  virtual void add (const NumericVector<T> & v) libmesh_override;
181 
182  virtual void add (const T a, const NumericVector<T> & v) libmesh_override;
183 
188  using NumericVector<T>::add_vector;
189 
190  virtual void add_vector (const NumericVector<T> & v,
191  const SparseMatrix<T> & A) libmesh_override;
192 
193  virtual void add_vector_transpose (const NumericVector<T> & v,
194  const SparseMatrix<T> & A) libmesh_override;
195 
196  virtual void scale (const T factor) libmesh_override;
197 
198  virtual void abs() libmesh_override;
199 
200  virtual T dot(const NumericVector<T> & v) const libmesh_override;
201 
202  virtual void localize (std::vector<T> & v_local) const libmesh_override;
203 
204  virtual void localize (NumericVector<T> & v_local) const libmesh_override;
205 
206  virtual void localize (NumericVector<T> & v_local,
207  const std::vector<numeric_index_type> & send_list) const libmesh_override;
208 
209  virtual void localize (std::vector<T> & v_local,
210  const std::vector<numeric_index_type> & indices) const libmesh_override;
211 
212  virtual void localize (const numeric_index_type first_local_idx,
213  const numeric_index_type last_local_idx,
214  const std::vector<numeric_index_type> & send_list) libmesh_override;
215 
216  virtual void localize_to_one (std::vector<T> & v_local,
217  const processor_id_type proc_id=0) const libmesh_override;
218 
219  virtual void pointwise_mult (const NumericVector<T> & vec1,
220  const NumericVector<T> & vec2) libmesh_override;
221 
222  virtual void swap (NumericVector<T> & v) libmesh_override;
223 
229  DataType & vec () { return _vec; }
230  const DataType & vec () const { return _vec; }
231 
232 private:
233 
237  DataType _vec;
238 
242  friend class EigenSparseMatrix<T>;
243  friend class EigenSparseLinearSolver<T>;
244 };
245 
246 
247 
248 // ---------------------------------------------------------
249 // EigenSparseVector inline methods
250 template <typename T>
251 inline
253  const ParallelType ptype)
254  : NumericVector<T>(comm_in, ptype)
255 {
256  this->_type = ptype;
257 }
258 
259 
260 
261 template <typename T>
262 inline
264  const numeric_index_type n,
265  const ParallelType ptype)
266  : NumericVector<T>(comm_in, ptype)
267 {
268  this->init(n, n, false, ptype);
269 }
270 
271 
272 
273 template <typename T>
274 inline
276  const numeric_index_type n,
277  const numeric_index_type n_local,
278  const ParallelType ptype)
279  : NumericVector<T>(comm_in, ptype)
280 {
281  this->init(n, n_local, false, ptype);
282 }
283 
284 
285 
286 template <typename T>
287 inline
289  const numeric_index_type N,
290  const numeric_index_type n_local,
291  const std::vector<numeric_index_type> & ghost,
292  const ParallelType ptype)
293  : NumericVector<T>(comm_in, ptype)
294 {
295  this->init(N, n_local, ghost, false, ptype);
296 }
297 
298 
299 
300 template <typename T>
301 inline
303 {
304  this->clear ();
305 }
306 
307 
308 
309 template <typename T>
310 inline
312  const numeric_index_type n_local,
313  const bool fast,
314  const ParallelType)
315 {
316  // Eigen vectors only for serial cases,
317  // but can provide a "parallel" vector on one processor.
318  if (n != n_local)
319  libmesh_error_msg("Error: EigenSparseVectors can only be used in serial!");
320 
321  this->_type = SERIAL;
322 
323  // Clear initialized vectors
324  if (this->initialized())
325  this->clear();
326 
327  _vec.resize(n);
328 
329  this->_is_initialized = true;
330 #ifndef NDEBUG
331  this->_is_closed = true;
332 #endif
333 
334  // Optionally zero out all components
335  if (fast == false)
336  this->zero ();
337 
338  return;
339 }
340 
341 
342 
343 template <typename T>
344 inline
346  const bool fast,
347  const ParallelType ptype)
348 {
349  this->init(n,n,fast,ptype);
350 }
351 
352 
353 template <typename T>
354 inline
356  const numeric_index_type n_local,
357  const std::vector<numeric_index_type> & libmesh_dbg_var(ghost),
358  const bool fast,
359  const ParallelType ptype)
360 {
361  libmesh_assert(ghost.empty());
362  this->init(n,n_local,fast,ptype);
363 }
364 
365 
366 
367 /* Default implementation for solver packages for which ghosted
368  vectors are not yet implemented. */
369 template <class T>
370 void EigenSparseVector<T>::init (const NumericVector<T> & other,
371  const bool fast)
372 {
373  this->init(other.size(),other.local_size(),fast,other.type());
374 }
375 
376 
377 
378 template <typename T>
379 inline
381 {
382  libmesh_assert (this->initialized());
383 
384 #ifndef NDEBUG
385  this->_is_closed = true;
386 #endif
387 }
388 
389 
390 
391 template <typename T>
392 inline
394 {
395  _vec.resize(0);
396 
397  this->_is_initialized = false;
398 #ifndef NDEBUG
399  this->_is_closed = false;
400 #endif
401 }
402 
403 
404 
405 template <typename T> inline
407 {
408  libmesh_assert (this->initialized());
409  libmesh_assert (this->closed());
410 
411  _vec.setZero();
412 }
413 
414 
415 
416 template <typename T>
417 inline
419 {
420  NumericVector<T> * cloned_vector = new EigenSparseVector<T>(this->comm());
421  cloned_vector->init(*this);
422  return UniquePtr<NumericVector<T>>(cloned_vector);
423 }
424 
425 
426 
427 template <typename T>
428 inline
430 {
431  NumericVector<T> * cloned_vector = new EigenSparseVector<T>(this->comm());
432  cloned_vector->init(*this, true);
433  *cloned_vector = *this;
434  return UniquePtr<NumericVector<T>>(cloned_vector);
435 }
436 
437 
438 
439 template <typename T>
440 inline
442 {
443  libmesh_assert (this->initialized());
444 
445  return static_cast<numeric_index_type>(_vec.size());
446 }
447 
448 
449 
450 template <typename T>
451 inline
453 {
454  libmesh_assert (this->initialized());
455 
456  return this->size();
457 }
458 
459 
460 
461 template <typename T>
462 inline
464 {
465  libmesh_assert (this->initialized());
466 
467  return 0;
468 }
469 
470 
471 
472 template <typename T>
473 inline
475 {
476  libmesh_assert (this->initialized());
477 
478  return this->size();
479 }
480 
481 
482 
483 template <typename T>
484 inline
485 void EigenSparseVector<T>::set (const numeric_index_type i, const T value)
486 {
487  libmesh_assert (this->initialized());
488  libmesh_assert_less (i, this->size());
489 
490  _vec[static_cast<eigen_idx_type>(i)] = value;
491 
492 #ifndef NDEBUG
493  this->_is_closed = false;
494 #endif
495 }
496 
497 
498 
499 template <typename T>
500 inline
501 void EigenSparseVector<T>::add (const numeric_index_type i, const T value)
502 {
503  libmesh_assert (this->initialized());
504  libmesh_assert_less (i, this->size());
505 
506  _vec[static_cast<eigen_idx_type>(i)] += value;
507 
508 #ifndef NDEBUG
509  this->_is_closed = false;
510 #endif
511 }
512 
513 
514 
515 template <typename T>
516 inline
518 {
519  libmesh_assert (this->initialized());
520  libmesh_assert ( ((i >= this->first_local_index()) &&
521  (i < this->last_local_index())) );
522 
523  return _vec[static_cast<eigen_idx_type>(i)];
524 }
525 
526 
527 
528 template <typename T>
529 inline
531 {
532  EigenSparseVector<T> & v = cast_ref<EigenSparseVector<T> &>(other);
533 
534  _vec.swap(v._vec);
535 
536  std::swap (this->_is_closed, v._is_closed);
537  std::swap (this->_is_initialized, v._is_initialized);
538  std::swap (this->_type, v._type);
539 }
540 
541 
542 } // namespace libMesh
543 
544 
545 #endif // #ifdef LIBMESH_HAVE_EIGEN
546 #endif // LIBMESH_EIGEN_SPARSE_VECTOR_H
virtual numeric_index_type n() const libmesh_override
double abs(double a)
EigenSV DataType
Convenient typedefs.
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
The Hex20 is an element in 3D composed of 20 nodes.
Definition: cell_hex20.h:66
virtual numeric_index_type size() const =0
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
uint8_t processor_id_type
Definition: id_types.h:99
This class provides an interface to Eigen iterative solvers that is compatible with the libMesh Linea...
Numeric vector.
Definition: dof_map.h:66
int32_t eigen_idx_type
bool _is_initialized
true once init() has been called.
virtual bool initialized() const
The libMesh namespace provides an interface to certain functionality in the library.
void sum(T &r, const Communicator &comm=Communicator_World)
const Number zero
.
Definition: libmesh.h:178
virtual numeric_index_type size() const libmesh_override
long double max(long double a, double b)
libmesh_assert(j)
virtual bool closed() const libmesh_override
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Generic sparse matrix.
Definition: dof_map.h:65
DataType _vec
Actual Eigen::SparseVector<> we are wrapping.
ParallelType
Defines an enum for parallel data structure types.
dof_id_type numeric_index_type
Definition: id_types.h:92
Eigen::Matrix< Number, Eigen::Dynamic, 1 > EigenSV
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
PetscErrorCode Vec Mat libmesh_dbg_var(j)
virtual void init() libmesh_override
Initialize this matrix using the sparsity structure computed by dof_map.
ParallelType _type
Type of vector.
virtual numeric_index_type local_size() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void clear() libmesh_override
Restores the NumericVector<T> to a pristine state.
X_input swap(X_system)
static PetscErrorCode Mat * A
const Parallel::Communicator & comm() const
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
virtual void zero() libmesh_override
Set all entries to zero.
virtual numeric_index_type last_local_index() const libmesh_override
static const bool value
Definition: xdr_io.C:108
const DataType & vec() const
ParallelType type() const
virtual numeric_index_type first_local_index() const libmesh_override
long double min(long double a, double b)
bool _is_closed
Flag which tracks whether the vector&#39;s values are consistent on all processors after insertion or add...