include/deal.II/lac/petsc_vector_base.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: petsc_vector_base.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2004, 2005, 2006, 2007, 2009, 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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                                  // forward declaration
00032 template <typename number> class Vector;
00033 
00034 
00042 namespace PETScWrappers
00043 {
00044                                    // forward declaration
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 // ------------------- inline and template functions --------------
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                                        // as explained in the class
00915                                        // documentation, this is not the copy
00916                                        // operator. so simply pass on to the
00917                                        // "correct" assignment operator
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                                        // as explained in the class
00930                                        // documentation, this is not the copy
00931                                        // operator. so simply pass on to the
00932                                        // "correct" assignment operator
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                                        // we have to do above actions in any
00980                                        // case to be consistent with the MPI
00981                                        // communication model (see the
00982                                        // comments in the documentation of
00983                                        // PETScWrappers::MPI::Vector), but we
00984                                        // can save some work if the addend is
00985                                        // zero
00986       if (value == PetscScalar())
00987         return *this;
00988 
00989                                        // use the PETSc function to add something
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                                        // we have to do above actions in any
01018                                        // case to be consistent with the MPI
01019                                        // communication model (see the
01020                                        // comments in the documentation of
01021                                        // PETScWrappers::MPI::Vector), but we
01022                                        // can save some work if the addend is
01023                                        // zero
01024       if (value == PetscScalar())
01025         return *this;
01026 
01027                                        // use the PETSc function to
01028                                        // add something
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                                        // we have to do above actions in any
01057                                        // case to be consistent with the MPI
01058                                        // communication model (see the
01059                                        // comments in the documentation of
01060                                        // PETScWrappers::MPI::Vector), but we
01061                                        // can save some work if the factor is
01062                                        // one
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                                        // we have to do above actions in any
01097                                        // case to be consistent with the MPI
01098                                        // communication model (see the
01099                                        // comments in the documentation of
01100                                        // PETScWrappers::MPI::Vector), but we
01101                                        // can save some work if the factor is
01102                                        // one
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 /*----------------------------   petsc_vector_base.h     ---------------------------*/
01195 
01196 #endif
01197 /*----------------------------   petsc_vector_base.h     ---------------------------*/
01198 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:39 by doxygen 1.7.3