Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
petsc_vector.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_petsc_vector_h
17 # define dealii_petsc_vector_h
18 
19 
20 # include <deal.II/base/config.h>
21 
22 # ifdef DEAL_II_WITH_PETSC
23 
24 # include <deal.II/base/index_set.h>
25 # include <deal.II/base/subscriptor.h>
26 
27 # include <deal.II/lac/exceptions.h>
28 # include <deal.II/lac/petsc_vector_base.h>
29 # include <deal.II/lac/vector.h>
30 # include <deal.II/lac/vector_operation.h>
31 # include <deal.II/lac/vector_type_traits.h>
32 
33 DEAL_II_NAMESPACE_OPEN
34 
35 
39 namespace PETScWrappers
40 {
48  namespace MPI
49  {
158  class Vector : public VectorBase
159  {
160  public:
165 
169  Vector();
170 
187  explicit Vector(const MPI_Comm &communicator,
188  const size_type n,
189  const size_type local_size);
190 
191 
202  template <typename Number>
203  explicit Vector(const MPI_Comm & communicator,
204  const ::Vector<Number> &v,
205  const size_type local_size);
206 
207 
220  DEAL_II_DEPRECATED
221  explicit Vector(const MPI_Comm & communicator,
222  const VectorBase &v,
223  const size_type local_size);
224 
248  Vector(const IndexSet &local,
249  const IndexSet &ghost,
250  const MPI_Comm &communicator);
251 
263  explicit Vector(const IndexSet &local, const MPI_Comm &communicator);
264 
268  Vector(const Vector &v);
269 
274  virtual void
275  clear() override;
276 
281  Vector &
282  operator=(const Vector &v);
283 
290  Vector &
291  operator=(const PetscScalar s);
292 
302  template <typename number>
303  Vector &
304  operator=(const ::Vector<number> &v);
305 
322  void
323  reinit(const MPI_Comm &communicator,
324  const size_type N,
325  const size_type local_size,
326  const bool omit_zeroing_entries = false);
327 
337  void
338  reinit(const Vector &v, const bool omit_zeroing_entries = false);
339 
347  void
348  reinit(const IndexSet &local,
349  const IndexSet &ghost,
350  const MPI_Comm &communicator);
351 
359  void
360  reinit(const IndexSet &local, const MPI_Comm &communicator);
361 
366  const MPI_Comm &
367  get_mpi_communicator() const override;
368 
380  void
381  print(std::ostream & out,
382  const unsigned int precision = 3,
383  const bool scientific = true,
384  const bool across = true) const;
385 
392  bool
393  all_zero() const;
394 
395  protected:
402  virtual void
404 
405 
406 
412  virtual void
413  create_vector(const size_type n,
414  const size_type local_size,
415  const IndexSet &ghostnodes);
416 
417 
418  private:
422  MPI_Comm communicator;
423  };
424 
425 
426  // ------------------ template and inline functions -------------
427 
428 
437  inline void
439  {
440  u.swap(v);
441  }
442 
443 
444 # ifndef DOXYGEN
445 
446  template <typename number>
447  Vector::Vector(const MPI_Comm & communicator,
448  const ::Vector<number> &v,
449  const size_type local_size)
450  : communicator(communicator)
451  {
452  Vector::create_vector(v.size(), local_size);
453 
454  *this = v;
455  }
456 
457 
458 
459  inline Vector &
460  Vector::operator=(const PetscScalar s)
461  {
463 
464  return *this;
465  }
466 
467 
468 
469  template <typename number>
470  inline Vector &
471  Vector::operator=(const ::Vector<number> &v)
472  {
473  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
474 
475  // FIXME: the following isn't necessarily fast, but this is due to
476  // the fact that PETSc doesn't offer an inlined access operator.
477  //
478  // if someone wants to contribute some code: to make this code
479  // faster, one could either first convert all values to PetscScalar,
480  // and then set them all at once using VecSetValues. This has the
481  // drawback that it could take quite some memory, if the vector is
482  // large, and it would in addition allocate memory on the heap, which
483  // is expensive. an alternative would be to split the vector into
484  // chunks of, say, 128 elements, convert a chunk at a time and set it
485  // in the output vector using VecSetValues. since 128 elements is
486  // small enough, this could easily be allocated on the stack (as a
487  // local variable) which would make the whole thing much more
488  // efficient.
489  //
490  // a second way to make things faster is for the special case that
491  // number==PetscScalar. we could then declare a specialization of
492  // this template, and omit the conversion. the problem with this is
493  // that the best we can do is to use VecSetValues, but this isn't
494  // very efficient either: it wants to see an array of indices, which
495  // in this case a) again takes up a whole lot of memory on the heap,
496  // and b) is totally dumb since its content would simply be the
497  // sequence 0,1,2,3,...,n. the best of all worlds would probably be a
498  // function in Petsc that would take a pointer to an array of
499  // PetscScalar values and simply copy n elements verbatim into the
500  // vector...
501  for (size_type i = 0; i < v.size(); ++i)
502  (*this)(i) = v(i);
503 
504  compress(::VectorOperation::insert);
505 
506  return *this;
507  }
508 
509 
510 
511  inline const MPI_Comm &
512  Vector::get_mpi_communicator() const
513  {
514  return communicator;
515  }
516 
517 # endif // DOXYGEN
518  } // namespace MPI
519 } // namespace PETScWrappers
520 
521 namespace internal
522 {
523  namespace LinearOperatorImplementation
524  {
525  template <typename>
526  class ReinitHelper;
527 
532  template <>
534  {
535  public:
536  template <typename Matrix>
537  static void
538  reinit_range_vector(const Matrix & matrix,
540  bool /*omit_zeroing_entries*/)
541  {
542  v.reinit(matrix.locally_owned_range_indices(),
543  matrix.get_mpi_communicator());
544  }
545 
546  template <typename Matrix>
547  static void
548  reinit_domain_vector(const Matrix & matrix,
550  bool /*omit_zeroing_entries*/)
551  {
552  v.reinit(matrix.locally_owned_domain_indices(),
553  matrix.get_mpi_communicator());
554  }
555  };
556 
557  } // namespace LinearOperatorImplementation
558 } /* namespace internal */
559 
568 template <>
569 struct is_serial_vector<PETScWrappers::MPI::Vector> : std::false_type
570 {};
571 
572 
573 DEAL_II_NAMESPACE_CLOSE
574 
575 # endif // DEAL_II_WITH_PETSC
576 
577 #endif
578 /*------------------------- petsc_vector.h -------------------------*/
virtual void clear() override
Vector< Number > & operator=(const Number s)
VectorBase & operator=(const VectorBase &)=delete
LinearAlgebra::distributed::Vector< Number > Vector
static void reinit_range_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
static void reinit_domain_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
void swap(BlockVector &u, BlockVector &v)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
types::global_dof_index size_type
Definition: petsc_vector.h:164
Vector & operator=(const Vector &v)
unsigned int global_dof_index
Definition: types.h:89
const MPI_Comm & get_mpi_communicator() const override
void reinit(const MPI_Comm &communicator, const size_type N, const size_type local_size, const bool omit_zeroing_entries=false)
virtual void create_vector(const size_type n, const size_type local_size)
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
void swap(Vector &u, Vector &v)
Definition: petsc_vector.h:438