numeric_vector.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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" // deprecated
27 #include "libmesh/id_types.h"
29 #include "libmesh/libmesh.h"
32 #include "libmesh/dense_vector.h"
33 
34 #ifdef LIBMESH_FORWARD_DECLARE_ENUMS
35 namespace libMesh
36 {
37 enum SolverPackage : int;
38 }
39 #else
41 #endif
42 
43 // C++ includes
44 #include <cstddef>
45 #include <set>
46 #include <vector>
47 #include <memory>
48 
49 namespace libMesh
50 {
51 
52 
53 // forward declarations
54 template <typename T> class NumericVector;
55 template <typename T> class DenseVector;
56 template <typename T> class DenseSubVector;
57 template <typename T> class SparseMatrix;
58 template <typename T> class ShellMatrix;
59 
74 template <typename T>
75 class NumericVector : public ReferenceCountedObject<NumericVector<T>>,
76  public ParallelObject
77 {
78 public:
79 
83  explicit
84  NumericVector (const Parallel::Communicator & comm_in,
85  const ParallelType ptype = AUTOMATIC);
86 
90  explicit
91  NumericVector (const Parallel::Communicator & comm_in,
92  const numeric_index_type n,
93  const ParallelType ptype = AUTOMATIC);
94 
99  NumericVector (const Parallel::Communicator & comm_in,
100  const numeric_index_type n,
101  const numeric_index_type n_local,
102  const ParallelType ptype = AUTOMATIC);
103 
109  NumericVector (const Parallel::Communicator & comm_in,
110  const numeric_index_type N,
111  const numeric_index_type n_local,
112  const std::vector<numeric_index_type> & ghost,
113  const ParallelType ptype = AUTOMATIC);
114 
126  virtual NumericVector<T> & operator= (const NumericVector<T> & v) = 0;
127 
132  NumericVector (NumericVector &&) = default;
133  NumericVector (const NumericVector &) = default;
134  NumericVector & operator= (NumericVector &&) = default;
135 
141  static std::unique_ptr<NumericVector<T>>
142  build(const Parallel::Communicator & comm,
143  const SolverPackage solver_package = libMesh::default_solver_package());
144 
149  virtual bool initialized() const { return _is_initialized; }
150 
154  ParallelType type() const { return _type; }
155 
159  ParallelType & type() { return _type; }
160 
165  virtual bool closed() const { return _is_closed; }
166 
171  virtual void close () = 0;
172 
176  virtual void clear ();
177 
182  virtual void zero () = 0;
183 
190  virtual std::unique_ptr<NumericVector<T>> zero_clone () const = 0;
191 
197  virtual std::unique_ptr<NumericVector<T>> clone () const = 0;
198 
209  virtual void init (const numeric_index_type n,
210  const numeric_index_type n_local,
211  const bool fast = false,
212  const ParallelType ptype = AUTOMATIC) = 0;
213 
217  virtual void init (const numeric_index_type n,
218  const bool fast = false,
219  const ParallelType ptype = AUTOMATIC) = 0;
220 
225  virtual void init (const numeric_index_type n,
226  const numeric_index_type n_local,
227  const std::vector<numeric_index_type> & ghost,
228  const bool fast = false,
229  const ParallelType ptype = AUTOMATIC) = 0;
230 
235  virtual void init (const NumericVector<T> & other,
236  const bool fast = false) = 0;
237 
243  virtual NumericVector<T> & operator= (const T s) = 0;
244 
250  virtual NumericVector<T> & operator= (const std::vector<T> & v) = 0;
251 
256  virtual Real min () const = 0;
257 
262  virtual Real max () const = 0;
263 
267  virtual T sum() const = 0;
268 
273  virtual Real l1_norm () const = 0;
274 
279  virtual Real l2_norm () const = 0;
280 
285  virtual Real linfty_norm () const = 0;
286 
293  virtual Real subset_l1_norm (const std::set<numeric_index_type> & indices) const;
294 
302  virtual Real subset_l2_norm (const std::set<numeric_index_type> & indices) const;
303 
310  virtual Real subset_linfty_norm (const std::set<numeric_index_type> & indices) const;
311 
315  virtual numeric_index_type size () const = 0;
316 
320  virtual numeric_index_type local_size() const = 0;
321 
328  virtual numeric_index_type first_local_index() const = 0;
329 
336  virtual numeric_index_type last_local_index() const = 0;
337 
341  virtual T operator() (const numeric_index_type i) const = 0;
342 
346  virtual T el(const numeric_index_type i) const { return (*this)(i); }
347 
354  virtual void get(const std::vector<numeric_index_type> & index,
355  T * values) const;
356 
363  void get(const std::vector<numeric_index_type> & index,
364  std::vector<T> & values) const;
365 
373  virtual NumericVector<T> & operator += (const NumericVector<T> & v) = 0;
374 
382  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) = 0;
383 
391  NumericVector<T> & operator *= (const T a) { this->scale(a); return *this; }
392 
400  NumericVector<T> & operator /= (const T a) { this->scale(1./a); return *this; }
401 
408  virtual NumericVector<T> & operator /= (const NumericVector<T> & /*v*/) = 0;
409 
414  virtual void reciprocal() = 0;
415 
419  virtual void conjugate() = 0;
420 
424  virtual void set (const numeric_index_type i, const T value) = 0;
425 
429  virtual void add (const numeric_index_type i, const T value) = 0;
430 
435  virtual void add (const T s) = 0;
436 
442  virtual void add (const NumericVector<T> & v) = 0;
443 
449  virtual void add (const T a, const NumericVector<T> & v) = 0;
450 
457  virtual void add_vector (const T * v,
458  const std::vector<numeric_index_type> & dof_indices);
459 
465  void add_vector (const std::vector<T> & v,
466  const std::vector<numeric_index_type> & dof_indices);
467 
473  virtual void add_vector (const NumericVector<T> & v,
474  const std::vector<numeric_index_type> & dof_indices);
475 
481  void add_vector (const DenseVector<T> & v,
482  const std::vector<numeric_index_type> & dof_indices);
483 
489  virtual void add_vector (const NumericVector<T> & v,
490  const SparseMatrix<T> & A) = 0;
491 
497  void add_vector (const NumericVector<T> & v,
498  const ShellMatrix<T> & A);
499 
505  virtual void add_vector_transpose (const NumericVector<T> & v,
506  const SparseMatrix<T> & A) = 0;
507 
511  virtual void insert (const T * v,
512  const std::vector<numeric_index_type> & dof_indices);
513 
517  void insert (const std::vector<T> & v,
518  const std::vector<numeric_index_type> & dof_indices);
519 
523  virtual void insert (const NumericVector<T> & v,
524  const std::vector<numeric_index_type> & dof_indices);
525 
529  void insert (const DenseVector<T> & v,
530  const std::vector<numeric_index_type> & dof_indices);
531 
535  void insert (const DenseSubVector<T> & v,
536  const std::vector<numeric_index_type> & dof_indices);
537 
541  virtual void scale (const T factor) = 0;
542 
546  virtual void abs() = 0;
547 
554  virtual T dot(const NumericVector<T> & v) const = 0;
555 
560  virtual void localize (std::vector<T> & v_local) const = 0;
561 
566  virtual void localize (NumericVector<T> & v_local) const = 0;
567 
572  virtual void localize (NumericVector<T> & v_local,
573  const std::vector<numeric_index_type> & send_list) const = 0;
574 
604  virtual void localize (std::vector<T> & v_local,
605  const std::vector<numeric_index_type> & indices) const = 0;
606 
611  virtual void localize (const numeric_index_type first_local_idx,
612  const numeric_index_type last_local_idx,
613  const std::vector<numeric_index_type> & send_list) = 0;
614 
620  virtual void localize_to_one (std::vector<T> & v_local,
621  const processor_id_type proc_id=0) const = 0;
622 
628  virtual int compare (const NumericVector<T> & other_vector,
629  const Real threshold = TOLERANCE) const;
630 
636  virtual int local_relative_compare (const NumericVector<T> & other_vector,
637  const Real threshold = TOLERANCE) const;
638 
644  virtual int global_relative_compare (const NumericVector<T> & other_vector,
645  const Real threshold = TOLERANCE) const;
646 
652  virtual void pointwise_mult (const NumericVector<T> & vec1,
653  const NumericVector<T> & vec2) = 0;
654 
659  virtual void print(std::ostream & os=libMesh::out) const;
660 
665  virtual void print_global(std::ostream & os=libMesh::out) const;
666 
670  friend std::ostream & operator << (std::ostream & os, const NumericVector<T> & v)
671  {
672  v.print_global(os);
673  return os;
674  }
675 
681  virtual void print_matlab(const std::string & /*name*/ = "") const
682  {
683  libmesh_not_implemented();
684  }
685 
693  const std::vector<numeric_index_type> &) const
694  {
695  libmesh_not_implemented();
696  }
697 
703  virtual void swap (NumericVector<T> & v);
704 
705 protected:
706 
713 
718 
723 };
724 
725 
726 /*----------------------- Inline functions ----------------------------------*/
727 
728 
729 
730 template <typename T>
731 inline
733  const ParallelType ptype) :
734  ParallelObject(comm_in),
735  _is_closed(false),
736  _is_initialized(false),
737  _type(ptype)
738 {
739 }
740 
741 
742 
743 template <typename T>
744 inline
746  const numeric_index_type /*n*/,
747  const ParallelType ptype) :
748  ParallelObject(comm_in),
749  _is_closed(false),
750  _is_initialized(false),
751  _type(ptype)
752 {
753  libmesh_not_implemented(); // Abstract base class!
754  // init(n, n, false, ptype);
755 }
756 
757 
758 
759 template <typename T>
760 inline
762  const numeric_index_type /*n*/,
763  const numeric_index_type /*n_local*/,
764  const ParallelType ptype) :
765  ParallelObject(comm_in),
766  _is_closed(false),
767  _is_initialized(false),
768  _type(ptype)
769 {
770  libmesh_not_implemented(); // Abstract base class!
771  // init(n, n_local, false, ptype);
772 }
773 
774 
775 
776 template <typename T>
777 inline
779  const numeric_index_type /*n*/,
780  const numeric_index_type /*n_local*/,
781  const std::vector<numeric_index_type> & /*ghost*/,
782  const ParallelType ptype) :
783  ParallelObject(comm_in),
784  _is_closed(false),
785  _is_initialized(false),
786  _type(ptype)
787 {
788  libmesh_not_implemented(); // Abstract base class!
789  // init(n, n_local, ghost, false, ptype);
790 }
791 
792 
793 
794 template <typename T>
795 inline
797 {
798  _is_closed = false;
799  _is_initialized = false;
800 }
801 
802 
803 
804 template <typename T>
805 inline
806 void NumericVector<T>::get(const std::vector<numeric_index_type> & index,
807  T * values) const
808 {
809  const std::size_t num = index.size();
810  for (std::size_t i=0; i<num; i++)
811  {
812  values[i] = (*this)(index[i]);
813  }
814 }
815 
816 
817 
818 template <typename T>
819 inline
820 void NumericVector<T>::get(const std::vector<numeric_index_type> & index,
821  std::vector<T> & values) const
822 {
823  const std::size_t num = index.size();
824  values.resize(num);
825  if (!num)
826  return;
827 
828  this->get(index, values.data());
829 }
830 
831 
832 
833 template <typename T>
834 inline
835 void NumericVector<T>::add_vector(const std::vector<T> & v,
836  const std::vector<numeric_index_type> & dof_indices)
837 {
838  libmesh_assert(v.size() == dof_indices.size());
839  if (!v.empty())
840  this->add_vector(v.data(), dof_indices);
841 }
842 
843 
844 
845 template <typename T>
846 inline
848  const std::vector<numeric_index_type> & dof_indices)
849 {
850  libmesh_assert(v.size() == dof_indices.size());
851  if (!v.empty())
852  this->add_vector(&v(0), dof_indices);
853 }
854 
855 
856 
857 template <typename T>
858 inline
859 void NumericVector<T>::insert(const std::vector<T> & v,
860  const std::vector<numeric_index_type> & dof_indices)
861 {
862  libmesh_assert(v.size() == dof_indices.size());
863  if (!v.empty())
864  this->insert(v.data(), dof_indices);
865 }
866 
867 
868 
869 template <typename T>
870 inline
872  const std::vector<numeric_index_type> & dof_indices)
873 {
874  libmesh_assert(v.size() == dof_indices.size());
875  if (!v.empty())
876  this->insert(&v(0), dof_indices);
877 }
878 
879 
880 
881 template <typename T>
882 inline
884  const std::vector<numeric_index_type> & dof_indices)
885 {
886  libmesh_assert(v.size() == dof_indices.size());
887  if (!v.empty())
888  this->insert(&v(0), dof_indices);
889 }
890 
891 
892 
893 // Full specialization of the print() member for complex
894 // variables. This must precede the non-specialized
895 // version, at least according to icc v7.1
896 template <>
897 inline
898 void NumericVector<Complex>::print(std::ostream & os) const
899 {
900  libmesh_assert (this->initialized());
901  os << "Size\tglobal = " << this->size()
902  << "\t\tlocal = " << this->local_size() << std::endl;
903 
904  // std::complex<>::operator<<() is defined, but use this form
905  os << "#\tReal part\t\tImaginary part" << std::endl;
906  for (numeric_index_type i=this->first_local_index(); i<this->last_local_index(); i++)
907  os << i << "\t"
908  << (*this)(i).real() << "\t\t"
909  << (*this)(i).imag() << std::endl;
910 }
911 
912 
913 
914 template <typename T>
915 inline
916 void NumericVector<T>::print(std::ostream & os) const
917 {
918  libmesh_assert (this->initialized());
919  os << "Size\tglobal = " << this->size()
920  << "\t\tlocal = " << this->local_size() << std::endl;
921 
922  os << "#\tValue" << std::endl;
923  for (numeric_index_type i=this->first_local_index(); i<this->last_local_index(); i++)
924  os << i << "\t" << (*this)(i) << std::endl;
925 }
926 
927 
928 
929 template <>
930 inline
931 void NumericVector<Complex>::print_global(std::ostream & os) const
932 {
933  libmesh_assert (this->initialized());
934 
935  std::vector<Complex> v(this->size());
936  this->localize(v);
937 
938  // Right now we only want one copy of the output
939  if (this->processor_id())
940  return;
941 
942  os << "Size\tglobal = " << this->size() << std::endl;
943  os << "#\tReal part\t\tImaginary part" << std::endl;
944  for (numeric_index_type i=0; i!=v.size(); i++)
945  os << i << "\t"
946  << v[i].real() << "\t\t"
947  << v[i].imag() << std::endl;
948 }
949 
950 
951 template <typename T>
952 inline
953 void NumericVector<T>::print_global(std::ostream & os) const
954 {
955  libmesh_assert (this->initialized());
956 
957  std::vector<T> v(this->size());
958  this->localize(v);
959 
960  // Right now we only want one copy of the output
961  if (this->processor_id())
962  return;
963 
964  os << "Size\tglobal = " << this->size() << std::endl;
965  os << "#\tValue" << std::endl;
966  for (numeric_index_type i=0; i!=v.size(); i++)
967  os << i << "\t" << v[i] << std::endl;
968 }
969 
970 
971 
972 template <typename T>
973 inline
975 {
976  std::swap(_is_closed, v._is_closed);
977  std::swap(_is_initialized, v._is_initialized);
978  std::swap(_type, v._type);
979 }
980 
981 
982 } // namespace libMesh
983 
984 
985 #endif // LIBMESH_NUMERIC_VECTOR_H
virtual unsigned int size() const override
Definition: dense_vector.h:92
virtual void create_subvector(NumericVector< T > &, const std::vector< numeric_index_type > &) const
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices)
virtual Real subset_l2_norm(const std::set< numeric_index_type > &indices) const
virtual void print_global(std::ostream &os=libMesh::out) const
virtual T el(const numeric_index_type i) const
virtual T sum() const =0
virtual void add_vector_transpose(const NumericVector< T > &v, const SparseMatrix< T > &A)=0
virtual bool initialized() const
virtual numeric_index_type size() const =0
virtual std::unique_ptr< NumericVector< T > > zero_clone() const =0
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
virtual NumericVector< T > & operator=(const NumericVector< T > &v)=0
uint8_t processor_id_type
Definition: id_types.h:99
ParallelType & type()
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
const Parallel::Communicator & comm() const
virtual std::unique_ptr< NumericVector< T > > clone() const =0
static const Real TOLERANCE
virtual Real subset_linfty_norm(const std::set< numeric_index_type > &indices) const
virtual T dot(const NumericVector< T > &v) const =0
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
virtual Real subset_l1_norm(const std::set< numeric_index_type > &indices) const
virtual void print_matlab(const std::string &="") const
virtual bool empty() const override
Definition: dense_vector.h:97
virtual void zero()=0
SolverPackage default_solver_package()
Definition: libmesh.C:971
virtual void scale(const T factor)=0
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const =0
virtual Real l2_norm() const =0
virtual NumericVector< T > & operator+=(const NumericVector< T > &v)=0
dof_id_type numeric_index_type
Definition: id_types.h:92
NumericVector(const Parallel::Communicator &comm_in, const ParallelType ptype=AUTOMATIC)
NumericVector< T > & operator/=(const T a)
NumericVector< T > & operator*=(const T a)
virtual int global_relative_compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
virtual void print(std::ostream &os=libMesh::out) const
virtual unsigned int size() const override
virtual Real max() const =0
virtual Real min() const =0
An object whose state is distributed along a set of processors.
virtual Real l1_norm() const =0
virtual void close()=0
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2)=0
virtual bool closed() const
virtual void abs()=0
virtual numeric_index_type first_local_index() const =0
ParallelType type() const
virtual int local_relative_compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual NumericVector< T > & operator-=(const NumericVector< T > &v)=0
virtual void swap(NumericVector< T > &v)
virtual int compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
static PetscErrorCode Mat * A
virtual bool empty() const override
void swap(Iterator &lhs, Iterator &rhs)
virtual numeric_index_type local_size() const =0
static const bool value
Definition: xdr_io.C:109
bool initialized()
Definition: libmesh.C:258
virtual void reciprocal()=0
virtual void add(const numeric_index_type i, const T value)=0
virtual T operator()(const numeric_index_type i) const =0
OStreamProxy out(std::cout)
virtual numeric_index_type last_local_index() const =0
virtual Real linfty_norm() const =0
const Elem & get(const ElemType type_in)
virtual void conjugate()=0
virtual void localize(std::vector< T > &v_local) const =0