numeric_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 #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"
28 #include "libmesh/id_types.h"
30 #include "libmesh/libmesh.h"
32 
33 // C++ includes
34 #include <cstddef>
35 #include <set>
36 #include <vector>
37 
38 namespace libMesh
39 {
40 
41 
42 // forward declarations
43 template <typename T> class NumericVector;
44 template <typename T> class DenseVector;
45 template <typename T> class DenseSubVector;
46 template <typename T> class SparseMatrix;
47 template <typename T> class ShellMatrix;
48 
49 
57 template <typename T>
58 class NumericVector : public ReferenceCountedObject<NumericVector<T> >,
59  public ParallelObject
60 {
61 public:
62 
66  explicit
67  NumericVector (const Parallel::Communicator &comm_in,
68  const ParallelType ptype = AUTOMATIC);
69 
73  explicit
74  NumericVector (const Parallel::Communicator &comm_in,
75  const numeric_index_type n,
76  const ParallelType ptype = AUTOMATIC);
77 
82  NumericVector (const Parallel::Communicator &comm_in,
83  const numeric_index_type n,
85  const ParallelType ptype = AUTOMATIC);
86 
92  NumericVector (const Parallel::Communicator &comm_in,
93  const numeric_index_type N,
95  const std::vector<numeric_index_type>& ghost,
96  const ParallelType ptype = AUTOMATIC);
97 
98 public:
99 
104  virtual ~NumericVector ();
105 
111  static AutoPtr<NumericVector<T> >
112  build(const Parallel::Communicator &comm,
113  const SolverPackage solver_package = libMesh::default_solver_package());
114 
115 #ifndef LIBMESH_DISABLE_COMMWORLD
116 
121  static AutoPtr<NumericVector<T> >
122  build(const SolverPackage solver_package = libMesh::default_solver_package());
123 #endif
124 
129  virtual bool initialized() const { return _is_initialized; }
130 
134  ParallelType type() const { return _type; }
135 
139  ParallelType & type() { return _type; }
140 
145  virtual bool closed() const { return _is_closed; }
146 
150  virtual void close () = 0;
151 
155  virtual void clear ();
156 
161  virtual void zero () = 0;
162 
169  virtual AutoPtr<NumericVector<T> > zero_clone () const = 0;
170 
175  virtual AutoPtr<NumericVector<T> > clone () const = 0;
176 
190  virtual void init (const numeric_index_type,
191  const numeric_index_type,
192  const bool = false,
193  const ParallelType = AUTOMATIC) = 0;
194 
198  virtual void init (const numeric_index_type,
199  const bool = false,
200  const ParallelType = AUTOMATIC) = 0;
201 
206  virtual void init (const numeric_index_type /*N*/,
207  const numeric_index_type /*n_local*/,
208  const std::vector<numeric_index_type>& /*ghost*/,
209  const bool /*fast*/ = false,
210  const ParallelType = AUTOMATIC) = 0;
211 
216  virtual void init (const NumericVector<T>& other,
217  const bool fast = false) = 0;
218 
222  virtual NumericVector<T> & operator= (const T s) = 0;
223 
227  virtual NumericVector<T> & operator= (const NumericVector<T> &V) = 0;
228 
232  virtual NumericVector<T> & operator= (const std::vector<T> &v) = 0;
233 
239  virtual Real min () const = 0;
240 
246  virtual Real max () const = 0;
247 
251  virtual T sum() const = 0;
252 
257  virtual Real l1_norm () const = 0;
258 
264  virtual Real l2_norm () const = 0;
265 
271  virtual Real linfty_norm () const = 0;
272 
281  virtual Real subset_l1_norm (const std::set<numeric_index_type> & indices) const;
282 
292  virtual Real subset_l2_norm (const std::set<numeric_index_type> & indices) const;
293 
302  virtual Real subset_linfty_norm (const std::set<numeric_index_type> & indices) const;
303 
311  virtual numeric_index_type size () const = 0;
312 
318  virtual numeric_index_type local_size() const = 0;
319 
325  virtual numeric_index_type first_local_index() const = 0;
326 
332  virtual numeric_index_type last_local_index() const = 0;
333 
337  virtual T operator() (const numeric_index_type i) const = 0;
338 
342  virtual T el(const numeric_index_type i) const { return (*this)(i); }
343 
350  virtual void get(const std::vector<numeric_index_type>& index, std::vector<T>& values) const;
351 
356  virtual NumericVector<T> & operator += (const NumericVector<T> &V) = 0;
357 
362  virtual NumericVector<T> & operator -= (const NumericVector<T> &V) = 0;
363 
368  NumericVector<T> & operator *= (const T a) { this->scale(a); return *this; }
369 
374  NumericVector<T> & operator /= (const T a) { this->scale(1./a); return *this; }
375 
379  virtual NumericVector<T> & operator /= (NumericVector<T> & /*v*/) = 0;
380 
384  virtual void reciprocal() = 0;
385 
390  virtual void conjugate() = 0;
391 
395  virtual void set (const numeric_index_type i, const T value) = 0;
396 
400  virtual void add (const numeric_index_type i, const T value) = 0;
401 
407  virtual void add (const T s) = 0;
408 
414  virtual void add (const NumericVector<T>& V) = 0;
415 
421  virtual void add (const T a, const NumericVector<T>& v) = 0;
422 
428  virtual void add_vector (const std::vector<T>& v,
429  const std::vector<numeric_index_type>& dof_indices) = 0;
430 
437  virtual void add_vector (const NumericVector<T>& V,
438  const std::vector<numeric_index_type>& dof_indices) = 0;
439 
444  virtual void add_vector (const NumericVector<T>&,
445  const SparseMatrix<T>&) = 0;
446 
451  void add_vector (const NumericVector<T>& v,
452  const ShellMatrix<T>& a);
453 
460  virtual void add_vector (const DenseVector<T>& V,
461  const std::vector<numeric_index_type>& dof_indices) = 0;
462 
467  virtual void add_vector_transpose (const NumericVector<T>&,
468  const SparseMatrix<T>&) = 0;
469 
474  virtual void insert (const std::vector<T>& v,
475  const std::vector<numeric_index_type>& dof_indices) = 0;
476 
483  virtual void insert (const NumericVector<T>& V,
484  const std::vector<numeric_index_type>& dof_indices) = 0;
485 
492  virtual void insert (const DenseVector<T>& V,
493  const std::vector<numeric_index_type>& dof_indices) = 0;
494 
500  virtual void insert (const DenseSubVector<T>& V,
501  const std::vector<numeric_index_type>& dof_indices) = 0;
502 
507  virtual void scale (const T factor) = 0;
508 
513  virtual void abs() = 0;
514 
518  virtual T dot(const NumericVector<T>&) const = 0;
519 
524  virtual void localize (std::vector<T>& v_local) const = 0;
525 
530  virtual void localize (NumericVector<T>& v_local) const = 0;
531 
537  virtual void localize (NumericVector<T>& v_local,
538  const std::vector<numeric_index_type>& send_list) const = 0;
539 
544  virtual void localize (const numeric_index_type first_local_idx,
545  const numeric_index_type last_local_idx,
546  const std::vector<numeric_index_type>& send_list) = 0;
547 
554  virtual void localize_to_one (std::vector<T>& v_local,
555  const processor_id_type proc_id=0) const = 0;
556 
565  virtual int compare (const NumericVector<T> &other_vector,
566  const Real threshold = TOLERANCE) const;
567 
576  virtual int local_relative_compare (const NumericVector<T> &other_vector,
577  const Real threshold = TOLERANCE) const;
578 
587  virtual int global_relative_compare (const NumericVector<T> &other_vector,
588  const Real threshold = TOLERANCE) const;
589 
594  virtual void pointwise_mult (const NumericVector<T>& vec1,
595  const NumericVector<T>& vec2) = 0;
596 
601  virtual void print(std::ostream& os=libMesh::out) const;
602 
607  virtual void print_global(std::ostream& os=libMesh::out) const;
608 
612  friend std::ostream& operator << (std::ostream& os, const NumericVector<T>& v)
613  {
614  v.print_global(os);
615  return os;
616  }
617 
624  virtual void print_matlab(const std::string name="NULL") const
625  {
626  libMesh::err << "ERROR: Not Implemented in base class yet!" << std::endl;
627  libMesh::err << "ERROR writing MATLAB file " << name << std::endl;
628  libmesh_error();
629  }
630 
638  const std::vector<numeric_index_type>& ) const
639  {
640  libMesh::err << "ERROR: Not Implemented in base class yet!" << std::endl;
641  libmesh_error();
642  }
643 
649  virtual void swap (NumericVector<T> &v);
650 
651 protected:
652 
658 
664 
669 };
670 
671 
672 /*----------------------- Inline functions ----------------------------------*/
673 
674 
675 
676 template <typename T>
677 inline
679  const ParallelType ptype) :
680  ParallelObject(comm_in),
681  _is_closed(false),
682  _is_initialized(false),
683  _type(ptype)
684 {
685 }
686 
687 
688 
689 template <typename T>
690 inline
692  const numeric_index_type /*n*/,
693  const ParallelType ptype) :
694  ParallelObject(comm_in),
695  _is_closed(false),
696  _is_initialized(false),
697  _type(ptype)
698 {
699  libmesh_error(); // Abstract base class!
700  // init(n, n, false, ptype);
701 }
702 
703 
704 
705 template <typename T>
706 inline
708  const numeric_index_type /*n*/,
709  const numeric_index_type /*n_local*/,
710  const ParallelType ptype) :
711  ParallelObject(comm_in),
712  _is_closed(false),
713  _is_initialized(false),
714  _type(ptype)
715 {
716  libmesh_error(); // Abstract base class!
717  // init(n, n_local, false, ptype);
718 }
719 
720 
721 
722 template <typename T>
723 inline
725  const numeric_index_type /*n*/,
726  const numeric_index_type /*n_local*/,
727  const std::vector<numeric_index_type>& /*ghost*/,
728  const ParallelType ptype) :
729  ParallelObject(comm_in),
730  _is_closed(false),
731  _is_initialized(false),
732  _type(ptype)
733 {
734  libmesh_error(); // Abstract base class!
735  // init(n, n_local, ghost, false, ptype);
736 }
737 
738 
739 
740 template <typename T>
741 inline
743 {
744  clear ();
745 }
746 
747 
748 
749 // These should be pure virtual, not bugs waiting to happen - RHS
750 /*
751 template <typename T>
752 inline
753 NumericVector<T> & NumericVector<T>::operator= (const T)
754 {
755  // libmesh_error();
756 
757  return *this;
758 }
759 
760 
761 
762 template <typename T>
763 inline
764 NumericVector<T> & NumericVector<T>::operator= (const NumericVector<T>&)
765 {
766  // libmesh_error();
767 
768  return *this;
769 }
770 
771 
772 
773 template <typename T>
774 inline
775 NumericVector<T> & NumericVector<T>::operator= (const std::vector<T>&)
776 {
777  // libmesh_error();
778 
779  return *this;
780 }
781 */
782 
783 
784 
785 template <typename T>
786 inline
788 {
789  _is_closed = false;
790  _is_initialized = false;
791 }
792 
793 
794 
795 template <typename T>
796 inline
797 void NumericVector<T>::get(const std::vector<numeric_index_type>& index, std::vector<T>& values) const
798 {
799  const std::size_t num = index.size();
800  values.resize(num);
801  for(std::size_t i=0; i<num; i++)
802  {
803  values[i] = (*this)(index[i]);
804  }
805 }
806 
807 
808 
809 // Full specialization of the print() member for complex
810 // variables. This must precede the non-specialized
811 // version, at least according to icc v7.1
812 template <>
813 inline
814 void NumericVector<Complex>::print(std::ostream& os) const
815 {
816  libmesh_assert (this->initialized());
817  os << "Size\tglobal = " << this->size()
818  << "\t\tlocal = " << this->local_size() << std::endl;
819 
820  // std::complex<>::operator<<() is defined, but use this form
821  os << "#\tReal part\t\tImaginary part" << std::endl;
822  for (numeric_index_type i=this->first_local_index(); i<this->last_local_index(); i++)
823  os << i << "\t"
824  << (*this)(i).real() << "\t\t"
825  << (*this)(i).imag() << std::endl;
826 }
827 
828 
829 
830 template <typename T>
831 inline
832 void NumericVector<T>::print(std::ostream& os) const
833 {
834  libmesh_assert (this->initialized());
835  os << "Size\tglobal = " << this->size()
836  << "\t\tlocal = " << this->local_size() << std::endl;
837 
838  os << "#\tValue" << std::endl;
839  for (numeric_index_type i=this->first_local_index(); i<this->last_local_index(); i++)
840  os << i << "\t" << (*this)(i) << std::endl;
841 }
842 
843 
844 
845 template <>
846 inline
847 void NumericVector<Complex>::print_global(std::ostream& os) const
848 {
849  libmesh_assert (this->initialized());
850 
851  std::vector<Complex> v(this->size());
852  this->localize(v);
853 
854  // Right now we only want one copy of the output
855  if (this->processor_id())
856  return;
857 
858  os << "Size\tglobal = " << this->size() << std::endl;
859  os << "#\tReal part\t\tImaginary part" << std::endl;
860  for (numeric_index_type i=0; i!=v.size(); i++)
861  os << i << "\t"
862  << v[i].real() << "\t\t"
863  << v[i].imag() << std::endl;
864 }
865 
866 
867 template <typename T>
868 inline
869 void NumericVector<T>::print_global(std::ostream& os) const
870 {
871  libmesh_assert (this->initialized());
872 
873  std::vector<T> v(this->size());
874  this->localize(v);
875 
876  // Right now we only want one copy of the output
877  if (this->processor_id())
878  return;
879 
880  os << "Size\tglobal = " << this->size() << std::endl;
881  os << "#\tValue" << std::endl;
882  for (numeric_index_type i=0; i!=v.size(); i++)
883  os << i << "\t" << v[i] << std::endl;
884 }
885 
886 
887 
888 template <typename T>
889 inline
891 {
892  std::swap(_is_closed, v._is_closed);
893  std::swap(_is_initialized, v._is_initialized);
894  std::swap(_type, v._type);
895 }
896 
897 
898 } // namespace libMesh
899 
900 
901 #endif // LIBMESH_NUMERIC_VECTOR_H

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

Hosted By:
SourceForge.net Logo