include/deal.II/lac/petsc_parallel_vector.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: petsc_parallel_vector.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2004, 2005, 2006, 2007, 2009, 2010, 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_parallel_vector_h
00013 #define __deal2__petsc_parallel_vector_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 #  include <deal.II/lac/vector.h>
00023 #  include <deal.II/lac/petsc_vector_base.h>
00024 #  include <deal.II/base/index_set.h>
00025 
00026 DEAL_II_NAMESPACE_OPEN
00027 
00028 
00029 
00030 // forward declaration
00031 template <typename> class Vector;
00032 class IndexSet;
00033 
00034 
00038 namespace PETScWrappers
00039 {
00047   namespace MPI
00048   {
00049 
00153     class Vector : public VectorBase
00154     {
00155       public:
00160         Vector ();
00161 
00186         explicit Vector (const MPI_Comm     &communicator,
00187                          const unsigned int  n,
00188                          const unsigned int  local_size);
00189 
00190 
00206         template <typename Number>
00207         explicit Vector (const MPI_Comm         &communicator,
00208                          const ::Vector<Number> &v,
00209                          const unsigned int      local_size);
00210 
00211 
00226         explicit Vector (const MPI_Comm     &communicator,
00227                          const VectorBase   &v,
00228                          const unsigned int  local_size);
00229 
00230 
00249         explicit Vector (const MPI_Comm     &communicator,
00250                          const IndexSet &  local,
00251                          const IndexSet & ghost = IndexSet(0));
00252 
00253 
00260         Vector & operator = (const Vector &v);
00261 
00262 
00298         Vector & operator = (const PETScWrappers::Vector &v);
00299 
00309         Vector & operator = (const PetscScalar s);
00310 
00326         template <typename number>
00327         Vector & operator = (const ::Vector<number> &v);
00328 
00355         void reinit (const MPI_Comm     &communicator,
00356                      const unsigned int  N,
00357                      const unsigned int  local_size,
00358                      const bool          fast = false);
00359 
00374         void reinit (const Vector &v,
00375                      const bool    fast = false);
00376 
00382         void reinit (const MPI_Comm     &communicator,
00383                      const IndexSet &  local,
00384                      const IndexSet & ghost = IndexSet(0));
00385 
00386 
00392         const MPI_Comm & get_mpi_communicator () const;
00393 
00394       protected:
00404         virtual void create_vector (const unsigned int n,
00405                                     const unsigned int local_size);
00406 
00407 
00408 
00417         virtual void create_vector (const unsigned int  n,
00418                                     const unsigned int  local_size,
00419                                     const IndexSet & ghostnodes);
00420 
00421 
00422       private:
00427         MPI_Comm communicator;
00428     };
00429 
00430 
00431 // ------------------ template and inline functions -------------
00432 
00433 
00442     inline
00443     void swap (Vector &u, Vector &v)
00444     {
00445       u.swap (v);
00446     }
00447 
00448 
00449 #ifndef DOXYGEN
00450 
00451     template <typename number>
00452     Vector::Vector (const MPI_Comm         &communicator,
00453                     const ::Vector<number> &v,
00454                     const unsigned int      local_size)
00455                     :
00456                     communicator (communicator)
00457     {
00458       Vector::create_vector (v.size(), local_size);
00459 
00460       *this = v;
00461     }
00462 
00463 
00464 
00465     inline
00466     Vector &
00467     Vector::operator = (const PetscScalar s)
00468     {
00469       VectorBase::operator = (s);
00470 
00471       return *this;
00472     }
00473 
00474 
00475 
00476     inline
00477     Vector &
00478     Vector::operator = (const Vector &v)
00479     {
00480                                        // if the vectors have different sizes,
00481                                        // then first resize the present one
00482       if (size() != v.size())
00483         reinit (v.communicator, v.size(), v.local_size(), true);
00484 
00485       const int ierr = VecCopy (v.vector, vector);
00486       AssertThrow (ierr == 0, ExcPETScError(ierr));
00487 
00488       return *this;
00489     }
00490 
00491 
00492 
00493     template <typename number>
00494     inline
00495     Vector &
00496     Vector::operator = (const ::Vector<number> &v)
00497     {
00498       Assert (size() == v.size(),
00499               ExcDimensionMismatch (size(), v.size()));
00500 
00501                                        // the following isn't necessarily fast,
00502                                        // but this is due to the fact that PETSc
00503                                        // doesn't offer an inlined access
00504                                        // operator.
00505                                        //
00506                                        // if someone wants to contribute some
00507                                        // code: to make this code faster, one
00508                                        // could either first convert all values
00509                                        // to PetscScalar, and then set them all
00510                                        // at once using VecSetValues. This has
00511                                        // the drawback that it could take quite
00512                                        // some memory, if the vector is large,
00513                                        // and it would in addition allocate
00514                                        // memory on the heap, which is
00515                                        // expensive. an alternative would be to
00516                                        // split the vector into chunks of, say,
00517                                        // 128 elements, convert a chunk at a
00518                                        // time and set it in the output vector
00519                                        // using VecSetValues. since 128 elements
00520                                        // is small enough, this could easily be
00521                                        // allocated on the stack (as a local
00522                                        // variable) which would make the whole
00523                                        // thing much more efficient.
00524                                        //
00525                                        // a second way to make things faster is
00526                                        // for the special case that
00527                                        // number==PetscScalar. we could then
00528                                        // declare a specialization of this
00529                                        // template, and omit the conversion. the
00530                                        // problem with this is that the best we
00531                                        // can do is to use VecSetValues, but
00532                                        // this isn't very efficient either: it
00533                                        // wants to see an array of indices,
00534                                        // which in this case a) again takes up a
00535                                        // whole lot of memory on the heap, and
00536                                        // b) is totally dumb since its content
00537                                        // would simply be the sequence
00538                                        // 0,1,2,3,...,n. the best of all worlds
00539                                        // would probably be a function in Petsc
00540                                        // that would take a pointer to an array
00541                                        // of PetscScalar values and simply copy
00542                                        // n elements verbatim into the vector...
00543       for (unsigned int i=0; i<v.size(); ++i)
00544         (*this)(i) = v(i);
00545 
00546       compress ();
00547 
00548       return *this;
00549     }
00550 
00551 
00552 
00553     inline
00554     const MPI_Comm &
00555     Vector::get_mpi_communicator () const
00556     {
00557       return communicator;
00558     }
00559 
00560 #endif // DOXYGEN
00561   }
00562 }
00563 
00566 DEAL_II_NAMESPACE_CLOSE
00567 
00568 #endif // DEAL_II_USE_PETSC
00569 
00570 /*----------------------------   petsc_parallel_vector.h     ---------------------------*/
00571 
00572 #endif
00573 /*----------------------------   petsc_parallel_vector.h     ---------------------------*/
00574 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

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