Reference documentation for deal.II version 8.1.0
petsc_parallel_vector.h
1 // ---------------------------------------------------------------------
2 // @f$Id: petsc_parallel_vector.h 30231 2013-08-05 22:07:43Z heister @f$
3 //
4 // Copyright (C) 2004 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__petsc_parallel_vector_h
18 #define __deal2__petsc_parallel_vector_h
19 
20 
21 #include <deal.II/base/config.h>
22 
23 #ifdef DEAL_II_WITH_PETSC
24 
25 # include <deal.II/base/subscriptor.h>
26 # include <deal.II/lac/exceptions.h>
27 # include <deal.II/lac/vector.h>
28 # include <deal.II/lac/petsc_vector_base.h>
29 # include <deal.II/base/index_set.h>
30 
32 
33 
34 
35 // forward declaration
36 template <typename> class Vector;
37 class IndexSet;
38 
39 
43 namespace PETScWrappers
44 {
52  namespace MPI
53  {
54 
158  class Vector : public VectorBase
159  {
160  public:
165 
177  static const bool supports_distributed_data = true;
178 
183  Vector ();
184 
209  explicit Vector (const MPI_Comm &communicator,
210  const size_type n,
211  const size_type local_size);
212 
213 
229  template <typename Number>
230  explicit Vector (const MPI_Comm &communicator,
231  const ::Vector<Number> &v,
232  const size_type local_size);
233 
234 
249  explicit Vector (const MPI_Comm &communicator,
250  const VectorBase &v,
251  const size_type local_size);
252 
253 
275  explicit Vector (const MPI_Comm &communicator,
276  const IndexSet &local,
277  const IndexSet &ghost) DEAL_II_DEPRECATED;
278 
300  Vector (const IndexSet &local,
301  const IndexSet &ghost,
302  const MPI_Comm &communicator);
303 
309  explicit Vector (const MPI_Comm &communicator,
310  const IndexSet &local) DEAL_II_DEPRECATED;
316  explicit Vector (const IndexSet &local,
317  const MPI_Comm &communicator);
318 
325  Vector &operator = (const Vector &v);
326 
327 
364 
374  Vector &operator = (const PetscScalar s);
375 
391  template <typename number>
392  Vector &operator = (const ::Vector<number> &v);
393 
420  void reinit (const MPI_Comm &communicator,
421  const size_type N,
422  const size_type local_size,
423  const bool fast = false);
424 
439  void reinit (const Vector &v,
440  const bool fast = false);
441 
447  void reinit (const MPI_Comm &communicator,
448  const IndexSet &local,
449  const IndexSet &ghost) DEAL_II_DEPRECATED;
455  void reinit (const IndexSet &local,
456  const IndexSet &ghost,
457  const MPI_Comm &communicator);
458 
464  void reinit (const MPI_Comm &communicator,
465  const IndexSet &local) DEAL_II_DEPRECATED;
471  void reinit (const IndexSet &local,
472  const MPI_Comm &communicator);
473 
479  const MPI_Comm &get_mpi_communicator () const;
480 
500  void print (std::ostream &out,
501  const unsigned int precision = 3,
502  const bool scientific = true,
503  const bool across = true) const;
504 
505  protected:
515  virtual void create_vector (const size_type n,
516  const size_type local_size);
517 
518 
519 
528  virtual void create_vector (const size_type n,
529  const size_type local_size,
530  const IndexSet &ghostnodes);
531 
532 
533  private:
538  MPI_Comm communicator;
539  };
540 
541 
542 // ------------------ template and inline functions -------------
543 
544 
553  inline
554  void swap (Vector &u, Vector &v)
555  {
556  u.swap (v);
557  }
558 
559 
560 #ifndef DOXYGEN
561 
562  template <typename number>
563  Vector::Vector (const MPI_Comm &communicator,
564  const ::Vector<number> &v,
565  const size_type local_size)
566  :
567  communicator (communicator)
568  {
569  Vector::create_vector (v.size(), local_size);
570 
571  *this = v;
572  }
573 
574 
575 
576  inline
577  Vector &
578  Vector::operator = (const PetscScalar s)
579  {
581 
582  return *this;
583  }
584 
585 
586 
587  inline
588  Vector &
589  Vector::operator = (const Vector &v)
590  {
591  if (v.size()==0)
592  {
593  // this happens if v has not been initialized to something useful:
594  // Vector x,v;x=v;
595  // we skip the code below and create a simple serial vector of
596  // length 0
597 
598  int ierr;
599 #if DEAL_II_PETSC_VERSION_LT(3,2,0)
600  ierr = VecDestroy (vector);
601 #else
602  ierr = VecDestroy (&vector);
603 #endif
604  AssertThrow (ierr == 0, ExcPETScError(ierr));
605 
606  const int n = 0;
607  ierr = VecCreateSeq (PETSC_COMM_SELF, n, &vector);
608  AssertThrow (ierr == 0, ExcPETScError(ierr));
609  ghosted = false;
611  return *this;
612  }
613 
614  // if the vectors have different sizes,
615  // then first resize the present one
616  if (size() != v.size())
617  {
618  if (v.has_ghost_elements())
619  reinit( v.locally_owned_elements(), v.ghost_indices, v.communicator);
620  else
621  reinit (v.communicator, v.size(), v.local_size(), true);
622  }
623 
624  const int ierr = VecCopy (v.vector, vector);
625  AssertThrow (ierr == 0, ExcPETScError(ierr));
626 
627  if (has_ghost_elements())
628  {
629  int ierr;
630 
631  ierr = VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD);
632  AssertThrow (ierr == 0, ExcPETScError(ierr));
633  ierr = VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD);
634  AssertThrow (ierr == 0, ExcPETScError(ierr));
635  }
636  return *this;
637  }
638 
639 
640 
641  template <typename number>
642  inline
643  Vector &
644  Vector::operator = (const ::Vector<number> &v)
645  {
646  Assert (size() == v.size(),
647  ExcDimensionMismatch (size(), v.size()));
648 
649  // the following isn't necessarily fast,
650  // but this is due to the fact that PETSc
651  // doesn't offer an inlined access
652  // operator.
653  //
654  // if someone wants to contribute some
655  // code: to make this code faster, one
656  // could either first convert all values
657  // to PetscScalar, and then set them all
658  // at once using VecSetValues. This has
659  // the drawback that it could take quite
660  // some memory, if the vector is large,
661  // and it would in addition allocate
662  // memory on the heap, which is
663  // expensive. an alternative would be to
664  // split the vector into chunks of, say,
665  // 128 elements, convert a chunk at a
666  // time and set it in the output vector
667  // using VecSetValues. since 128 elements
668  // is small enough, this could easily be
669  // allocated on the stack (as a local
670  // variable) which would make the whole
671  // thing much more efficient.
672  //
673  // a second way to make things faster is
674  // for the special case that
675  // number==PetscScalar. we could then
676  // declare a specialization of this
677  // template, and omit the conversion. the
678  // problem with this is that the best we
679  // can do is to use VecSetValues, but
680  // this isn't very efficient either: it
681  // wants to see an array of indices,
682  // which in this case a) again takes up a
683  // whole lot of memory on the heap, and
684  // b) is totally dumb since its content
685  // would simply be the sequence
686  // 0,1,2,3,...,n. the best of all worlds
687  // would probably be a function in Petsc
688  // that would take a pointer to an array
689  // of PetscScalar values and simply copy
690  // n elements verbatim into the vector...
691  for (size_type i=0; i<v.size(); ++i)
692  (*this)(i) = v(i);
693 
694  compress (::VectorOperation::insert);
695 
696  return *this;
697  }
698 
699 
700 
701  inline
702  const MPI_Comm &
704  {
705  return communicator;
706  }
707 
708 #endif // DOXYGEN
709  }
710 }
711 
714 DEAL_II_NAMESPACE_CLOSE
715 
716 #endif // DEAL_II_WITH_PETSC
717 
718 /*---------------------------- petsc_parallel_vector.h ---------------------------*/
719 
720 #endif
721 /*---------------------------- petsc_parallel_vector.h ---------------------------*/
IndexSet locally_owned_elements() const
types::global_dof_index size_type
Vector & operator=(const Vector &v)
bool has_ghost_elements() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
VectorBase & operator=(const PetscScalar s)
void clear()
Definition: index_set.h:662
void compress() DEAL_II_DEPRECATED
size_type size() const
unsigned int global_dof_index
Definition: types.h:100
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t size() const
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
static const bool supports_distributed_data
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
void reinit(const MPI_Comm &communicator, const size_type N, const size_type local_size, const bool fast=false)
const MPI_Comm & get_mpi_communicator() const
virtual void create_vector(const size_type n, const size_type local_size)
size_type local_size() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void swap(Vector &u, Vector &v)