parallel_block_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2016 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__parallel_block_vector_h
#define dealii__parallel_block_vector_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/lac/block_indices.h>
#include <deal.II/lac/block_vector_base.h>
#include <deal.II/lac/parallel_vector.h>

#include <deal.II/lac/petsc_parallel_block_vector.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>


#include <cstdio>
#include <vector>

DEAL_II_NAMESPACE_OPEN


namespace parallel
{
  namespace distributed
  {

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

      typedef typename BaseClass::BlockType  BlockType;

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

      explicit BlockVector (const size_type num_blocks = 0,
                            const size_type block_size = 0);

      BlockVector (const BlockVector<Number> &V);


#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG

      template <typename OtherNumber>
      explicit
      BlockVector (const BlockVector<OtherNumber> &v);
#endif

      BlockVector (const std::vector<size_type> &block_sizes);

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

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

      ~BlockVector ();

      BlockVector &operator = (const value_type s);

      BlockVector &
      operator= (const BlockVector &V);

      template <class Number2>
      BlockVector &
      operator= (const BlockVector<Number2> &V);

      BlockVector &
      operator= (const Vector<Number> &V);

#ifdef DEAL_II_WITH_PETSC

      BlockVector<Number> &
      operator = (const PETScWrappers::MPI::BlockVector &petsc_vec);
#endif

#ifdef DEAL_II_WITH_TRILINOS

      BlockVector<Number> &
      operator = (const TrilinosWrappers::MPI::BlockVector &trilinos_vec);
#endif

      void reinit (const size_type num_blocks,
                   const size_type block_size = 0,
                   const bool omit_zeroing_entries = false);

      void reinit (const std::vector<size_type> &N,
                   const bool                    omit_zeroing_entries=false);

      template <typename Number2>
      void reinit (const BlockVector<Number2> &V,
                   const bool                 omit_zeroing_entries=false);

      void compress (::VectorOperation::values operation);

      void update_ghost_values () const;

      void zero_out_ghosts ();

      bool has_ghost_elements() const;

      bool all_zero () const;

      bool is_non_negative () const;

      template <typename Number2>
      bool operator == (const BlockVector<Number2> &v) const;

      template <typename Number2>
      bool operator != (const BlockVector<Number2> &v) const;

      template <typename Number2>
      Number operator * (const BlockVector<Number2> &V) const;

      real_type norm_sqr () const;

      Number mean_value () const;

      real_type l1_norm () const;

      real_type l2_norm () const;

      real_type lp_norm (const real_type p) const;

      real_type linfty_norm () const;

      Number add_and_dot (const Number               a,
                          const BlockVector<Number> &V,
                          const BlockVector<Number> &W);

      template <class BlockVector2>
      void scale (const BlockVector2 &v);

      void swap (BlockVector<Number> &v);

      DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize);
    };

#ifndef DOXYGEN
    /*----------------------- Inline functions ----------------------------------*/


    template <typename Number>
    inline
    BlockVector<Number>::BlockVector (const size_type n_blocks,
                                      const size_type block_size)
    {
      reinit (n_blocks, block_size);
    }



    template <typename Number>
    inline
    BlockVector<Number>::BlockVector (const std::vector<size_type> &n)
    {
      reinit (n, false);
    }


    template <typename Number>
    inline
    BlockVector<Number>::BlockVector (const std::vector<IndexSet> &local_ranges,
                                      const std::vector<IndexSet> &ghost_indices,
                                      const MPI_Comm  communicator)
    {
      std::vector<size_type> sizes(local_ranges.size());
      for (unsigned int i=0; i<local_ranges.size(); ++i)
        sizes[i] = local_ranges[i].size();

      this->block_indices.reinit(sizes);
      this->components.resize(this->n_blocks());

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


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

      this->block_indices.reinit(sizes);
      this->components.resize(this->n_blocks());

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



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

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


#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG

    template <typename Number>
    template <typename OtherNumber>
    inline
    BlockVector<Number>::BlockVector (const BlockVector<OtherNumber> &v)
    {
      reinit (v, true);
      *this = v;
    }

#endif



    template <typename Number>
    inline
    void BlockVector<Number>::reinit (const size_type n_bl,
                                      const size_type bl_sz,
                                      const bool         omit_zeroing_entries)
    {
      std::vector<size_type> n(n_bl, bl_sz);
      reinit(n, omit_zeroing_entries);
    }


    template <typename Number>
    inline
    void BlockVector<Number>::reinit (const std::vector<size_type> &n,
                                      const bool                    omit_zeroing_entries)
    {
      this->block_indices.reinit (n);
      if (this->components.size() != this->n_blocks())
        this->components.resize(this->n_blocks());

      for (size_type i=0; i<this->n_blocks(); ++i)
        this->components[i].reinit(n[i], omit_zeroing_entries);
    }



    template <typename Number>
    template <typename Number2>
    inline
    void BlockVector<Number>::reinit (const BlockVector<Number2> &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)
        this->block(i).reinit(v.block(i), omit_zeroing_entries);
    }



    template <typename Number>
    inline
    BlockVector<Number>::~BlockVector ()
    {}



    template <typename Number>
    inline
    BlockVector<Number> &
    BlockVector<Number>::operator = (const value_type s)
    {

      AssertIsFinite(s);

      BaseClass::operator = (s);
      return *this;
    }



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

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

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

      this->collect_sizes();
      return *this;
    }



    template <typename Number>
    inline
    BlockVector<Number> &
    BlockVector<Number>::operator = (const Vector<Number> &v)
    {
      BaseClass::operator = (v);
      return *this;
    }



    template <typename Number>
    template <typename Number2>
    inline
    BlockVector<Number> &
    BlockVector<Number>::operator = (const BlockVector<Number2> &v)
    {
      reinit (v, true);
      BaseClass::operator = (v);
      return *this;
    }



#ifdef DEAL_II_WITH_PETSC

    template <typename Number>
    inline
    BlockVector<Number> &
    BlockVector<Number>::operator = (const PETScWrappers::MPI::BlockVector &petsc_vec)
    {
      AssertDimension(this->n_blocks(), petsc_vec.n_blocks());
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        this->block(i) = petsc_vec.block(i);

      return *this;
    }

#endif



#ifdef DEAL_II_WITH_TRILINOS

    template <typename Number>
    inline
    BlockVector<Number> &
    BlockVector<Number>::operator = (const TrilinosWrappers::MPI::BlockVector &trilinos_vec)
    {
      AssertDimension(this->n_blocks(), trilinos_vec.n_blocks());
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        this->block(i) = trilinos_vec.block(i);

      return *this;
    }

#endif



    template <typename Number>
    inline
    void
    BlockVector<Number>::compress (::VectorOperation::values operation)
    {
      // start all requests for all blocks before finishing the transfers as
      // this saves repeated synchronizations
      for (unsigned int block=0; block<this->n_blocks(); ++block)
        this->block(block).compress_start(block*10 + 8273, operation);
      for (unsigned int block=0; block<this->n_blocks(); ++block)
        this->block(block).compress_finish(operation);
    }



    template <typename Number>
    inline
    void
    BlockVector<Number>::update_ghost_values () const
    {
      for (unsigned int block=0; block<this->n_blocks(); ++block)
        this->block(block).update_ghost_values_start(block*10 + 9923);
      for (unsigned int block=0; block<this->n_blocks(); ++block)
        this->block(block).update_ghost_values_finish();
    }



    template <typename Number>
    inline
    void
    BlockVector<Number>::zero_out_ghosts ()
    {
      for (unsigned int block=0; block<this->n_blocks(); ++block)
        this->block(block).zero_out_ghosts();
    }



    template <typename Number>
    inline
    bool
    BlockVector<Number>::has_ghost_elements () const
    {
      bool has_ghost_elements = false;
      for (unsigned int block=0; block<this->n_blocks(); ++block)
        if (this->block(block).has_ghost_elements() == true)
          has_ghost_elements = true;
      return has_ghost_elements;
    }



    template <typename Number>
    inline
    bool
    BlockVector<Number>::all_zero () const
    {
      Assert (this->n_blocks() > 0, ExcEmptyObject());

      // use int instead of bool. in order to make global reduction operations
      // work also when MPI_Init was not called, only call MPI_Allreduce
      // commands when there is more than one processor (note that reinit()
      // functions handle this case correctly through the job_supports_mpi()
      // query). this is the same in all the functions below
      int local_result = -1;
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        local_result = std::max(local_result,
                                -static_cast<int>(this->block(i).all_zero_local()));

      if (this->block(0).partitioner->n_mpi_processes() > 1)
        return -Utilities::MPI::max(local_result,
                                    this->block(0).partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    inline
    bool
    BlockVector<Number>::is_non_negative () const
    {
      Assert (this->n_blocks() > 0, ExcEmptyObject());
      int local_result = -1;
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        local_result = std::max(local_result,
                                -static_cast<int>(this->block(i).is_non_negative_local()));
      if (this->block(0).partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::max(local_result,
                                   this->block(0).partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    template <typename Number2>
    inline
    bool
    BlockVector<Number>::operator == (const BlockVector<Number2> &v) const
    {
      Assert (this->n_blocks() > 0, ExcEmptyObject());
      AssertDimension (this->n_blocks(), v.n_blocks());

      // MPI does not support bools, so use unsigned int instead. Two vectors
      // are equal if the check for non-equal fails on all processors
      unsigned int local_result = 0;
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        local_result = std::max(local_result,
                                static_cast<unsigned int>(!this->block(i).vectors_equal_local(v.block(i))));
      unsigned int result =
        this->block(0).partitioner->n_mpi_processes() > 1
        ?
        Utilities::MPI::max(local_result, this->block(0).partitioner->get_communicator())
        :
        local_result;
      return result==0;
    }



    template <typename Number>
    template <typename Number2>
    inline
    bool
    BlockVector<Number>::operator != (const BlockVector<Number2> &v) const
    {
      return !(operator == (v));
    }



    template <typename Number>
    template <typename Number2>
    inline
    Number
    BlockVector<Number>::operator * (const BlockVector<Number2> &v) const
    {
      Assert (this->n_blocks() > 0, ExcEmptyObject());
      AssertDimension (this->n_blocks(), v.n_blocks());

      Number local_result = Number();
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        local_result += this->block(i).inner_product_local(v.block(i));

      if (this->block(0).partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::sum (local_result,
                                    this->block(0).partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    inline
    typename BlockVector<Number>::real_type
    BlockVector<Number>::norm_sqr () const
    {
      Assert (this->n_blocks() > 0, ExcEmptyObject());

      real_type local_result = real_type();
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        local_result += this->block(i).norm_sqr_local();

      if (this->block(0).partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::sum (local_result,
                                    this->block(0).partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    inline
    Number
    BlockVector<Number>::mean_value () const
    {
      Assert (this->n_blocks() > 0, ExcEmptyObject());

      Number local_result = Number();
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        local_result += this->block(i).mean_value_local()*(real_type)this->block(i).partitioner->local_size();

      if (this->block(0).partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::sum (local_result,
                                    this->block(0).partitioner->get_communicator())/
               (real_type)this->size();
      else
        return local_result/(real_type)this->size();
    }



    template <typename Number>
    inline
    typename BlockVector<Number>::real_type
    BlockVector<Number>::l1_norm () const
    {
      Assert (this->n_blocks() > 0, ExcEmptyObject());

      real_type local_result = real_type();
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        local_result += this->block(i).l1_norm_local();

      if (this->block(0).partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::sum (local_result,
                                    this->block(0).partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    inline
    typename BlockVector<Number>::real_type
    BlockVector<Number>::l2_norm () const
    {
      return std::sqrt(norm_sqr());
    }



    template <typename Number>
    inline
    typename BlockVector<Number>::real_type
    BlockVector<Number>::lp_norm (const real_type p) const
    {
      Assert (this->n_blocks() > 0, ExcEmptyObject());

      real_type local_result = real_type();
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        local_result += std::pow(this->block(i).lp_norm_local(p), p);

      if (this->block(0).partitioner->n_mpi_processes() > 1)
        return std::pow (Utilities::MPI::sum(local_result,
                                             this->block(0).partitioner->get_communicator()),
                         static_cast<real_type>(1.0/p));
      else
        return std::pow (local_result, static_cast<real_type>(1.0/p));
    }



    template <typename Number>
    inline
    typename BlockVector<Number>::real_type
    BlockVector<Number>::linfty_norm () const
    {
      Assert (this->n_blocks() > 0, ExcEmptyObject());

      real_type local_result = real_type();
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        local_result = std::max(local_result, this->block(i).linfty_norm_local());

      if (this->block(0).partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::max (local_result,
                                    this->block(0).partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    inline
    Number
    BlockVector<Number>::add_and_dot (const Number               a,
                                      const BlockVector<Number> &V,
                                      const BlockVector<Number> &W)
    {
      Number local_result = Number();
      for (unsigned int i=0; i<this->n_blocks(); ++i)
        local_result += this->block(i).add_and_dot_local(a, V.block(i), W.block(i));

      if (this->block(0).partitioner->n_mpi_processes() > 1)
        return Utilities::MPI::sum (local_result,
                                    this->block(0).partitioner->get_communicator());
      else
        return local_result;
    }



    template <typename Number>
    inline
    void
    BlockVector<Number>::swap (BlockVector<Number> &v)
    {
      Assert (this->n_blocks() == v.n_blocks(),
              ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));

      for (size_type i=0; i<this->n_blocks(); ++i)
        ::swap (this->components[i], v.components[i]);
      ::swap (this->block_indices, v.block_indices);
    }



    template <typename Number>
    template <class BlockVector2>
    void BlockVector<Number>::scale (const BlockVector2 &v)
    {
      BaseClass::scale (v);
    }

#endif // DOXYGEN

  } // end of namespace distributed

} // end of namespace parallel

template <typename Number>
inline
void swap (parallel::distributed::BlockVector<Number> &u,
           parallel::distributed::BlockVector<Number> &v)
{
  u.swap (v);
}

DEAL_II_NAMESPACE_CLOSE

#endif