libMesh
numeric_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 #ifndef LIBMESH_NUMERIC_VECTOR_H
21 #define LIBMESH_NUMERIC_VECTOR_H
22 
23 // Local includes
24 #include "libmesh/libmesh_common.h"
25 #include "libmesh/auto_ptr.h"
26 #include "libmesh/enum_parallel_type.h"
27 #include "libmesh/enum_solver_package.h"
28 #include "libmesh/id_types.h"
29 #include "libmesh/reference_counted_object.h"
30 #include "libmesh/libmesh.h"
31 #include "libmesh/parallel_object.h"
32 #include "libmesh/dense_subvector.h"
33 #include "libmesh/dense_vector.h"
34 
35 // C++ includes
36 #include <cstddef>
37 #include <set>
38 #include <vector>
39 
40 namespace libMesh
41 {
42 
43 
44 // forward declarations
45 template <typename T> class NumericVector;
46 template <typename T> class DenseVector;
47 template <typename T> class DenseSubVector;
48 template <typename T> class SparseMatrix;
49 template <typename T> class ShellMatrix;
50 
51 
60 template <typename T>
61 class NumericVector : public ReferenceCountedObject<NumericVector<T>>,
62  public ParallelObject
63 {
64 public:
65 
69  explicit
70  NumericVector (const Parallel::Communicator & comm_in,
71  const ParallelType ptype = AUTOMATIC);
72 
76  explicit
77  NumericVector (const Parallel::Communicator & comm_in,
78  const numeric_index_type n,
79  const ParallelType ptype = AUTOMATIC);
80 
85  NumericVector (const Parallel::Communicator & comm_in,
86  const numeric_index_type n,
87  const numeric_index_type n_local,
88  const ParallelType ptype = AUTOMATIC);
89 
95  NumericVector (const Parallel::Communicator & comm_in,
96  const numeric_index_type N,
97  const numeric_index_type n_local,
98  const std::vector<numeric_index_type> & ghost,
99  const ParallelType ptype = AUTOMATIC);
100 
101 public:
102 
107  virtual ~NumericVector ();
108 
114  static UniquePtr<NumericVector<T>>
115  build(const Parallel::Communicator & comm,
116  const SolverPackage solver_package = libMesh::default_solver_package());
117 
118 #ifndef LIBMESH_DISABLE_COMMWORLD
119 
127 #ifdef LIBMESH_ENABLE_DEPRECATED
128  static UniquePtr<NumericVector<T>>
129  build(const SolverPackage solver_package = libMesh::default_solver_package());
130 #endif
131 #endif
132 
137  virtual bool initialized() const { return _is_initialized; }
138 
142  ParallelType type() const { return _type; }
143 
147  ParallelType & type() { return _type; }
148 
153  virtual bool closed() const { return _is_closed; }
154 
159  virtual void close () = 0;
160 
164  virtual void clear ();
165 
170  virtual void zero () = 0;
171 
178  virtual UniquePtr<NumericVector<T>> zero_clone () const = 0;
179 
185  virtual UniquePtr<NumericVector<T>> clone () const = 0;
186 
197  virtual void init (const numeric_index_type,
198  const numeric_index_type,
199  const bool = false,
200  const ParallelType = AUTOMATIC) = 0;
201 
205  virtual void init (const numeric_index_type,
206  const bool = false,
207  const ParallelType = AUTOMATIC) = 0;
208 
213  virtual void init (const numeric_index_type /*N*/,
214  const numeric_index_type /*n_local*/,
215  const std::vector<numeric_index_type> & /*ghost*/,
216  const bool /*fast*/ = false,
217  const ParallelType = AUTOMATIC) = 0;
218 
223  virtual void init (const NumericVector<T> & other,
224  const bool fast = false) = 0;
225 
231  virtual NumericVector<T> & operator= (const T s) = 0;
232 
238  virtual NumericVector<T> & operator= (const NumericVector<T> & v) = 0;
239 
245  virtual NumericVector<T> & operator= (const std::vector<T> & v) = 0;
246 
251  virtual Real min () const = 0;
252 
257  virtual Real max () const = 0;
258 
262  virtual T sum() const = 0;
263 
268  virtual Real l1_norm () const = 0;
269 
274  virtual Real l2_norm () const = 0;
275 
280  virtual Real linfty_norm () const = 0;
281 
288  virtual Real subset_l1_norm (const std::set<numeric_index_type> & indices) const;
289 
297  virtual Real subset_l2_norm (const std::set<numeric_index_type> & indices) const;
298 
305  virtual Real subset_linfty_norm (const std::set<numeric_index_type> & indices) const;
306 
310  virtual numeric_index_type size () const = 0;
311 
315  virtual numeric_index_type local_size() const = 0;
316 
323  virtual numeric_index_type first_local_index() const = 0;
324 
331  virtual numeric_index_type last_local_index() const = 0;
332 
336  virtual T operator() (const numeric_index_type i) const = 0;
337 
341  virtual T el(const numeric_index_type i) const { return (*this)(i); }
342 
349  virtual void get(const std::vector<numeric_index_type> & index,
350  T * values) const;
351 
358  void get(const std::vector<numeric_index_type> & index,
359  std::vector<T> & values) const;
360 
368  virtual NumericVector<T> & operator += (const NumericVector<T> & v) = 0;
369 
377  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) = 0;
378 
386  NumericVector<T> & operator *= (const T a) { this->scale(a); return *this; }
387 
395  NumericVector<T> & operator /= (const T a) { this->scale(1./a); return *this; }
396 
403  virtual NumericVector<T> & operator /= (NumericVector<T> & /*v*/) = 0;
404 
409  virtual void reciprocal() = 0;
410 
414  virtual void conjugate() = 0;
415 
419  virtual void set (const numeric_index_type i, const T value) = 0;
420 
424  virtual void add (const numeric_index_type i, const T value) = 0;
425 
430  virtual void add (const T s) = 0;
431 
437  virtual void add (const NumericVector<T> & v) = 0;
438 
444  virtual void add (const T a, const NumericVector<T> & v) = 0;
445 
452  virtual void add_vector (const T * v,
453  const std::vector<numeric_index_type> & dof_indices);
454 
460  void add_vector (const std::vector<T> & v,
461  const std::vector<numeric_index_type> & dof_indices);
462 
468  virtual void add_vector (const NumericVector<T> & v,
469  const std::vector<numeric_index_type> & dof_indices);
470 
476  void add_vector (const DenseVector<T> & v,
477  const std::vector<numeric_index_type> & dof_indices);
478 
484  virtual void add_vector (const NumericVector<T> & v,
485  const SparseMatrix<T> & A) = 0;
486 
492  void add_vector (const NumericVector<T> & v,
493  const ShellMatrix<T> & A);
494 
500  virtual void add_vector_transpose (const NumericVector<T> & v,
501  const SparseMatrix<T> & A) = 0;
502 
506  virtual void insert (const T * v,
507  const std::vector<numeric_index_type> & dof_indices);
508 
512  void insert (const std::vector<T> & v,
513  const std::vector<numeric_index_type> & dof_indices);
514 
518  virtual void insert (const NumericVector<T> & v,
519  const std::vector<numeric_index_type> & dof_indices);
520 
524  void insert (const DenseVector<T> & v,
525  const std::vector<numeric_index_type> & dof_indices);
526 
530  void insert (const DenseSubVector<T> & v,
531  const std::vector<numeric_index_type> & dof_indices);
532 
536  virtual void scale (const T factor) = 0;
537 
541  virtual void abs() = 0;
542 
549  virtual T dot(const NumericVector<T> & v) const = 0;
550 
555  virtual void localize (std::vector<T> & v_local) const = 0;
556 
561  virtual void localize (NumericVector<T> & v_local) const = 0;
562 
567  virtual void localize (NumericVector<T> & v_local,
568  const std::vector<numeric_index_type> & send_list) const = 0;
569 
599  virtual void localize (std::vector<T> & v_local,
600  const std::vector<numeric_index_type> & indices) const = 0;
601 
606  virtual void localize (const numeric_index_type first_local_idx,
607  const numeric_index_type last_local_idx,
608  const std::vector<numeric_index_type> & send_list) = 0;
609 
615  virtual void localize_to_one (std::vector<T> & v_local,
616  const processor_id_type proc_id=0) const = 0;
617 
623  virtual int compare (const NumericVector<T> & other_vector,
624  const Real threshold = TOLERANCE) const;
625 
631  virtual int local_relative_compare (const NumericVector<T> & other_vector,
632  const Real threshold = TOLERANCE) const;
633 
639  virtual int global_relative_compare (const NumericVector<T> & other_vector,
640  const Real threshold = TOLERANCE) const;
641 
647  virtual void pointwise_mult (const NumericVector<T> & vec1,
648  const NumericVector<T> & vec2) = 0;
649 
654  virtual void print(std::ostream & os=libMesh::out) const;
655 
660  virtual void print_global(std::ostream & os=libMesh::out) const;
661 
665  friend std::ostream & operator << (std::ostream & os, const NumericVector<T> & v)
666  {
667  v.print_global(os);
668  return os;
669  }
670 
676  virtual void print_matlab(const std::string & /*name*/ = "") const
677  {
678  libmesh_not_implemented();
679  }
680 
688  const std::vector<numeric_index_type> &) const
689  {
690  libmesh_not_implemented();
691  }
692 
698  virtual void swap (NumericVector<T> & v);
699 
700 protected:
701 
708 
713 
718 };
719 
720 
721 /*----------------------- Inline functions ----------------------------------*/
722 
723 
724 
725 template <typename T>
726 inline
728  const ParallelType ptype) :
729  ParallelObject(comm_in),
730  _is_closed(false),
731  _is_initialized(false),
732  _type(ptype)
733 {
734 }
735 
736 
737 
738 template <typename T>
739 inline
741  const numeric_index_type /*n*/,
742  const ParallelType ptype) :
743  ParallelObject(comm_in),
744  _is_closed(false),
745  _is_initialized(false),
746  _type(ptype)
747 {
748  libmesh_not_implemented(); // Abstract base class!
749  // init(n, n, false, ptype);
750 }
751 
752 
753 
754 template <typename T>
755 inline
757  const numeric_index_type /*n*/,
758  const numeric_index_type /*n_local*/,
759  const ParallelType ptype) :
760  ParallelObject(comm_in),
761  _is_closed(false),
762  _is_initialized(false),
763  _type(ptype)
764 {
765  libmesh_not_implemented(); // Abstract base class!
766  // init(n, n_local, false, ptype);
767 }
768 
769 
770 
771 template <typename T>
772 inline
774  const numeric_index_type /*n*/,
775  const numeric_index_type /*n_local*/,
776  const std::vector<numeric_index_type> & /*ghost*/,
777  const ParallelType ptype) :
778  ParallelObject(comm_in),
779  _is_closed(false),
780  _is_initialized(false),
781  _type(ptype)
782 {
783  libmesh_not_implemented(); // Abstract base class!
784  // init(n, n_local, ghost, false, ptype);
785 }
786 
787 
788 
789 template <typename T>
790 inline
792 {
793  clear ();
794 }
795 
796 
797 
798 template <typename T>
799 inline
801 {
802  _is_closed = false;
803  _is_initialized = false;
804 }
805 
806 
807 
808 template <typename T>
809 inline
810 void NumericVector<T>::get(const std::vector<numeric_index_type> & index,
811  T * values) const
812 {
813  const std::size_t num = index.size();
814  for (std::size_t i=0; i<num; i++)
815  {
816  values[i] = (*this)(index[i]);
817  }
818 }
819 
820 
821 
822 template <typename T>
823 inline
824 void NumericVector<T>::get(const std::vector<numeric_index_type> & index,
825  std::vector<T> & values) const
826 {
827  const std::size_t num = index.size();
828  values.resize(num);
829  if (!num)
830  return;
831 
832  this->get(index, &values[0]);
833 }
834 
835 
836 
837 template <typename T>
838 inline
839 void NumericVector<T>::add_vector(const std::vector<T> & v,
840  const std::vector<numeric_index_type> & dof_indices)
841 {
842  libmesh_assert(v.size() == dof_indices.size());
843  if (!v.empty())
844  this->add_vector(&v[0], dof_indices);
845 }
846 
847 
848 
849 template <typename T>
850 inline
852  const std::vector<numeric_index_type> & dof_indices)
853 {
854  libmesh_assert(v.size() == dof_indices.size());
855  if (!v.empty())
856  this->add_vector(&v(0), dof_indices);
857 }
858 
859 
860 
861 template <typename T>
862 inline
863 void NumericVector<T>::insert(const std::vector<T> & v,
864  const std::vector<numeric_index_type> & dof_indices)
865 {
866  libmesh_assert(v.size() == dof_indices.size());
867  if (!v.empty())
868  this->insert(&v[0], dof_indices);
869 }
870 
871 
872 
873 template <typename T>
874 inline
876  const std::vector<numeric_index_type> & dof_indices)
877 {
878  libmesh_assert(v.size() == dof_indices.size());
879  if (!v.empty())
880  this->insert(&v(0), dof_indices);
881 }
882 
883 
884 
885 template <typename T>
886 inline
888  const std::vector<numeric_index_type> & dof_indices)
889 {
890  libmesh_assert(v.size() == dof_indices.size());
891  if (!v.empty())
892  this->insert(&v(0), dof_indices);
893 }
894 
895 
896 
897 // Full specialization of the print() member for complex
898 // variables. This must precede the non-specialized
899 // version, at least according to icc v7.1
900 template <>
901 inline
902 void NumericVector<Complex>::print(std::ostream & os) const
903 {
904  libmesh_assert (this->initialized());
905  os << "Size\tglobal = " << this->size()
906  << "\t\tlocal = " << this->local_size() << std::endl;
907 
908  // std::complex<>::operator<<() is defined, but use this form
909  os << "#\tReal part\t\tImaginary part" << std::endl;
910  for (numeric_index_type i=this->first_local_index(); i<this->last_local_index(); i++)
911  os << i << "\t"
912  << (*this)(i).real() << "\t\t"
913  << (*this)(i).imag() << std::endl;
914 }
915 
916 
917 
918 template <typename T>
919 inline
920 void NumericVector<T>::print(std::ostream & os) const
921 {
922  libmesh_assert (this->initialized());
923  os << "Size\tglobal = " << this->size()
924  << "\t\tlocal = " << this->local_size() << std::endl;
925 
926  os << "#\tValue" << std::endl;
927  for (numeric_index_type i=this->first_local_index(); i<this->last_local_index(); i++)
928  os << i << "\t" << (*this)(i) << std::endl;
929 }
930 
931 
932 
933 template <>
934 inline
935 void NumericVector<Complex>::print_global(std::ostream & os) const
936 {
937  libmesh_assert (this->initialized());
938 
939  std::vector<Complex> v(this->size());
940  this->localize(v);
941 
942  // Right now we only want one copy of the output
943  if (this->processor_id())
944  return;
945 
946  os << "Size\tglobal = " << this->size() << std::endl;
947  os << "#\tReal part\t\tImaginary part" << std::endl;
948  for (numeric_index_type i=0; i!=v.size(); i++)
949  os << i << "\t"
950  << v[i].real() << "\t\t"
951  << v[i].imag() << std::endl;
952 }
953 
954 
955 template <typename T>
956 inline
957 void NumericVector<T>::print_global(std::ostream & os) const
958 {
959  libmesh_assert (this->initialized());
960 
961  std::vector<T> v(this->size());
962  this->localize(v);
963 
964  // Right now we only want one copy of the output
965  if (this->processor_id())
966  return;
967 
968  os << "Size\tglobal = " << this->size() << std::endl;
969  os << "#\tValue" << std::endl;
970  for (numeric_index_type i=0; i!=v.size(); i++)
971  os << i << "\t" << v[i] << std::endl;
972 }
973 
974 
975 
976 template <typename T>
977 inline
979 {
980  std::swap(_is_closed, v._is_closed);
981  std::swap(_is_initialized, v._is_initialized);
982  std::swap(_type, v._type);
983 }
984 
985 
986 } // namespace libMesh
987 
988 
989 #endif // LIBMESH_NUMERIC_VECTOR_H
virtual bool closed() const
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices)
Inserts the entries of v in *this at the locations specified by v.
virtual void create_subvector(NumericVector< T > &, const std::vector< numeric_index_type > &) const
Fills in subvector from this vector using the indices in rows.
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const =0
Creates a local copy of the global vector in v_local only on processor proc_id.
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
virtual int compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
virtual numeric_index_type size() const =0
virtual numeric_index_type last_local_index() const =0
virtual Real subset_linfty_norm(const std::set< numeric_index_type > &indices) const
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.
virtual int global_relative_compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
virtual Real max() const =0
virtual Real l2_norm() const =0
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A)=0
Computes , i.e.
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
SolverPackage
Defines an enum for various linear solver packages.
virtual T el(const numeric_index_type i) const
uint8_t processor_id_type
Definition: id_types.h:99
ParallelType & type()
Numeric vector.
Definition: dof_map.h:66
static const Real TOLERANCE
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.
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
virtual Real l1_norm() const =0
virtual bool empty() const libmesh_override
Definition: dense_vector.h:92
virtual T operator()(const numeric_index_type i) const =0
virtual void zero()=0
Set all entries to zero.
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual Real min() const =0
Defines a dense subvector for use in finite element computations.
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
Generic sparse matrix.
Definition: dof_map.h:65
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual void scale(const T factor)=0
Scale each element of the vector by the given factor.
virtual NumericVector< T > & operator+=(const NumericVector< T > &v)=0
Adds v to *this, .
virtual void print_matlab(const std::string &="") const
Print the contents of the vector in Matlab&#39;s sparse matrix format.
ParallelType
Defines an enum for parallel data structure types.
dof_id_type numeric_index_type
Definition: id_types.h:92
NumericVector(const Parallel::Communicator &comm_in, const ParallelType ptype=AUTOMATIC)
Dummy-Constructor.
NumericVector< T > & operator/=(const T a)
Scales the vector by 1/a, .
NumericVector< T > & operator*=(const T a)
Scales the vector by a, .
virtual Real subset_l1_norm(const std::set< numeric_index_type > &indices) const
virtual T sum() const =0
virtual int local_relative_compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
virtual void print_global(std::ostream &os=libMesh::out) const
Prints the global contents of the vector, by default to libMesh::out.
virtual NumericVector< T > & operator=(const T s)=0
Sets all entries of the vector to the value s.
This class forms the base class for all other classes that are expected to be implemented in parallel...
virtual Real linfty_norm() const =0
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
ParallelType _type
Type of vector.
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2)=0
Computes (summation not implied) i.e.
virtual numeric_index_type local_size() const =0
virtual void abs()=0
Sets for each entry in the vector.
virtual bool empty() const libmesh_override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual NumericVector< T > & operator-=(const NumericVector< T > &v)=0
Subtracts v from *this, .
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
static PetscErrorCode Mat * A
virtual UniquePtr< NumericVector< T > > zero_clone() const =0
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
const Parallel::Communicator & comm() const
OStreamProxy out
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
static const bool value
Definition: xdr_io.C:108
ParallelType type() const
virtual numeric_index_type first_local_index() const =0
virtual void print(std::ostream &os=libMesh::out) const
Prints the local contents of the vector, by default to libMesh::out.
virtual T dot(const NumericVector< T > &v) const =0
virtual void clear()
Restores the NumericVector<T> to a pristine state.
Defines a dense vector for use in Finite Element-type computations.
virtual void reciprocal()=0
Computes the pointwise reciprocal, .
Generic shell matrix, i.e.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
virtual ~NumericVector()
Destructor, deallocates memory.
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
virtual unsigned int size() const libmesh_override
virtual Real subset_l2_norm(const std::set< numeric_index_type > &indices) const
bool _is_closed
Flag which tracks whether the vector&#39;s values are consistent on all processors after insertion or add...
virtual void conjugate()=0
Negates the imaginary component of each entry in the vector.
processor_id_type processor_id() const
virtual UniquePtr< NumericVector< T > > clone() const =0