Reference documentation for deal.II version GIT b8135fa6eb 2023-03-25 15:55:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
petsc_block_vector.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2021 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_block_vector_h
17 #define dealii_petsc_block_vector_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #ifdef DEAL_II_WITH_PETSC
23 
26 # include <deal.II/lac/exceptions.h>
29 
31 
32 
33 namespace PETScWrappers
34 {
35  // forward declaration
36  class BlockVector;
37 
38  namespace MPI
39  {
61  class BlockVector : public BlockVectorBase<Vector>
62  {
63  public:
68 
73 
85 
89  BlockVector();
90 
97  explicit BlockVector(const unsigned int n_blocks,
98  const MPI_Comm & communicator,
99  const size_type block_size,
101 
106  BlockVector(const BlockVector &V);
107 
115  BlockVector(const std::vector<size_type> &block_sizes,
116  const MPI_Comm & communicator,
117  const std::vector<size_type> &local_elements);
118 
123  explicit BlockVector(const std::vector<IndexSet> &parallel_partitioning,
124  const MPI_Comm &communicator = MPI_COMM_WORLD);
125 
129  BlockVector(const std::vector<IndexSet> &parallel_partitioning,
130  const std::vector<IndexSet> &ghost_indices,
131  const MPI_Comm & communicator);
132 
139  explicit BlockVector(Vec v);
140 
144  ~BlockVector() override;
145 
150  BlockVector &
151  operator=(const value_type s);
152 
156  BlockVector &
157  operator=(const BlockVector &V);
158 
164  void
165  reinit(Vec v);
166 
176  void
177  reinit(const unsigned int n_blocks,
178  const MPI_Comm & communicator,
179  const size_type block_size,
181  const bool omit_zeroing_entries = false);
182 
203  void
204  reinit(const std::vector<size_type> &block_sizes,
205  const MPI_Comm & communicator,
206  const std::vector<size_type> &locally_owned_sizes,
207  const bool omit_zeroing_entries = false);
208 
223  void
224  reinit(const BlockVector &V, const bool omit_zeroing_entries = false);
225 
230  void
231  reinit(const std::vector<IndexSet> &parallel_partitioning,
232  const MPI_Comm & communicator);
233 
237  void
238  reinit(const std::vector<IndexSet> &parallel_partitioning,
239  const std::vector<IndexSet> &ghost_entries,
240  const MPI_Comm & communicator);
241 
249  void
250  collect_sizes();
251 
258  void
259  reinit(const unsigned int num_blocks);
260 
264  bool
265  has_ghost_elements() const;
266 
271  const MPI_Comm &
272  get_mpi_communicator() const;
273 
281  operator const Vec &() const;
282 
288  Vec &
289  petsc_vector();
290 
308  void
309  swap(BlockVector &v);
310 
314  void
315  print(std::ostream & out,
316  const unsigned int precision = 3,
317  const bool scientific = true,
318  const bool across = true) const;
319 
328 
329  private:
330  void
331  setup_nest_vec();
332 
334  };
335 
338  /*--------------------- Inline functions --------------------------------*/
339 
341  : petsc_nest_vector(nullptr)
342  {}
343 
344 
345 
346  inline BlockVector::BlockVector(const unsigned int n_blocks,
347  const MPI_Comm & communicator,
348  const size_type block_size,
349  const size_type locally_owned_size)
350  : petsc_nest_vector(nullptr)
351  {
352  reinit(n_blocks, communicator, block_size, locally_owned_size);
353  }
354 
355 
356 
358  const std::vector<size_type> &block_sizes,
359  const MPI_Comm & communicator,
360  const std::vector<size_type> &local_elements)
361  : petsc_nest_vector(nullptr)
362  {
363  reinit(block_sizes, communicator, local_elements, false);
364  }
365 
366 
369  , petsc_nest_vector(nullptr)
370  {
371  this->block_indices = v.block_indices;
372 
373  this->components.resize(this->n_blocks());
374  for (unsigned int i = 0; i < this->n_blocks(); ++i)
375  this->components[i] = v.components[i];
376 
377  this->collect_sizes();
378  }
379 
381  const std::vector<IndexSet> &parallel_partitioning,
382  const MPI_Comm & communicator)
383  : petsc_nest_vector(nullptr)
384  {
385  reinit(parallel_partitioning, communicator);
386  }
387 
389  const std::vector<IndexSet> &parallel_partitioning,
390  const std::vector<IndexSet> &ghost_indices,
391  const MPI_Comm & communicator)
392  : petsc_nest_vector(nullptr)
393  {
394  reinit(parallel_partitioning, ghost_indices, communicator);
395  }
396 
399  , petsc_nest_vector(nullptr)
400  {
401  this->reinit(v);
402  }
403 
404  inline BlockVector &
406  {
408  return *this;
409  }
410 
411  inline BlockVector &
413  {
414  // we only allow assignment to vectors with the same number of blocks
415  // or to an empty BlockVector
416  Assert(this->n_blocks() == 0 || this->n_blocks() == v.n_blocks(),
417  ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
418 
419  if (this->n_blocks() != v.n_blocks())
420  this->block_indices = v.block_indices;
421 
422  this->components.resize(this->n_blocks());
423  for (unsigned int i = 0; i < this->n_blocks(); ++i)
424  this->components[i] = v.components[i];
425 
426  this->collect_sizes();
427 
428  return *this;
429  }
430 
431 
432 
433  inline void
434  BlockVector::reinit(const unsigned int n_blocks,
435  const MPI_Comm & communicator,
436  const size_type block_size,
437  const size_type locally_owned_size,
438  const bool omit_zeroing_entries)
439  {
440  reinit(std::vector<size_type>(n_blocks, block_size),
441  communicator,
442  std::vector<size_type>(n_blocks, locally_owned_size),
443  omit_zeroing_entries);
444  }
445 
446 
447 
448  inline void
449  BlockVector::reinit(const std::vector<size_type> &block_sizes,
450  const MPI_Comm & communicator,
451  const std::vector<size_type> &locally_owned_sizes,
452  const bool omit_zeroing_entries)
453  {
454  this->block_indices.reinit(block_sizes);
455 
456  this->components.resize(this->n_blocks());
457  for (unsigned int i = 0; i < this->n_blocks(); ++i)
458  this->components[i].reinit(communicator,
459  block_sizes[i],
460  locally_owned_sizes[i],
461  omit_zeroing_entries);
462 
463  this->collect_sizes();
464  }
465 
466 
467  inline void
468  BlockVector::reinit(const BlockVector &v, const bool omit_zeroing_entries)
469  {
470  if (this->n_blocks() != v.n_blocks())
471  this->block_indices = v.get_block_indices();
472 
473  this->components.resize(this->n_blocks());
474  for (unsigned int i = 0; i < this->n_blocks(); ++i)
475  this->components[i].reinit(v.components[i], omit_zeroing_entries);
476 
477  this->collect_sizes();
478  }
479 
480  inline void
481  BlockVector::reinit(const std::vector<IndexSet> &parallel_partitioning,
482  const MPI_Comm & communicator)
483  {
484  // update the number of blocks
485  this->block_indices.reinit(parallel_partitioning.size(), 0);
486 
487  // initialize each block
488  this->components.resize(this->n_blocks());
489  for (unsigned int i = 0; i < this->n_blocks(); ++i)
490  this->components[i].reinit(parallel_partitioning[i], communicator);
491 
492  // update block_indices content
493  this->collect_sizes();
494  }
495 
496  inline void
497  BlockVector::reinit(const std::vector<IndexSet> &parallel_partitioning,
498  const std::vector<IndexSet> &ghost_entries,
499  const MPI_Comm & communicator)
500  {
501  AssertDimension(parallel_partitioning.size(), ghost_entries.size());
502 
503  // update the number of blocks
504  this->block_indices.reinit(parallel_partitioning.size(), 0);
505 
506  // initialize each block
507  this->components.resize(this->n_blocks());
508  for (unsigned int i = 0; i < this->n_blocks(); ++i)
509  this->components[i].reinit(parallel_partitioning[i],
510  ghost_entries[i],
511  communicator);
512 
513  // update block_indices content
514  this->collect_sizes();
515  }
516 
517 
518 
519  inline const MPI_Comm &
521  {
522  static MPI_Comm comm = PETSC_COMM_SELF;
523  MPI_Comm pcomm =
524  PetscObjectComm(reinterpret_cast<PetscObject>(petsc_nest_vector));
525  if (pcomm != MPI_COMM_NULL)
526  comm = pcomm;
527  return comm;
528  }
529 
530  inline bool
532  {
533  bool ghosted = block(0).has_ghost_elements();
534 # ifdef DEBUG
535  for (unsigned int i = 0; i < this->n_blocks(); ++i)
536  Assert(block(i).has_ghost_elements() == ghosted, ExcInternalError());
537 # endif
538  return ghosted;
539  }
540 
541 
542  inline void
544  {
545  std::swap(this->components, v.components);
547 
549  }
550 
551 
552 
553  inline void
554  BlockVector::print(std::ostream & out,
555  const unsigned int precision,
556  const bool scientific,
557  const bool across) const
558  {
559  for (unsigned int i = 0; i < this->n_blocks(); ++i)
560  {
561  if (across)
562  out << 'C' << i << ':';
563  else
564  out << "Component " << i << std::endl;
565  this->components[i].print(out, precision, scientific, across);
566  }
567  }
568 
569 
570 
578  inline void
580  {
581  u.swap(v);
582  }
583  } // namespace MPI
584 
585 } // namespace PETScWrappers
586 
587 namespace internal
588 {
589  namespace LinearOperatorImplementation
590  {
591  template <typename>
592  class ReinitHelper;
593 
598  template <>
600  {
601  public:
602  template <typename Matrix>
603  static void
604  reinit_range_vector(const Matrix & matrix,
606  bool /*omit_zeroing_entries*/)
607  {
608  v.reinit(matrix.locally_owned_range_indices(),
609  matrix.get_mpi_communicator());
610  }
611 
612  template <typename Matrix>
613  static void
614  reinit_domain_vector(const Matrix & matrix,
616  bool /*omit_zeroing_entries*/)
617  {
618  v.reinit(matrix.locally_owned_domain_indices(),
619  matrix.get_mpi_communicator());
620  }
621  };
622 
623  } // namespace LinearOperatorImplementation
624 } /* namespace internal */
625 
626 
630 template <>
632 {};
633 
634 
636 
637 #endif // DEAL_II_WITH_PETSC
638 
639 #endif
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
const BlockIndices & get_block_indices() const
::internal::BlockVectorIterators::Iterator< BlockVectorBase, false > iterator
unsigned int n_blocks() const
const value_type * const_pointer
typename BlockType::const_reference const_reference
types::global_dof_index size_type
typename BlockType::value_type value_type
BlockVectorBase & operator=(const value_type s)
std::vector< Vector > components
::internal::BlockVectorIterators::Iterator< BlockVectorBase, true > const_iterator
typename BlockType::reference reference
BlockType & block(const unsigned int i)
std::size_t locally_owned_size() const
BaseClass::value_type value_type
void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const bool across=true) const
BaseClass::const_reference const_reference
void swap(BlockVector &u, BlockVector &v)
BlockVector & operator=(const value_type s)
BaseClass::const_pointer const_pointer
const MPI_Comm & get_mpi_communicator() const
Definition: vector.h:109
bool has_ghost_elements() const
static void reinit_range_vector(const Matrix &matrix, PETScWrappers::MPI::BlockVector &v, bool)
static void reinit_domain_vector(const Matrix &matrix, PETScWrappers::MPI::BlockVector &v, bool)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
#define DeclException0(Exception0)
Definition: exceptions.h:465
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
static ::ExceptionBase & ExcNonMatchingBlockVectors()
static ::ExceptionBase & ExcIteratorRangeDoesNotMatchVectorSize()
constexpr int block_size
Definition: cuda_size.h:29
@ matrix
Contents is actually a matrix.
static const char V
BlockVector< double > BlockVector
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
void swap(BlockVector &u, BlockVector &v)
const MPI_Comm & comm