petsc_parallel_block_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_petsc_parallel_block_vector_h
#define dealii_petsc_parallel_block_vector_h


#include <deal.II/base/config.h>

#ifdef DEAL_II_WITH_PETSC

#  include <deal.II/lac/petsc_parallel_vector.h>
#  include <deal.II/lac/block_indices.h>
#  include <deal.II/lac/block_vector_base.h>
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/vector_type_traits.h>

DEAL_II_NAMESPACE_OPEN


namespace PETScWrappers
{
  // forward declaration
  class BlockVector;

  namespace MPI
  {

    class BlockVector : public BlockVectorBase<Vector>
    {
    public:
      typedef BlockVectorBase<Vector> BaseClass;

      typedef BaseClass::BlockType  BlockType;

      typedef BaseClass::value_type      value_type;
      typedef BaseClass::pointer         pointer;
      typedef BaseClass::const_pointer   const_pointer;
      typedef BaseClass::reference       reference;
      typedef BaseClass::const_reference const_reference;
      typedef BaseClass::size_type       size_type;
      typedef BaseClass::iterator        iterator;
      typedef BaseClass::const_iterator  const_iterator;

      BlockVector () = default;

      explicit BlockVector (const unsigned int  n_blocks,
                            const MPI_Comm     &communicator,
                            const size_type     block_size,
                            const size_type     local_size);

      BlockVector (const BlockVector  &V);

      BlockVector (const std::vector<size_type> &block_sizes,
                   const MPI_Comm               &communicator,
                   const std::vector<size_type> &local_elements);

      explicit BlockVector (const std::vector<IndexSet> &parallel_partitioning,
                            const MPI_Comm &communicator = MPI_COMM_WORLD);

      BlockVector (const std::vector<IndexSet> &parallel_partitioning,
                   const std::vector<IndexSet> &ghost_indices,
                   const MPI_Comm &communicator);



      ~BlockVector () = default;

      BlockVector &operator= (const value_type s);

      BlockVector &
      operator= (const BlockVector &V);

      void reinit (const unsigned int  n_blocks,
                   const MPI_Comm     &communicator,
                   const size_type     block_size,
                   const size_type     local_size,
                   const bool omit_zeroing_entries = false);

      void reinit (const std::vector<size_type> &block_sizes,
                   const MPI_Comm               &communicator,
                   const std::vector<size_type> &local_sizes,
                   const bool                    omit_zeroing_entries=false);

      void reinit (const BlockVector &V,
                   const bool         omit_zeroing_entries=false);

      void reinit (const std::vector<IndexSet> &parallel_partitioning,
                   const MPI_Comm              &communicator);

      void reinit (const std::vector<IndexSet> &parallel_partitioning,
                   const std::vector<IndexSet> &ghost_entries,
                   const MPI_Comm              &communicator);

      void reinit (const unsigned int num_blocks);

      bool has_ghost_elements() const;

      const MPI_Comm &get_mpi_communicator () const;

      void swap (BlockVector &v);

      void print (std::ostream       &out,
                  const unsigned int  precision = 3,
                  const bool          scientific = true,
                  const bool          across = true) const;

      DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize);
      DeclException0 (ExcNonMatchingBlockVectors);
    };

    /*----------------------- Inline functions ----------------------------------*/

    inline
    BlockVector::BlockVector (const unsigned int  n_blocks,
                              const MPI_Comm     &communicator,
                              const size_type     block_size,
                              const size_type     local_size)
    {
      reinit (n_blocks, communicator, block_size, local_size);
    }



    inline
    BlockVector::BlockVector (const std::vector<size_type> &block_sizes,
                              const MPI_Comm     &communicator,
                              const std::vector<size_type> &local_elements)
    {
      reinit (block_sizes, communicator, local_elements, false);
    }


    inline
    BlockVector::BlockVector (const BlockVector &v)
      :
      BlockVectorBase<Vector > ()
    {
      this->components.resize (v.n_blocks());
      this->block_indices = v.block_indices;

      for (unsigned int i=0; i<this->n_blocks(); ++i)
        this->components[i] = v.components[i];
    }

    inline
    BlockVector::BlockVector (const std::vector<IndexSet> &parallel_partitioning,
                              const MPI_Comm              &communicator)
    {
      reinit(parallel_partitioning, communicator);
    }

    inline
    BlockVector::BlockVector (const std::vector<IndexSet> &parallel_partitioning,
                              const std::vector<IndexSet> &ghost_indices,
                              const MPI_Comm &communicator)
    {
      reinit(parallel_partitioning, ghost_indices, communicator);
    }

    inline
    BlockVector &
    BlockVector::operator= (const value_type s)
    {
      BaseClass::operator= (s);
      return *this;
    }

    inline
    BlockVector &
    BlockVector::operator= (const BlockVector &v)
    {
      // we only allow assignment to vectors with the same number of blocks
      // or to an empty BlockVector
      Assert (n_blocks() == 0 || n_blocks() == v.n_blocks(),
              ExcDimensionMismatch(n_blocks(), v.n_blocks()));

      if (this->n_blocks() != v.n_blocks())
        reinit(v.n_blocks());

      for (size_type i=0; i<this->n_blocks(); ++i)
        this->components[i] = v.block(i);

      collect_sizes();

      return *this;
    }



    inline
    void
    BlockVector::reinit (const unsigned int  n_blocks,
                         const MPI_Comm     &communicator,
                         const size_type     block_size,
                         const size_type     local_size,
                         const bool omit_zeroing_entries)
    {
      reinit(std::vector<size_type>(n_blocks, block_size),
             communicator,
             std::vector<size_type>(n_blocks, local_size),
             omit_zeroing_entries);
    }



    inline
    void
    BlockVector::reinit (const std::vector<size_type> &block_sizes,
                         const MPI_Comm               &communicator,
                         const std::vector<size_type> &local_sizes,
                         const bool                    omit_zeroing_entries)
    {
      this->block_indices.reinit (block_sizes);
      if (this->components.size() != this->n_blocks())
        this->components.resize(this->n_blocks());

      for (unsigned int i=0; i<this->n_blocks(); ++i)
        this->components[i].reinit(communicator, block_sizes[i],
                                   local_sizes[i], omit_zeroing_entries);
    }


    inline
    void
    BlockVector::reinit (const BlockVector &v,
                         const bool omit_zeroing_entries)
    {
      this->block_indices = v.get_block_indices();
      if (this->components.size() != this->n_blocks())
        this->components.resize(this->n_blocks());

      for (unsigned int i=0; i<this->n_blocks(); ++i)
        block(i).reinit(v.block(i), omit_zeroing_entries);
    }

    inline
    void
    BlockVector::reinit (const std::vector<IndexSet> &parallel_partitioning,
                         const MPI_Comm              &communicator)
    {
      std::vector<size_type> sizes(parallel_partitioning.size());
      for (unsigned int i=0; i<parallel_partitioning.size(); ++i)
        sizes[i] = parallel_partitioning[i].size();

      this->block_indices.reinit(sizes);
      if (this->components.size() != this->n_blocks())
        this->components.resize(this->n_blocks());

      for (unsigned int i=0; i<this->n_blocks(); ++i)
        block(i).reinit(parallel_partitioning[i], communicator);
    }

    inline
    void
    BlockVector::reinit (const std::vector<IndexSet> &parallel_partitioning,
                         const std::vector<IndexSet> &ghost_entries,
                         const MPI_Comm              &communicator)
    {
      std::vector<types::global_dof_index> sizes(parallel_partitioning.size());
      for (unsigned int i=0; i<parallel_partitioning.size(); ++i)
        sizes[i] = parallel_partitioning[i].size();

      this->block_indices.reinit(sizes);
      if (this->components.size() != this->n_blocks())
        this->components.resize(this->n_blocks());

      for (unsigned int i=0; i<this->n_blocks(); ++i)
        block(i).reinit(parallel_partitioning[i], ghost_entries[i], communicator);
    }



    inline
    const MPI_Comm &
    BlockVector::get_mpi_communicator () const
    {
      return block(0).get_mpi_communicator();
    }

    inline
    bool
    BlockVector::has_ghost_elements() const
    {
      bool ghosted=block(0).has_ghost_elements();
#ifdef DEBUG
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        Assert(block(i).has_ghost_elements()==ghosted, ExcInternalError());
#endif
      return ghosted;
    }


    inline
    void
    BlockVector::swap (BlockVector &v)
    {
      std::swap(this->components, v.components);

      ::swap (this->block_indices, v.block_indices);
    }



    inline
    void
    BlockVector::print (std::ostream       &out,
                        const unsigned int  precision,
                        const bool          scientific,
                        const bool          across) const
    {
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        {
          if (across)
            out << 'C' << i << ':';
          else
            out << "Component " << i << std::endl;
          this->components[i].print(out, precision, scientific, across);
        }
    }



    inline
    void swap (BlockVector &u,
               BlockVector &v)
    {
      u.swap (v);
    }

  }

}

namespace internal
{
  namespace LinearOperatorImplementation
  {
    template <typename> class ReinitHelper;

    template <>
    class ReinitHelper<PETScWrappers::MPI::BlockVector>
    {
    public:
      template <typename Matrix>
      static
      void reinit_range_vector (const Matrix &matrix,
                                PETScWrappers::MPI::BlockVector &v,
                                bool /*omit_zeroing_entries*/)
      {
        v.reinit(matrix.locally_owned_range_indices(), matrix.get_mpi_communicator());
      }

      template <typename Matrix>
      static
      void reinit_domain_vector(const Matrix &matrix,
                                PETScWrappers::MPI::BlockVector &v,
                                bool /*omit_zeroing_entries*/)
      {
        v.reinit(matrix.locally_owned_domain_indices(), matrix.get_mpi_communicator());
      }
    };

  } /* namespace LinearOperator */
} /* namespace internal */


template <>
struct is_serial_vector< PETScWrappers::MPI::BlockVector > : std::false_type
{
};


DEAL_II_NAMESPACE_CLOSE

#endif  // DEAL_II_WITH_PETSC

#endif