include/deal.II/lac/petsc_vector.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: petsc_vector.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2004, 2005, 2006, 2007, 2009, 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_h
00013 #define __deal2__petsc_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/petsc_vector_base.h>
00023 #  include <deal.II/lac/petsc_parallel_vector.h>
00024 #  include <deal.II/lac/vector.h>
00025 
00026 DEAL_II_NAMESPACE_OPEN
00027 
00028 
00029 namespace PETScWrappers
00030 {
00047   class Vector : public VectorBase
00048   {
00049     public:
00054       Vector ();
00055 
00070       explicit Vector (const unsigned int n);
00071 
00078       template <typename Number>
00079       explicit Vector (const ::Vector<Number> &v);
00080 
00091       explicit Vector (const Vec & v);
00092 
00097       Vector (const Vector &v);
00098 
00113       explicit Vector (const MPI::Vector &v);
00114 
00119       Vector & operator = (const Vector &v);
00120 
00134       Vector & operator = (const MPI::Vector &v);
00135 
00155       Vector & operator = (const PetscScalar s);
00156 
00163       template <typename number>
00164       Vector & operator = (const ::Vector<number> &v);
00165 
00183       void reinit (const unsigned int N,
00184                    const bool         fast = false);
00185 
00196       void reinit (const Vector &v,
00197                    const bool    fast = false);
00198 
00199     protected:
00207       void create_vector (const unsigned int n);
00208   };
00209 
00212 // ------------------ template and inline functions -------------
00213 
00214 
00223   inline
00224   void swap (Vector &u, Vector &v)
00225   {
00226     u.swap (v);
00227   }
00228 
00229 
00230 #ifndef DOXYGEN
00231 
00232   template <typename number>
00233   Vector::Vector (const ::Vector<number> &v)
00234   {
00235     Vector::create_vector (v.size());
00236 
00237     *this = v;
00238   }
00239 
00240 
00241 
00242   inline
00243   Vector::Vector (const Vec & v)
00244                   :
00245                   VectorBase(v)
00246   {}
00247 
00248 
00249   inline
00250   Vector &
00251   Vector::operator = (const PetscScalar s)
00252   {
00253     VectorBase::operator = (s);
00254 
00255     return *this;
00256   }
00257 
00258 
00259   inline
00260   Vector &
00261   Vector::operator = (const Vector &v)
00262   {
00263                                      // if the vectors have different sizes,
00264                                      // then first resize the present one
00265     if (size() != v.size())
00266       reinit (v.size(), true);
00267 
00268     const int ierr = VecCopy (v.vector, vector);
00269     AssertThrow (ierr == 0, ExcPETScError(ierr));
00270 
00271     return *this;
00272   }
00273 
00274 
00275 
00276   inline
00277   Vector &
00278   Vector::operator = (const MPI::Vector &v)
00279   {
00280     int ierr;
00281     if (attained_ownership)
00282       {
00283                                      // the petsc function we call wants to
00284                                      // generate the vector itself, so destroy
00285                                      // the old one first
00286 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
00287         ierr = VecDestroy (vector);
00288 #else
00289         ierr = VecDestroy (&vector);
00290 #endif
00291         AssertThrow (ierr == 0, ExcPETScError(ierr));
00292       }
00293 
00294     attained_ownership = true;
00295 
00296                                      // then do the gather
00297                                      // operation. <rant>petsc has changed its
00298                                      // interface again, and replaced a single
00299                                      // function call by several calls that
00300                                      // are hard to understand. gets me all
00301                                      // annoyed at their development
00302                                      // model</rant>
00303 #if DEAL_II_PETSC_VERSION_LT(2,2,0)
00304     ierr = VecConvertMPIToSeqAll (static_cast<const Vec &>(v),
00305                                   &vector);
00306     AssertThrow (ierr == 0, ExcPETScError(ierr));
00307 
00308 #else
00309 
00310     VecScatter ctx;
00311 
00312     ierr = VecScatterCreateToAll (static_cast<const Vec &>(v), &ctx, &vector);
00313     AssertThrow (ierr == 0, ExcPETScError(ierr));
00314 
00315 #if ((PETSC_VERSION_MAJOR == 2) && \
00316      ((PETSC_VERSION_MINOR < 3) || \
00317       ((PETSC_VERSION_MINOR == 3) &&            \
00318        (PETSC_VERSION_SUBMINOR < 3))))
00319     ierr = VecScatterBegin (static_cast<const Vec &>(v), vector,
00320                             INSERT_VALUES, SCATTER_FORWARD, ctx);
00321     AssertThrow (ierr == 0, ExcPETScError(ierr));
00322 
00323     ierr = VecScatterEnd (static_cast<const Vec &>(v), vector,
00324                           INSERT_VALUES, SCATTER_FORWARD, ctx);
00325 
00326 #else
00327 
00328     ierr = VecScatterBegin (ctx,static_cast<const Vec &>(v), vector,
00329                             INSERT_VALUES, SCATTER_FORWARD);
00330     AssertThrow (ierr == 0, ExcPETScError(ierr));
00331 
00332     ierr = VecScatterEnd (ctx, static_cast<const Vec &>(v), vector,
00333                           INSERT_VALUES, SCATTER_FORWARD);
00334 
00335 #endif
00336     AssertThrow (ierr == 0, ExcPETScError(ierr));
00337 
00338 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
00339     ierr = VecScatterDestroy (ctx);
00340 #else
00341     ierr = VecScatterDestroy (&ctx);
00342 #endif
00343     AssertThrow (ierr == 0, ExcPETScError(ierr));
00344 #endif
00345 
00346     return *this;
00347   }
00348 
00349 
00350 
00351   template <typename number>
00352   inline
00353   Vector &
00354   Vector::operator = (const ::Vector<number> &v)
00355   {
00356     reinit (v.size(), true);
00357                                      // the following isn't necessarily fast,
00358                                      // but this is due to the fact that PETSc
00359                                      // doesn't offer an inlined access
00360                                      // operator.
00361                                      //
00362                                      // if someone wants to contribute some
00363                                      // code: to make this code faster, one
00364                                      // could either first convert all values
00365                                      // to PetscScalar, and then set them all
00366                                      // at once using VecSetValues. This has
00367                                      // the drawback that it could take quite
00368                                      // some memory, if the vector is large,
00369                                      // and it would in addition allocate
00370                                      // memory on the heap, which is
00371                                      // expensive. an alternative would be to
00372                                      // split the vector into chunks of, say,
00373                                      // 128 elements, convert a chunk at a
00374                                      // time and set it in the output vector
00375                                      // using VecSetValues. since 128 elements
00376                                      // is small enough, this could easily be
00377                                      // allocated on the stack (as a local
00378                                      // variable) which would make the whole
00379                                      // thing much more efficient.
00380                                      //
00381                                      // a second way to make things faster is
00382                                      // for the special case that
00383                                      // number==PetscScalar. we could then
00384                                      // declare a specialization of this
00385                                      // template, and omit the conversion. the
00386                                      // problem with this is that the best we
00387                                      // can do is to use VecSetValues, but
00388                                      // this isn't very efficient either: it
00389                                      // wants to see an array of indices,
00390                                      // which in this case a) again takes up a
00391                                      // whole lot of memory on the heap, and
00392                                      // b) is totally dumb since its content
00393                                      // would simply be the sequence
00394                                      // 0,1,2,3,...,n. the best of all worlds
00395                                      // would probably be a function in Petsc
00396                                      // that would take a pointer to an array
00397                                      // of PetscScalar values and simply copy
00398                                      // n elements verbatim into the vector...
00399     for (unsigned int i=0; i<v.size(); ++i)
00400       (*this)(i) = v(i);
00401 
00402     compress ();
00403 
00404     return *this;
00405   }
00406 #endif // DOXYGEN
00407 
00408 }
00409 
00410 DEAL_II_NAMESPACE_CLOSE
00411 
00412 #endif // DEAL_II_USE_PETSC
00413 
00414 /*----------------------------   petsc_vector.h     ---------------------------*/
00415 
00416 #endif
00417 /*----------------------------   petsc_vector.h     ---------------------------*/
00418 
 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