eigen_sparse_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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
32 #include "libmesh/numeric_vector.h"
33 
34 // C++ includes
35 
36 // Eigen includes
37 
38 namespace libMesh
39 {
40 
41 
42 // Forward declarations
43 template <typename T> class EigenSparseMatrix;
44 template <typename T> class EigenSparseLinearSolver;
45 template <typename T> class SparseMatrix;
46 
47 
55 template <typename T>
56 class EigenSparseVector : public NumericVector<T>
57 {
58  public:
59 
63  explicit
64  EigenSparseVector (const Parallel::Communicator &comm,
65  const ParallelType = AUTOMATIC);
66 
70  explicit
71  EigenSparseVector (const Parallel::Communicator &comm,
72  const numeric_index_type n,
73  const ParallelType = AUTOMATIC);
74 
79  EigenSparseVector (const Parallel::Communicator &comm,
80  const numeric_index_type n,
82  const ParallelType = AUTOMATIC);
83 
89  EigenSparseVector (const Parallel::Communicator &comm,
90  const numeric_index_type N,
92  const std::vector<numeric_index_type>& ghost,
93  const ParallelType = AUTOMATIC);
94 
100 
104  typedef EigenSV DataType;
105 
109  void close ();
110 
114  void clear ();
115 
120  void zero ();
121 
127  virtual AutoPtr<NumericVector<T> > zero_clone () const;
128 
132  AutoPtr<NumericVector<T> > clone () const;
133 
147  void init (const numeric_index_type N,
149  const bool fast=false,
150  const ParallelType ptype=AUTOMATIC);
151 
155  void init (const numeric_index_type N,
156  const bool fast=false,
157  const ParallelType ptype=AUTOMATIC);
158 
163  void init (const numeric_index_type /*N*/,
164  const numeric_index_type /*n_local*/,
165  const std::vector<numeric_index_type>& /*ghost*/,
166  const bool /*fast*/ = false,
167  const ParallelType = AUTOMATIC);
168 
173  virtual void init (const NumericVector<T>& other,
174  const bool fast = false);
175 
179  NumericVector<T> & operator= (const T s);
180 
185 
190 
194  NumericVector<T> & operator= (const std::vector<T> &v);
195 
201  Real min () const;
202 
208  Real max () const;
209 
213  T sum () const;
214 
219  Real l1_norm () const;
220 
226  Real l2_norm () const;
227 
233  Real linfty_norm () const;
234 
242  numeric_index_type size () const;
243 
249 
255 
261 
265  T operator() (const numeric_index_type i) const;
266 
272 
278 
283 
287  virtual void reciprocal();
288 
293  virtual void conjugate();
294 
298  void set (const numeric_index_type i, const T value);
299 
303  void add (const numeric_index_type i, const T value);
304 
310  void add (const T s);
311 
317  void add (const NumericVector<T>& V);
318 
324  void add (const T a, const NumericVector<T>& v);
325 
331  void add_vector (const std::vector<T>& v,
332  const std::vector<numeric_index_type>& dof_indices);
333 
340  void add_vector (const NumericVector<T>& V,
341  const std::vector<numeric_index_type>& dof_indices);
342 
347  void add_vector (const NumericVector<T> &,
348  const SparseMatrix<T> &);
349 
356  void add_vector (const DenseVector<T>& V,
357  const std::vector<numeric_index_type>& dof_indices);
358 
364  const SparseMatrix<T> &);
365 
370  virtual void insert (const std::vector<T>& v,
371  const std::vector<numeric_index_type>& dof_indices);
372 
379  virtual void insert (const NumericVector<T>& V,
380  const std::vector<numeric_index_type>& dof_indices);
381 
387  virtual void insert (const DenseVector<T>& V,
388  const std::vector<numeric_index_type>& dof_indices);
389 
395  virtual void insert (const DenseSubVector<T>& V,
396  const std::vector<numeric_index_type>& dof_indices);
397 
402  void scale (const T factor);
403 
408  virtual void abs();
409 
413  virtual T dot(const NumericVector<T>& V) const;
414 
419  void localize (std::vector<T>& v_local) const;
420 
425  void localize (NumericVector<T>& v_local) const;
426 
432  void localize (NumericVector<T>& v_local,
433  const std::vector<numeric_index_type>& send_list) const;
434 
439  void localize (const numeric_index_type first_local_idx,
440  const numeric_index_type last_local_idx,
441  const std::vector<numeric_index_type>& send_list);
442 
443 
450  void localize_to_one (std::vector<T>& v_local,
451  const processor_id_type proc_id=0) const;
452 
457  virtual void pointwise_mult (const NumericVector<T>& vec1,
458  const NumericVector<T>& vec2);
459 
463  virtual void swap (NumericVector<T> &v);
464 
469  DataType & vec () { return _vec; }
470  const DataType & vec () const { return _vec; }
471 
472  private:
473 
478 
482  friend class EigenSparseMatrix<T>;
483  friend class EigenSparseLinearSolver<T>;
484 };
485 
486 
487 
488 //----------------------- ----------------------------------
489 // EigenSparseVector inline methods
490 template <typename T>
491 inline
493  const ParallelType ptype)
494  : NumericVector<T>(comm, ptype)
495 {
496  this->_type = ptype;
497 }
498 
499 
500 
501 template <typename T>
502 inline
504  const numeric_index_type n,
505  const ParallelType ptype)
506  : NumericVector<T>(comm, ptype)
507 {
508  this->init(n, n, false, ptype);
509 }
510 
511 
512 
513 template <typename T>
514 inline
516  const numeric_index_type n,
518  const ParallelType ptype)
519  : NumericVector<T>(comm, ptype)
520 {
521  this->init(n, n_local, false, ptype);
522 }
523 
524 
525 
526 template <typename T>
527 inline
529  const numeric_index_type N,
531  const std::vector<numeric_index_type>& ghost,
532  const ParallelType ptype)
533  : NumericVector<T>(comm, ptype)
534 {
535  this->init(N, n_local, ghost, false, ptype);
536 }
537 
538 
539 
540 template <typename T>
541 inline
543 {
544  this->clear ();
545 }
546 
547 
548 
549 template <typename T>
550 inline
552  const numeric_index_type libmesh_dbg_var(n_local),
553  const bool fast,
554  const ParallelType)
555 {
556  // Eigen vectors only for serial cases,
557  // but can provide a "parallel" vector on one processor.
558  libmesh_assert_equal_to (n, n_local);
559 
560  this->_type = SERIAL;
561 
562  // Clear initialized vectors
563  if (this->initialized())
564  this->clear();
565 
566  _vec.resize(n);
567 
568  this->_is_initialized = true;
569 #ifndef NDEBUG
570  this->_is_closed = true;
571 #endif
572 
573  // Optionally zero out all components
574  if (fast == false)
575  this->zero ();
576 
577  return;
578 }
579 
580 
581 
582 template <typename T>
583 inline
585  const bool fast,
586  const ParallelType ptype)
587 {
588  this->init(n,n,fast,ptype);
589 }
590 
591 
592 template <typename T>
593 inline
596  const std::vector<numeric_index_type>& libmesh_dbg_var(ghost),
597  const bool fast,
598  const ParallelType ptype)
599 {
600  libmesh_assert(ghost.empty());
601  this->init(n,n_local,fast,ptype);
602 }
603 
604 
605 
606 /* Default implementation for solver packages for which ghosted
607  vectors are not yet implemented. */
608 template <class T>
610  const bool fast)
611 {
612  this->init(other.size(),other.local_size(),fast,other.type());
613 }
614 
615 
616 
617 template <typename T>
618 inline
620 {
621  libmesh_assert (this->initialized());
622 
623 #ifndef NDEBUG
624  this->_is_closed = true;
625 #endif
626 }
627 
628 
629 
630 template <typename T>
631 inline
633 {
634  _vec.resize(0);
635 
636  this->_is_initialized = false;
637 #ifndef NDEBUG
638  this->_is_closed = false;
639 #endif
640 }
641 
642 
643 
644 template <typename T> inline
646 {
647  libmesh_assert (this->initialized());
648  libmesh_assert (this->closed());
649 
650  _vec.setZero();
651 }
652 
653 
654 
655 template <typename T>
656 inline
658 {
659  AutoPtr<NumericVector<T> > cloned_vector
660  (new EigenSparseVector<T>(this->comm()));
661 
662  cloned_vector->init(*this);
663 
664  return cloned_vector;
665 }
666 
667 
668 
669 template <typename T>
670 inline
672 {
673  AutoPtr<NumericVector<T> > cloned_vector
674  (new EigenSparseVector<T>(this->comm()));
675 
676  cloned_vector->init(*this, true);
677 
678  *cloned_vector = *this;
679 
680  return cloned_vector;
681 }
682 
683 
684 
685 template <typename T>
686 inline
688 {
689  libmesh_assert (this->initialized());
690 
691  return static_cast<numeric_index_type>(_vec.size());
692 }
693 
694 
695 
696 template <typename T>
697 inline
699 {
700  libmesh_assert (this->initialized());
701 
702  return this->size();
703 }
704 
705 
706 
707 template <typename T>
708 inline
710 {
711  libmesh_assert (this->initialized());
712 
713  return 0;
714 }
715 
716 
717 
718 template <typename T>
719 inline
721 {
722  libmesh_assert (this->initialized());
723 
724  return this->size();
725 }
726 
727 
728 
729 template <typename T>
730 inline
731 void EigenSparseVector<T>::set (const numeric_index_type i, const T value)
732 {
733  libmesh_assert (this->initialized());
734  libmesh_assert_less (i, this->size());
735 
736  _vec[static_cast<eigen_idx_type>(i)] = value;
737 
738 #ifndef NDEBUG
739  this->_is_closed = false;
740 #endif
741 }
742 
743 
744 
745 template <typename T>
746 inline
747 void EigenSparseVector<T>::add (const numeric_index_type i, const T value)
748 {
749  libmesh_assert (this->initialized());
750  libmesh_assert_less (i, this->size());
751 
752  _vec[static_cast<eigen_idx_type>(i)] += value;
753 
754 #ifndef NDEBUG
755  this->_is_closed = false;
756 #endif
757 }
758 
759 
760 
761 template <typename T>
762 inline
764 {
765  libmesh_assert (this->initialized());
766  libmesh_assert ( ((i >= this->first_local_index()) &&
767  (i < this->last_local_index())) );
768 
769  return _vec[static_cast<eigen_idx_type>(i)];
770 }
771 
772 
773 
774 template <typename T>
775 inline
777 {
778  EigenSparseVector<T>& v = libmesh_cast_ref<EigenSparseVector<T>&>(other);
779 
780  _vec.swap(v._vec);
781 
782  std::swap (this->_is_closed, v._is_closed);
783  std::swap (this->_is_initialized, v._is_initialized);
784  std::swap (this->_type, v._type);
785 }
786 
787 
788 } // namespace libMesh
789 
790 
791 #endif // #ifdef LIBMESH_HAVE_EIGEN
792 #endif // LIBMESH_EIGEN_SPARSE_VECTOR_H

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo