petsc_block_vector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2017 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_block_vector_h
#define dealii__petsc_block_vector_h


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

#ifdef DEAL_II_WITH_PETSC

#  include <deal.II/lac/petsc_vector.h>
#  include <deal.II/lac/petsc_parallel_block_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
{
  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;

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

    BlockVector (const BlockVector  &V);

    explicit BlockVector (const MPI::BlockVector &v);

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

    template <typename InputIterator>
    BlockVector (const std::vector<size_type> &n,
                 const InputIterator           first,
                 const InputIterator           end);

    ~BlockVector ();

    BlockVector &operator = (const value_type s);

    BlockVector &
    operator= (const BlockVector &V);

    BlockVector &
    operator = (const MPI::BlockVector &v);

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

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

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

    void reinit (const unsigned int num_blocks);

    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);
  } DEAL_II_DEPRECATED;

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



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



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


  inline
  BlockVector::BlockVector (const BlockVector &v)
    :
    BlockVectorBase<Vector > ()
  {
    this->components.resize (v.n_blocks());
    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 MPI::BlockVector &v)
    :
    BlockVectorBase<Vector > ()
  {
    this->components.resize (v.get_block_indices().size());
    block_indices = v.get_block_indices();

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



  template <typename InputIterator>
  BlockVector::BlockVector (const std::vector<size_type> &n,
                            const InputIterator           first,
                            const InputIterator           end)
  {
    // first set sizes of blocks, but
    // don't initialize them as we will
    // copy elements soon
    (void)end;
    reinit (n, true);
    InputIterator start = first;
    for (unsigned int b=0; b<n.size(); ++b)
      {
        InputIterator end = start;
        std::advance (end, static_cast<signed int>(n[b]));

        for (size_type i=0; i<n[b]; ++i, ++start)
          this->block(b)(i) = *start;
      }
    Assert (start == end, ExcIteratorRangeDoesNotMatchVectorSize());
  }



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



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



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



  inline
  BlockVector::~BlockVector ()
  {}


  inline
  void
  BlockVector::reinit (const unsigned int 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);
  }



  inline
  void
  BlockVector::reinit (const std::vector<size_type> &n,
                       const bool                    omit_zeroing_entries)
  {
    block_indices.reinit (n);
    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(n[i], omit_zeroing_entries);
  }


  inline
  void
  BlockVector::reinit (const BlockVector &v,
                       const bool omit_zeroing_entries)
  {
    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 unsigned int num_blocks)
  {
    reinit (num_blocks, 0, true);
  }



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

    for (unsigned int i=0; i<this->n_blocks(); ++i)
      this->components[i].swap (v.components[i]);
    ::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 LinearOperator
  {
    template <typename> class ReinitHelper;

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

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

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


template <>
struct is_serial_vector< PETScWrappers::BlockVector > : std_cxx11::true_type
{
};


DEAL_II_NAMESPACE_CLOSE

#endif  // DEAL_II_WITH_PETSC

#endif