00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__petsc_vector_base_h
00013 #define __deal2__petsc_vector_base_h
00014
00015
00016 #include <deal.II/base/config.h>
00017
00018 #ifdef DEAL_II_USE_PETSC
00019
00020 # include <deal.II/base/subscriptor.h>
00021 # include <deal.II/lac/exceptions.h>
00022
00023 # include <vector>
00024 # include <utility>
00025
00026 # include <petscvec.h>
00027 # include <deal.II/base/index_set.h>
00028
00029 DEAL_II_NAMESPACE_OPEN
00030
00031
00032 template <typename number> class Vector;
00033
00034
00042 namespace PETScWrappers
00043 {
00044
00045 class VectorBase;
00046
00056 namespace internal
00057 {
00077 class VectorReference
00078 {
00079 private:
00085 VectorReference (const VectorBase &vector,
00086 const unsigned int index);
00087
00088 public:
00109 const VectorReference & operator = (const VectorReference &r) const;
00110
00119 VectorReference & operator = (const VectorReference &r);
00120
00125 const VectorReference & operator = (const PetscScalar &s) const;
00126
00131 const VectorReference & operator += (const PetscScalar &s) const;
00132
00137 const VectorReference & operator -= (const PetscScalar &s) const;
00138
00143 const VectorReference & operator *= (const PetscScalar &s) const;
00144
00149 const VectorReference & operator /= (const PetscScalar &s) const;
00150
00157 operator PetscScalar () const;
00158
00162 DeclException1 (ExcPETScError,
00163 int,
00164 << "An error with error number " << arg1
00165 << " occured while calling a PETSc function");
00169 DeclException3 (ExcAccessToNonlocalElement,
00170 int, int, int,
00171 << "You tried to access element " << arg1
00172 << " of a distributed vector, but only elements "
00173 << arg2 << " through " << arg3
00174 << " are stored locally and can be accessed.");
00178 DeclException2 (ExcWrongMode,
00179 int, int,
00180 << "You tried to do a "
00181 << (arg1 == 1 ?
00182 "'set'" :
00183 (arg1 == 2 ?
00184 "'add'" : "???"))
00185 << " operation but the vector is currently in "
00186 << (arg2 == 1 ?
00187 "'set'" :
00188 (arg2 == 2 ?
00189 "'add'" : "???"))
00190 << " mode. You first have to call 'compress()'.");
00191
00192 private:
00197 const VectorBase &vector;
00198
00203 const unsigned int index;
00204
00210 friend class ::PETScWrappers::VectorBase;
00211 };
00212 }
00242 class VectorBase : public Subscriptor
00243 {
00244 public:
00252 typedef PetscScalar value_type;
00253 typedef PetscReal real_type;
00254 typedef std::size_t size_type;
00255 typedef internal::VectorReference reference;
00256 typedef const internal::VectorReference const_reference;
00257
00263 VectorBase ();
00264
00270 VectorBase (const VectorBase &v);
00271
00279 explicit VectorBase (const Vec & v);
00280
00284 virtual ~VectorBase ();
00285
00299 void compress ();
00300
00321 VectorBase & operator = (const PetscScalar s);
00322
00331 bool operator == (const VectorBase &v) const;
00332
00341 bool operator != (const VectorBase &v) const;
00342
00347 unsigned int size () const;
00348
00362 unsigned int local_size () const;
00363
00380 std::pair<unsigned int, unsigned int>
00381 local_range () const;
00382
00388 bool in_local_range (const unsigned int index) const;
00389
00394 bool has_ghost_elements() const;
00395
00400 reference
00401 operator () (const unsigned int index);
00402
00407 PetscScalar
00408 operator () (const unsigned int index) const;
00409
00416 reference
00417 operator [] (const unsigned int index);
00418
00426 PetscScalar
00427 operator [] (const unsigned int index) const;
00428
00439 void set (const std::vector<unsigned int> &indices,
00440 const std::vector<PetscScalar> &values);
00441
00448 void add (const std::vector<unsigned int> &indices,
00449 const std::vector<PetscScalar> &values);
00450
00458 void add (const std::vector<unsigned int> &indices,
00459 const ::Vector<PetscScalar> &values);
00460
00470 void add (const unsigned int n_elements,
00471 const unsigned int *indices,
00472 const PetscScalar *values);
00473
00479 PetscScalar operator * (const VectorBase &vec) const;
00480
00484 real_type norm_sqr () const;
00485
00490 PetscScalar mean_value () const;
00491
00496 real_type l1_norm () const;
00497
00503 real_type l2_norm () const;
00504
00511 real_type lp_norm (const real_type p) const;
00512
00517 real_type linfty_norm () const;
00518
00525 real_type normalize () const;
00526
00531 real_type min () const;
00532
00537 real_type max () const;
00538
00539
00545 VectorBase & abs ();
00546
00550 VectorBase & conjugate ();
00551
00558 VectorBase & mult ();
00559
00565 VectorBase & mult (const VectorBase &v);
00566
00575 bool all_zero () const;
00576
00585 bool is_non_negative () const;
00586
00591 VectorBase & operator *= (const PetscScalar factor);
00592
00597 VectorBase & operator /= (const PetscScalar factor);
00598
00603 VectorBase & operator += (const VectorBase &V);
00604
00609 VectorBase & operator -= (const VectorBase &V);
00610
00616 void add (const PetscScalar s);
00617
00622 void add (const VectorBase &V);
00623
00628 void add (const PetscScalar a, const VectorBase &V);
00629
00634 void add (const PetscScalar a, const VectorBase &V,
00635 const PetscScalar b, const VectorBase &W);
00636
00642 void sadd (const PetscScalar s,
00643 const VectorBase &V);
00644
00649 void sadd (const PetscScalar s,
00650 const PetscScalar a,
00651 const VectorBase &V);
00652
00656 void sadd (const PetscScalar s,
00657 const PetscScalar a,
00658 const VectorBase &V,
00659 const PetscScalar b,
00660 const VectorBase &W);
00661
00666 void sadd (const PetscScalar s,
00667 const PetscScalar a,
00668 const VectorBase &V,
00669 const PetscScalar b,
00670 const VectorBase &W,
00671 const PetscScalar c,
00672 const VectorBase &X);
00673
00683 void scale (const VectorBase &scaling_factors);
00684
00688 void equ (const PetscScalar a, const VectorBase &V);
00689
00693 void equ (const PetscScalar a, const VectorBase &V,
00694 const PetscScalar b, const VectorBase &W);
00695
00712 void ratio (const VectorBase &a,
00713 const VectorBase &b);
00714
00721 void update_ghost_values() const;
00722
00737 void print (std::ostream &out,
00738 const unsigned int precision = 3,
00739 const bool scientific = true,
00740 const bool across = true) const;
00741
00765 void swap (VectorBase &v);
00766
00778 operator const Vec & () const;
00779
00785 std::size_t memory_consumption () const;
00786
00787 protected:
00793 Vec vector;
00794
00802 bool ghosted;
00803
00811 IndexSet ghost_indices;
00812
00834 struct LastAction
00835 {
00836 enum Values { none, insert, add };
00837 };
00838
00839
00848 mutable LastAction::Values last_action;
00849
00853 friend class internal::VectorReference;
00854
00862 bool attained_ownership;
00863
00872 void do_set_add_operation (const unsigned int n_elements,
00873 const unsigned int *indices,
00874 const PetscScalar *values,
00875 const bool add_values);
00876
00877
00878 };
00879
00880
00881
00882
00883
00892 inline
00893 void swap (VectorBase &u, VectorBase &v)
00894 {
00895 u.swap (v);
00896 }
00897
00898 #ifndef DOXYGEN
00899 namespace internal
00900 {
00901 inline
00902 VectorReference::VectorReference (const VectorBase &vector,
00903 const unsigned int index)
00904 :
00905 vector (vector),
00906 index (index)
00907 {}
00908
00909
00910 inline
00911 const VectorReference &
00912 VectorReference::operator = (const VectorReference &r) const
00913 {
00914
00915
00916
00917
00918 *this = static_cast<PetscScalar> (r);
00919
00920 return *this;
00921 }
00922
00923
00924
00925 inline
00926 VectorReference &
00927 VectorReference::operator = (const VectorReference &r)
00928 {
00929
00930
00931
00932
00933 *this = static_cast<PetscScalar> (r);
00934
00935 return *this;
00936 }
00937
00938
00939
00940 inline
00941 const VectorReference &
00942 VectorReference::operator = (const PetscScalar &value) const
00943 {
00944 Assert ((vector.last_action == VectorBase::LastAction::insert)
00945 ||
00946 (vector.last_action == VectorBase::LastAction::none),
00947 ExcWrongMode (VectorBase::LastAction::insert,
00948 vector.last_action));
00949
00950 #ifdef PETSC_USE_64BIT_INDICES
00951 const PetscInt petsc_i = index;
00952 #else
00953 const signed int petsc_i = index;
00954 #endif
00955
00956 const int ierr
00957 = VecSetValues (vector, 1, &petsc_i, &value, INSERT_VALUES);
00958 AssertThrow (ierr == 0, ExcPETScError(ierr));
00959
00960 vector.last_action = VectorBase::LastAction::insert;
00961
00962 return *this;
00963 }
00964
00965
00966
00967 inline
00968 const VectorReference &
00969 VectorReference::operator += (const PetscScalar &value) const
00970 {
00971 Assert ((vector.last_action == VectorBase::LastAction::add)
00972 ||
00973 (vector.last_action == VectorBase::LastAction::none),
00974 ExcWrongMode (VectorBase::LastAction::add,
00975 vector.last_action));
00976
00977 vector.last_action = VectorBase::LastAction::add;
00978
00979
00980
00981
00982
00983
00984
00985
00986 if (value == PetscScalar())
00987 return *this;
00988
00989
00990 #ifdef PETSC_USE_64BIT_INDICES
00991 const PetscInt petsc_i = index;
00992 #else
00993 const signed int petsc_i = index;
00994 #endif
00995 const int ierr
00996 = VecSetValues (vector, 1, &petsc_i, &value, ADD_VALUES);
00997 AssertThrow (ierr == 0, ExcPETScError(ierr));
00998
00999
01000 return *this;
01001 }
01002
01003
01004
01005 inline
01006 const VectorReference &
01007 VectorReference::operator -= (const PetscScalar &value) const
01008 {
01009 Assert ((vector.last_action == VectorBase::LastAction::add)
01010 ||
01011 (vector.last_action == VectorBase::LastAction::none),
01012 ExcWrongMode (VectorBase::LastAction::add,
01013 vector.last_action));
01014
01015 vector.last_action = VectorBase::LastAction::add;
01016
01017
01018
01019
01020
01021
01022
01023
01024 if (value == PetscScalar())
01025 return *this;
01026
01027
01028
01029 #ifdef PETSC_USE_64BIT_INDICES
01030 const PetscInt petsc_i = index;
01031 #else
01032 const signed int petsc_i = index;
01033 #endif
01034 const PetscScalar subtractand = -value;
01035 const int ierr
01036 = VecSetValues (vector, 1, &petsc_i, &subtractand, ADD_VALUES);
01037 AssertThrow (ierr == 0, ExcPETScError(ierr));
01038
01039 return *this;
01040 }
01041
01042
01043
01044 inline
01045 const VectorReference &
01046 VectorReference::operator *= (const PetscScalar &value) const
01047 {
01048 Assert ((vector.last_action == VectorBase::LastAction::insert)
01049 ||
01050 (vector.last_action == VectorBase::LastAction::none),
01051 ExcWrongMode (VectorBase::LastAction::insert,
01052 vector.last_action));
01053
01054 vector.last_action = VectorBase::LastAction::insert;
01055
01056
01057
01058
01059
01060
01061
01062
01063 if (value == 1.)
01064 return *this;
01065
01066 #ifdef PETSC_USE_64BIT_INDICES
01067 const PetscInt petsc_i = index;
01068 #else
01069 const signed int petsc_i = index;
01070 #endif
01071
01072 const PetscScalar new_value
01073 = static_cast<PetscScalar>(*this) * value;
01074
01075 const int ierr
01076 = VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES);
01077 AssertThrow (ierr == 0, ExcPETScError(ierr));
01078
01079 return *this;
01080 }
01081
01082
01083
01084 inline
01085 const VectorReference &
01086 VectorReference::operator /= (const PetscScalar &value) const
01087 {
01088 Assert ((vector.last_action == VectorBase::LastAction::insert)
01089 ||
01090 (vector.last_action == VectorBase::LastAction::none),
01091 ExcWrongMode (VectorBase::LastAction::insert,
01092 vector.last_action));
01093
01094 vector.last_action = VectorBase::LastAction::insert;
01095
01096
01097
01098
01099
01100
01101
01102
01103 if (value == 1.)
01104 return *this;
01105
01106 #ifdef PETSC_USE_64BIT_INDICES
01107 const PetscInt petsc_i = index;
01108 #else
01109 const signed int petsc_i = index;
01110 #endif
01111
01112 const PetscScalar new_value
01113 = static_cast<PetscScalar>(*this) / value;
01114
01115 const int ierr
01116 = VecSetValues (vector, 1, &petsc_i, &new_value, INSERT_VALUES);
01117 AssertThrow (ierr == 0, ExcPETScError(ierr));
01118
01119 return *this;
01120 }
01121 }
01122
01123
01124
01125 inline
01126 bool
01127 VectorBase::in_local_range (const unsigned int index) const
01128 {
01129 #ifdef PETSC_USE_64BIT_INDICES
01130 PetscInt begin, end;
01131 #else
01132 int begin, end;
01133 #endif
01134
01135 const int ierr = VecGetOwnershipRange (static_cast<const Vec &>(vector),
01136 &begin, &end);
01137 AssertThrow (ierr == 0, ExcPETScError(ierr));
01138
01139 return ((index >= static_cast<unsigned int>(begin)) &&
01140 (index < static_cast<unsigned int>(end)));
01141 }
01142
01143 inline
01144 bool
01145 VectorBase::has_ghost_elements() const
01146 {
01147 return ghosted;
01148 }
01149
01150
01151
01152 inline
01153 internal::VectorReference
01154 VectorBase::operator () (const unsigned int index)
01155 {
01156 return internal::VectorReference (*this, index);
01157 }
01158
01159
01160
01161 inline
01162 PetscScalar
01163 VectorBase::operator () (const unsigned int index) const
01164 {
01165 return static_cast<PetscScalar>(internal::VectorReference (*this, index));
01166 }
01167
01168
01169
01170 inline
01171 internal::VectorReference
01172 VectorBase::operator [] (const unsigned int index)
01173 {
01174 return operator()(index);
01175 }
01176
01177
01178
01179 inline
01180 PetscScalar
01181 VectorBase::operator [] (const unsigned int index) const
01182 {
01183 return operator()(index);
01184 }
01185
01186
01187 #endif // DOXYGEN
01188 }
01189
01190 DEAL_II_NAMESPACE_CLOSE
01191
01192 #endif // DEAL_II_USE_PETSC
01193
01194
01195
01196 #endif
01197
01198