block_vector_base.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2004 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_block_vector_base_h
#define dealii_block_vector_base_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/numbers.h>
 
#include <deal.II/lac/block_indices.h>
#include <deal.II/lac/exceptions.h>
#include <deal.II/lac/read_vector.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_operation.h>
 
#include <cmath>
#include <cstddef>
#include <iterator>
#include <type_traits>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
 
namespace internal
{
  template <typename T>
  using has_block_t = decltype(std::declval<const T>().block(0));
 
  template <typename T>
  constexpr bool has_block = internal::is_supported_operation<has_block_t, T>;
 
  template <typename T>
  using has_n_blocks_t = decltype(std::declval<const T>().n_blocks());
 
  template <typename T>
  constexpr bool has_n_blocks =
    internal::is_supported_operation<has_n_blocks_t, T>;
 
  template <typename T>
  constexpr bool is_block_vector = has_block<T> && has_n_blocks<T>;
} // namespace internal
 
template <typename VectorType>
struct IsBlockVector
{
public:
  static const bool value = internal::is_block_vector<VectorType>;
};
 
 
// instantiation of the static member
template <typename VectorType>
const bool IsBlockVector<VectorType>::value;
 
 
 
namespace internal
{
  namespace BlockVectorIterators
  {
    template <typename BlockVectorType, bool Constness>
    class Iterator
    {
    public:
      using size_type = types::global_dof_index;
 
      using value_type =
        std::conditional_t<Constness,
                           const typename BlockVectorType::value_type,
                           typename BlockVectorType::value_type>;
 
      using iterator_category = std::random_access_iterator_tag;
      using difference_type   = std::ptrdiff_t;
      using reference         = typename BlockVectorType::reference;
      using pointer           = value_type *;
 
      using dereference_type =
        std::conditional_t<Constness,
                           value_type,
                           typename BlockVectorType::BlockType::reference>;
 
      using BlockVector =
        std::conditional_t<Constness, const BlockVectorType, BlockVectorType>;
 
      Iterator(BlockVector &parent, const size_type global_index);
 
      Iterator(const Iterator<BlockVectorType, !Constness> &c);
 
 
      Iterator(const Iterator &c);
 
    private:
      Iterator(BlockVector    &parent,
               const size_type global_index,
               const size_type current_block,
               const size_type index_within_block,
               const size_type next_break_forward,
               const size_type next_break_backward);
 
    public:
      Iterator &
      operator=(const Iterator &c);
 
      dereference_type
      operator*() const;
 
      dereference_type
      operator[](const difference_type d) const;
 
      Iterator &
      operator++();
 
      Iterator
      operator++(int);
 
      Iterator &
      operator--();
 
      Iterator
      operator--(int);
 
      template <bool OtherConstness>
      bool
      operator==(const Iterator<BlockVectorType, OtherConstness> &i) const;
 
      template <bool OtherConstness>
      bool
      operator!=(const Iterator<BlockVectorType, OtherConstness> &i) const;
 
      template <bool OtherConstness>
      bool
      operator<(const Iterator<BlockVectorType, OtherConstness> &i) const;
 
      template <bool OtherConstness>
      bool
      operator<=(const Iterator<BlockVectorType, OtherConstness> &i) const;
 
      template <bool OtherConstness>
      bool
      operator>(const Iterator<BlockVectorType, OtherConstness> &i) const;
 
      template <bool OtherConstness>
      bool
      operator>=(const Iterator<BlockVectorType, OtherConstness> &i) const;
 
      template <bool OtherConstness>
      difference_type
      operator-(const Iterator<BlockVectorType, OtherConstness> &i) const;
 
      Iterator
      operator+(const difference_type &d) const;
 
      Iterator
      operator-(const difference_type &d) const;
 
      Iterator &
      operator+=(const difference_type &d);
 
      Iterator &
      operator-=(const difference_type &d);
 
      DeclExceptionMsg(ExcPointerToDifferentVectors,
                       "Your program tried to compare iterators pointing to "
                       "different block vectors. There is no reasonable way "
                       "to do this.");
 
    private:
      BlockVector *parent;
 
      size_type global_index;
 
      unsigned int current_block;
      size_type    index_within_block;
 
      size_type next_break_forward;
      size_type next_break_backward;
 
      void
      move_forward();
 
      void
      move_backward();
 
      // Mark all other instances of this template as friends.
      template <typename, bool>
      friend class Iterator;
    };
  } // namespace BlockVectorIterators
} // namespace internal
 
 
template <typename VectorType>
class BlockVectorBase : public ReadVector<typename VectorType::value_type>
{
public:
  using BlockType = VectorType;
 
  /*
   * Declare standard types used in
   * all containers. These types
   * parallel those in the
   * <tt>C++</tt> standard
   * libraries
   * <tt>std::vector<...></tt>
   * class. This includes iterator
   * types.
   */
  using value_type    = typename BlockType::value_type;
  using pointer       = value_type *;
  using const_pointer = const value_type *;
  using iterator =
    ::internal::BlockVectorIterators::Iterator<BlockVectorBase, false>;
  using const_iterator =
    ::internal::BlockVectorIterators::Iterator<BlockVectorBase, true>;
  using reference       = typename BlockType::reference;
  using const_reference = typename BlockType::const_reference;
  using size_type       = types::global_dof_index;
 
  using real_type = typename BlockType::real_type;
 
  BlockVectorBase() = default;
 
  BlockVectorBase(const BlockVectorBase & /*V*/) = default;
 
  BlockVectorBase(BlockVectorBase && /*V*/) noexcept = default;
 
  void
  collect_sizes();
 
  void
  compress(VectorOperation::values operation);
 
  BlockType &
  block(const unsigned int i);
 
  const BlockType &
  block(const unsigned int i) const;
 
  const BlockIndices &
  get_block_indices() const;
 
  unsigned int
  n_blocks() const;
 
  virtual size_type
  size() const override;
 
  std::size_t
  locally_owned_size() const;
 
  IndexSet
  locally_owned_elements() const;
 
  iterator
  begin();
 
  const_iterator
  begin() const;
 
  iterator
  end();
 
  const_iterator
  end() const;
 
  value_type
  operator()(const size_type i) const;
 
  reference
  operator()(const size_type i);
 
  value_type
  operator[](const size_type i) const;
 
  reference
  operator[](const size_type i);
 
  template <typename OtherNumber>
  void
  extract_subvector_to(const std::vector<size_type> &indices,
                       std::vector<OtherNumber>     &values) const;
 
  virtual void
  extract_subvector_to(const ArrayView<const types::global_dof_index> &indices,
                       const ArrayView<value_type> &entries) const override;
 
  template <typename ForwardIterator, typename OutputIterator>
  void
  extract_subvector_to(ForwardIterator       indices_begin,
                       const ForwardIterator indices_end,
                       OutputIterator        values_begin) const;
 
  BlockVectorBase &
  operator=(const value_type s);
 
  BlockVectorBase &
  operator=(const BlockVectorBase &V);
 
  BlockVectorBase &
  operator=(BlockVectorBase && /*V*/) = default; // NOLINT
 
  template <typename VectorType2>
  BlockVectorBase &
  operator=(const BlockVectorBase<VectorType2> &V);
 
  BlockVectorBase &
  operator=(const VectorType &v);
 
  template <typename VectorType2>
  bool
  operator==(const BlockVectorBase<VectorType2> &v) const;
 
  value_type
  operator*(const BlockVectorBase &V) const;
 
  real_type
  norm_sqr() const;
 
  value_type
  mean_value() const;
 
  real_type
  l1_norm() const;
 
  real_type
  l2_norm() const;
 
  real_type
  linfty_norm() const;
 
  value_type
  add_and_dot(const value_type       a,
              const BlockVectorBase &V,
              const BlockVectorBase &W);
 
  bool
  in_local_range(const size_type global_index) const;
 
  bool
  all_zero() const;
 
  bool
  is_non_negative() const;
 
  BlockVectorBase &
  operator+=(const BlockVectorBase &V);
 
  BlockVectorBase &
  operator-=(const BlockVectorBase &V);
 
 
  template <typename Number>
  void
  add(const std::vector<size_type> &indices, const std::vector<Number> &values);
 
  template <typename Number>
  void
  add(const std::vector<size_type> &indices, const Vector<Number> &values);
 
  template <typename Number>
  void
  add(const size_type  n_elements,
      const size_type *indices,
      const Number    *values);
 
  void
  add(const value_type s);
 
  void
  add(const value_type a, const BlockVectorBase &V);
 
  void
  add(const value_type       a,
      const BlockVectorBase &V,
      const value_type       b,
      const BlockVectorBase &W);
 
  void
  sadd(const value_type s, const BlockVectorBase &V);
 
  void
  sadd(const value_type s, const value_type a, const BlockVectorBase &V);
 
  void
  sadd(const value_type       s,
       const value_type       a,
       const BlockVectorBase &V,
       const value_type       b,
       const BlockVectorBase &W);
 
  void
  sadd(const value_type       s,
       const value_type       a,
       const BlockVectorBase &V,
       const value_type       b,
       const BlockVectorBase &W,
       const value_type       c,
       const BlockVectorBase &X);
 
  BlockVectorBase &
  operator*=(const value_type factor);
 
  BlockVectorBase &
  operator/=(const value_type factor);
 
  template <class BlockVector2>
  void
  scale(const BlockVector2 &v);
 
  template <class BlockVector2>
  void
  equ(const value_type a, const BlockVector2 &V);
 
  void
  update_ghost_values() const;
 
  MPI_Comm
  get_mpi_communicator() const;
 
  std::size_t
  memory_consumption() const;
 
protected:
  std::vector<VectorType> components;
 
  BlockIndices block_indices;
 
  // Make the iterator class a friend.
  template <typename N, bool C>
  friend class ::internal::BlockVectorIterators::Iterator;
 
  template <typename>
  friend class BlockVectorBase;
};
 
 
/*----------------------- Inline functions ----------------------------------*/
 
 
#ifndef DOXYGEN
namespace internal
{
  namespace BlockVectorIterators
  {
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness>::Iterator(
      const Iterator<BlockVectorType, Constness> &c)
      : parent(c.parent)
      , global_index(c.global_index)
      , current_block(c.current_block)
      , index_within_block(c.index_within_block)
      , next_break_forward(c.next_break_forward)
      , next_break_backward(c.next_break_backward)
    {}
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness>::Iterator(
      const Iterator<BlockVectorType, !Constness> &c)
      : parent(c.parent)
      , global_index(c.global_index)
      , current_block(c.current_block)
      , index_within_block(c.index_within_block)
      , next_break_forward(c.next_break_forward)
      , next_break_backward(c.next_break_backward)
    {
      // Only permit copy-constructing const iterators from non-const
      // iterators, and not vice versa (i.e., Constness must always be
      // true).
      static_assert(Constness == true,
                    "Constructing a non-const iterator from a const iterator "
                    "does not make sense.");
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness>::Iterator(
      BlockVector    &parent,
      const size_type global_index,
      const size_type current_block,
      const size_type index_within_block,
      const size_type next_break_forward,
      const size_type next_break_backward)
      : parent(&parent)
      , global_index(global_index)
      , current_block(current_block)
      , index_within_block(index_within_block)
      , next_break_forward(next_break_forward)
      , next_break_backward(next_break_backward)
    {}
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness> &
    Iterator<BlockVectorType, Constness>::operator=(const Iterator &c)
    {
      parent              = c.parent;
      global_index        = c.global_index;
      index_within_block  = c.index_within_block;
      current_block       = c.current_block;
      next_break_forward  = c.next_break_forward;
      next_break_backward = c.next_break_backward;
 
      return *this;
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline typename Iterator<BlockVectorType, Constness>::dereference_type
    Iterator<BlockVectorType, Constness>::operator*() const
    {
      return parent->block(current_block)(index_within_block);
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline typename Iterator<BlockVectorType, Constness>::dereference_type
    Iterator<BlockVectorType, Constness>::operator[](
      const difference_type d) const
    {
      // if the index pointed to is
      // still within the block we
      // currently point into, then we
      // can save the computation of
      // the block
      if ((global_index + d >= next_break_backward) &&
          (global_index + d <= next_break_forward))
        return parent->block(current_block)(index_within_block + d);
 
      // if the index is not within the
      // block of the block vector into
      // which we presently point, then
      // there is no way: we have to
      // search for the block. this can
      // be done through the parent
      // class as well.
      return (*parent)(global_index + d);
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness> &
    Iterator<BlockVectorType, Constness>::operator++()
    {
      move_forward();
      return *this;
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness>
    Iterator<BlockVectorType, Constness>::operator++(int)
    {
      const Iterator old_value = *this;
      move_forward();
      return old_value;
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness> &
    Iterator<BlockVectorType, Constness>::operator--()
    {
      move_backward();
      return *this;
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness>
    Iterator<BlockVectorType, Constness>::operator--(int)
    {
      const Iterator old_value = *this;
      move_backward();
      return old_value;
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline bool
    Iterator<BlockVectorType, Constness>::operator==(
      const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert(parent == i.parent, ExcPointerToDifferentVectors());
 
      return (global_index == i.global_index);
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline bool
    Iterator<BlockVectorType, Constness>::operator!=(
      const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert(parent == i.parent, ExcPointerToDifferentVectors());
 
      return (global_index != i.global_index);
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline bool
    Iterator<BlockVectorType, Constness>::operator<(
      const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert(parent == i.parent, ExcPointerToDifferentVectors());
 
      return (global_index < i.global_index);
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline bool
    Iterator<BlockVectorType, Constness>::operator<=(
      const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert(parent == i.parent, ExcPointerToDifferentVectors());
 
      return (global_index <= i.global_index);
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline bool
    Iterator<BlockVectorType, Constness>::operator>(
      const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert(parent == i.parent, ExcPointerToDifferentVectors());
 
      return (global_index > i.global_index);
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline bool
    Iterator<BlockVectorType, Constness>::operator>=(
      const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert(parent == i.parent, ExcPointerToDifferentVectors());
 
      return (global_index >= i.global_index);
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline typename Iterator<BlockVectorType, Constness>::difference_type
    Iterator<BlockVectorType, Constness>::operator-(
      const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert(parent == i.parent, ExcPointerToDifferentVectors());
 
      return (static_cast<signed int>(global_index) -
              static_cast<signed int>(i.global_index));
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness>
    Iterator<BlockVectorType, Constness>::operator+(
      const difference_type &d) const
    {
      // if the index pointed to is
      // still within the block we
      // currently point into, then we
      // can save the computation of
      // the block
      if ((global_index + d >= next_break_backward) &&
          (global_index + d <= next_break_forward))
        return Iterator(*parent,
                        global_index + d,
                        current_block,
                        index_within_block + d,
                        next_break_forward,
                        next_break_backward);
      else
        // outside present block, so
        // have to seek new block
        // anyway
        return Iterator(*parent, global_index + d);
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness>
    Iterator<BlockVectorType, Constness>::operator-(
      const difference_type &d) const
    {
      // if the index pointed to is
      // still within the block we
      // currently point into, then we
      // can save the computation of
      // the block
      if ((global_index - d >= next_break_backward) &&
          (global_index - d <= next_break_forward))
        return Iterator(*parent,
                        global_index - d,
                        current_block,
                        index_within_block - d,
                        next_break_forward,
                        next_break_backward);
      else
        // outside present block, so
        // have to seek new block
        // anyway
        return Iterator(*parent, global_index - d);
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness> &
    Iterator<BlockVectorType, Constness>::operator+=(const difference_type &d)
    {
      // if the index pointed to is
      // still within the block we
      // currently point into, then we
      // can save the computation of
      // the block
      if ((global_index + d >= next_break_backward) &&
          (global_index + d <= next_break_forward))
        {
          global_index += d;
          index_within_block += d;
        }
      else
        // outside present block, so
        // have to seek new block
        // anyway
        *this = Iterator(*parent, global_index + d);
 
      return *this;
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    inline Iterator<BlockVectorType, Constness> &
    Iterator<BlockVectorType, Constness>::operator-=(const difference_type &d)
    {
      // if the index pointed to is
      // still within the block we
      // currently point into, then we
      // can save the computation of
      // the block
      if ((global_index - d >= next_break_backward) &&
          (global_index - d <= next_break_forward))
        {
          global_index -= d;
          index_within_block -= d;
        }
      else
        // outside present block, so
        // have to seek new block
        // anyway
        *this = Iterator(*parent, global_index - d);
 
      return *this;
    }
 
 
    template <typename BlockVectorType, bool Constness>
    Iterator<BlockVectorType, Constness>::Iterator(BlockVector    &parent,
                                                   const size_type global_index)
      : parent(&parent)
      , global_index(global_index)
    {
      // find which block we are
      // in. for this, take into
      // account that it happens at
      // times that people want to
      // initialize iterators
      // past-the-end
      if (global_index < parent.size())
        {
          const std::pair<size_type, size_type> indices =
            parent.block_indices.global_to_local(global_index);
          current_block      = indices.first;
          index_within_block = indices.second;
 
          next_break_backward =
            parent.block_indices.local_to_global(current_block, 0);
          next_break_forward =
            (parent.block_indices.local_to_global(current_block, 0) +
             parent.block_indices.block_size(current_block) - 1);
        }
      else
        // past the end. only have one
        // value for this
        {
          this->global_index  = parent.size();
          current_block       = parent.n_blocks();
          index_within_block  = 0;
          next_break_backward = global_index;
          next_break_forward  = numbers::invalid_size_type;
        };
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    void
    Iterator<BlockVectorType, Constness>::move_forward()
    {
      if (global_index != next_break_forward)
        ++index_within_block;
      else
        {
          // ok, we traverse a boundary
          // between blocks:
          index_within_block = 0;
          ++current_block;
 
          // break backwards is now old
          // break forward
          next_break_backward = next_break_forward + 1;
 
          // compute new break forward
          if (current_block < parent->block_indices.size())
            next_break_forward +=
              parent->block_indices.block_size(current_block);
          else
            // if we are beyond the end,
            // then move the next
            // boundary arbitrarily far
            // away
            next_break_forward = numbers::invalid_size_type;
        };
 
      ++global_index;
    }
 
 
 
    template <typename BlockVectorType, bool Constness>
    void
    Iterator<BlockVectorType, Constness>::move_backward()
    {
      if (global_index != next_break_backward)
        --index_within_block;
      else if (current_block != 0)
        {
          // ok, we traverse a boundary
          // between blocks:
          --current_block;
          index_within_block =
            parent->block_indices.block_size(current_block) - 1;
 
          // break forwards is now old
          // break backward
          next_break_forward = next_break_backward - 1;
 
          // compute new break forward
          next_break_backward -=
            parent->block_indices.block_size(current_block);
        }
      else
        // current block was 0, we now
        // get into unspecified terrain
        {
          --current_block;
          index_within_block  = numbers::invalid_size_type;
          next_break_forward  = 0;
          next_break_backward = 0;
        };
 
      --global_index;
    }
 
 
  } // namespace BlockVectorIterators
 
} // namespace internal
 
 
 
template <typename VectorType>
inline typename BlockVectorBase<VectorType>::size_type
BlockVectorBase<VectorType>::size() const
{
  return block_indices.total_size();
}
 
 
 
template <typename VectorType>
inline std::size_t
BlockVectorBase<VectorType>::locally_owned_size() const
{
  std::size_t local_size = 0;
  for (unsigned int b = 0; b < n_blocks(); ++b)
    local_size += block(b).locally_owned_size();
  return local_size;
}
 
 
 
template <typename VectorType>
inline IndexSet
BlockVectorBase<VectorType>::locally_owned_elements() const
{
  IndexSet is(size());
 
  // copy index sets from blocks into the global one, shifted
  // by the appropriate amount for each block
  for (unsigned int b = 0; b < n_blocks(); ++b)
    {
      IndexSet x = block(b).locally_owned_elements();
      is.add_indices(x, block_indices.block_start(b));
    }
 
  is.compress();
 
  return is;
}
 
 
 
template <typename VectorType>
inline unsigned int
BlockVectorBase<VectorType>::n_blocks() const
{
  return block_indices.size();
}
 
 
template <typename VectorType>
inline typename BlockVectorBase<VectorType>::BlockType &
BlockVectorBase<VectorType>::block(const unsigned int i)
{
  AssertIndexRange(i, n_blocks());
 
  return components[i];
}
 
 
 
template <typename VectorType>
inline const typename BlockVectorBase<VectorType>::BlockType &
BlockVectorBase<VectorType>::block(const unsigned int i) const
{
  AssertIndexRange(i, n_blocks());
 
  return components[i];
}
 
 
 
template <typename VectorType>
inline const BlockIndices &
BlockVectorBase<VectorType>::get_block_indices() const
{
  return block_indices;
}
 
 
template <typename VectorType>
inline void
BlockVectorBase<VectorType>::collect_sizes()
{
  std::vector<size_type> sizes(n_blocks());
 
  for (size_type i = 0; i < n_blocks(); ++i)
    sizes[i] = block(i).size();
 
  block_indices.reinit(sizes);
}
 
 
 
template <typename VectorType>
inline void
BlockVectorBase<VectorType>::compress(VectorOperation::values operation)
{
  for (unsigned int i = 0; i < n_blocks(); ++i)
    block(i).compress(operation);
}
 
 
 
template <typename VectorType>
inline typename BlockVectorBase<VectorType>::iterator
BlockVectorBase<VectorType>::begin()
{
  return iterator(*this, 0U);
}
 
 
 
template <typename VectorType>
inline typename BlockVectorBase<VectorType>::const_iterator
BlockVectorBase<VectorType>::begin() const
{
  return const_iterator(*this, 0U);
}
 
 
template <typename VectorType>
inline typename BlockVectorBase<VectorType>::iterator
BlockVectorBase<VectorType>::end()
{
  return iterator(*this, size());
}
 
 
 
template <typename VectorType>
inline typename BlockVectorBase<VectorType>::const_iterator
BlockVectorBase<VectorType>::end() const
{
  return const_iterator(*this, size());
}
 
 
template <typename VectorType>
inline bool
BlockVectorBase<VectorType>::in_local_range(const size_type global_index) const
{
  const std::pair<size_type, size_type> local_index =
    block_indices.global_to_local(global_index);
 
  return components[local_index.first].in_local_range(global_index);
}
 
 
template <typename VectorType>
bool
BlockVectorBase<VectorType>::all_zero() const
{
  for (size_type i = 0; i < n_blocks(); ++i)
    if (components[i].all_zero() == false)
      return false;
 
  return true;
}
 
 
 
template <typename VectorType>
bool
BlockVectorBase<VectorType>::is_non_negative() const
{
  for (size_type i = 0; i < n_blocks(); ++i)
    if (components[i].is_non_negative() == false)
      return false;
 
  return true;
}
 
 
 
template <typename VectorType>
typename BlockVectorBase<VectorType>::value_type
BlockVectorBase<VectorType>::operator*(
  const BlockVectorBase<VectorType> &v) const
{
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
 
  value_type sum = 0.;
  for (size_type i = 0; i < n_blocks(); ++i)
    sum += components[i] * v.components[i];
 
  return sum;
}
 
 
template <typename VectorType>
typename BlockVectorBase<VectorType>::real_type
BlockVectorBase<VectorType>::norm_sqr() const
{
  real_type sum = 0.;
  for (size_type i = 0; i < n_blocks(); ++i)
    sum += components[i].norm_sqr();
 
  return sum;
}
 
 
 
template <typename VectorType>
typename BlockVectorBase<VectorType>::value_type
BlockVectorBase<VectorType>::mean_value() const
{
  value_type sum = 0.;
  // need to do static_cast as otherwise it won't work with
  // value_type=complex<T>
  for (size_type i = 0; i < n_blocks(); ++i)
    sum += components[i].mean_value() *
           (typename numbers::NumberTraits<value_type>::real_type(
             components[i].size()));
 
  return sum / (typename numbers::NumberTraits<value_type>::real_type(size()));
}
 
 
 
template <typename VectorType>
typename BlockVectorBase<VectorType>::real_type
BlockVectorBase<VectorType>::l1_norm() const
{
  real_type sum = 0.;
  for (size_type i = 0; i < n_blocks(); ++i)
    sum += components[i].l1_norm();
 
  return sum;
}
 
 
 
template <typename VectorType>
typename BlockVectorBase<VectorType>::real_type
BlockVectorBase<VectorType>::l2_norm() const
{
  return std::sqrt(norm_sqr());
}
 
 
 
template <typename VectorType>
typename BlockVectorBase<VectorType>::real_type
BlockVectorBase<VectorType>::linfty_norm() const
{
  real_type sum = 0.;
  for (size_type i = 0; i < n_blocks(); ++i)
    {
      value_type newval = components[i].linfty_norm();
      if (sum < newval)
        sum = newval;
    }
  return sum;
}
 
 
 
template <typename VectorType>
typename BlockVectorBase<VectorType>::value_type
BlockVectorBase<VectorType>::add_and_dot(
  const typename BlockVectorBase<VectorType>::value_type a,
  const BlockVectorBase<VectorType>                     &V,
  const BlockVectorBase<VectorType>                     &W)
{
  AssertDimension(n_blocks(), V.n_blocks());
  AssertDimension(n_blocks(), W.n_blocks());
 
  value_type sum = 0.;
  for (size_type i = 0; i < n_blocks(); ++i)
    sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
 
  return sum;
}
 
 
 
template <typename VectorType>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator+=(const BlockVectorBase<VectorType> &v)
{
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
 
  for (size_type i = 0; i < n_blocks(); ++i)
    {
      components[i] += v.components[i];
    }
 
  return *this;
}
 
 
 
template <typename VectorType>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator-=(const BlockVectorBase<VectorType> &v)
{
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
 
  for (size_type i = 0; i < n_blocks(); ++i)
    {
      components[i] -= v.components[i];
    }
  return *this;
}
 
 
 
template <typename VectorType>
template <typename Number>
inline void
BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
                                 const std::vector<Number>    &values)
{
  Assert(indices.size() == values.size(),
         ExcDimensionMismatch(indices.size(), values.size()));
  add(indices.size(), indices.data(), values.data());
}
 
 
 
template <typename VectorType>
template <typename Number>
inline void
BlockVectorBase<VectorType>::add(const std::vector<size_type> &indices,
                                 const Vector<Number>         &values)
{
  Assert(indices.size() == values.size(),
         ExcDimensionMismatch(indices.size(), values.size()));
  const size_type n_indices = indices.size();
  for (size_type i = 0; i < n_indices; ++i)
    (*this)(indices[i]) += values(i);
}
 
 
 
template <typename VectorType>
template <typename Number>
inline void
BlockVectorBase<VectorType>::add(const size_type  n_indices,
                                 const size_type *indices,
                                 const Number    *values)
{
  for (size_type i = 0; i < n_indices; ++i)
    (*this)(indices[i]) += values[i];
}
 
 
 
template <typename VectorType>
void
BlockVectorBase<VectorType>::add(const value_type a)
{
  AssertIsFinite(a);
 
  for (size_type i = 0; i < n_blocks(); ++i)
    {
      components[i].add(a);
    }
}
 
 
 
template <typename VectorType>
void
BlockVectorBase<VectorType>::add(const value_type                   a,
                                 const BlockVectorBase<VectorType> &v)
{
  AssertIsFinite(a);
 
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
 
  for (size_type i = 0; i < n_blocks(); ++i)
    {
      components[i].add(a, v.components[i]);
    }
}
 
 
 
template <typename VectorType>
void
BlockVectorBase<VectorType>::add(const value_type                   a,
                                 const BlockVectorBase<VectorType> &v,
                                 const value_type                   b,
                                 const BlockVectorBase<VectorType> &w)
{
  AssertIsFinite(a);
  AssertIsFinite(b);
 
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
  Assert(n_blocks() == w.n_blocks(),
         ExcDimensionMismatch(n_blocks(), w.n_blocks()));
 
 
  for (size_type i = 0; i < n_blocks(); ++i)
    {
      components[i].add(a, v.components[i], b, w.components[i]);
    }
}
 
 
 
template <typename VectorType>
void
BlockVectorBase<VectorType>::sadd(const value_type                   x,
                                  const BlockVectorBase<VectorType> &v)
{
  AssertIsFinite(x);
 
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
 
  for (size_type i = 0; i < n_blocks(); ++i)
    {
      components[i].sadd(x, v.components[i]);
    }
}
 
 
 
template <typename VectorType>
void
BlockVectorBase<VectorType>::sadd(const value_type                   x,
                                  const value_type                   a,
                                  const BlockVectorBase<VectorType> &v)
{
  AssertIsFinite(x);
  AssertIsFinite(a);
 
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
 
  for (size_type i = 0; i < n_blocks(); ++i)
    {
      components[i].sadd(x, a, v.components[i]);
    }
}
 
 
 
template <typename VectorType>
void
BlockVectorBase<VectorType>::sadd(const value_type                   x,
                                  const value_type                   a,
                                  const BlockVectorBase<VectorType> &v,
                                  const value_type                   b,
                                  const BlockVectorBase<VectorType> &w)
{
  AssertIsFinite(x);
  AssertIsFinite(a);
  AssertIsFinite(b);
 
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
  Assert(n_blocks() == w.n_blocks(),
         ExcDimensionMismatch(n_blocks(), w.n_blocks()));
 
  for (size_type i = 0; i < n_blocks(); ++i)
    {
      components[i].sadd(x, a, v.components[i], b, w.components[i]);
    }
}
 
 
 
template <typename VectorType>
void
BlockVectorBase<VectorType>::sadd(const value_type                   x,
                                  const value_type                   a,
                                  const BlockVectorBase<VectorType> &v,
                                  const value_type                   b,
                                  const BlockVectorBase<VectorType> &w,
                                  const value_type                   c,
                                  const BlockVectorBase<VectorType> &y)
{
  AssertIsFinite(x);
  AssertIsFinite(a);
  AssertIsFinite(b);
  AssertIsFinite(c);
 
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
  Assert(n_blocks() == w.n_blocks(),
         ExcDimensionMismatch(n_blocks(), w.n_blocks()));
  Assert(n_blocks() == y.n_blocks(),
         ExcDimensionMismatch(n_blocks(), y.n_blocks()));
 
  for (size_type i = 0; i < n_blocks(); ++i)
    {
      components[i].sadd(
        x, a, v.components[i], b, w.components[i], c, y.components[i]);
    }
}
 
 
 
template <typename VectorType>
template <class BlockVector2>
void
BlockVectorBase<VectorType>::scale(const BlockVector2 &v)
{
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
  for (size_type i = 0; i < n_blocks(); ++i)
    components[i].scale(v.block(i));
}
 
 
 
template <typename VectorType>
std::size_t
BlockVectorBase<VectorType>::memory_consumption() const
{
  return (MemoryConsumption::memory_consumption(this->block_indices) +
          MemoryConsumption::memory_consumption(this->components));
}
 
 
 
template <typename VectorType>
template <class BlockVector2>
void
BlockVectorBase<VectorType>::equ(const value_type a, const BlockVector2 &v)
{
  AssertIsFinite(a);
 
  Assert(n_blocks() == v.n_blocks(),
         ExcDimensionMismatch(n_blocks(), v.n_blocks()));
 
  for (size_type i = 0; i < n_blocks(); ++i)
    components[i].equ(a, v.components[i]);
}
 
 
 
template <typename VectorType>
void
BlockVectorBase<VectorType>::update_ghost_values() const
{
  for (size_type i = 0; i < n_blocks(); ++i)
    block(i).update_ghost_values();
}
 
 
 
template <typename VectorType>
MPI_Comm
BlockVectorBase<VectorType>::get_mpi_communicator() const
{
  if (n_blocks() > 0)
    return block(0).get_mpi_communicator();
  else
    return MPI_COMM_SELF;
}
 
 
 
template <typename VectorType>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator=(const value_type s)
{
  AssertIsFinite(s);
 
  for (size_type i = 0; i < n_blocks(); ++i)
    components[i] = s;
 
  return *this;
}
 
 
template <typename VectorType>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator=(const BlockVectorBase<VectorType> &v)
{
  AssertDimension(n_blocks(), v.n_blocks());
 
  for (size_type i = 0; i < n_blocks(); ++i)
    components[i] = v.components[i];
 
  return *this;
}
 
 
template <typename VectorType>
template <typename VectorType2>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator=(const BlockVectorBase<VectorType2> &v)
{
  AssertDimension(n_blocks(), v.n_blocks());
 
  for (size_type i = 0; i < n_blocks(); ++i)
    components[i] = v.components[i];
 
  return *this;
}
 
 
 
template <typename VectorType>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator=(const VectorType &v)
{
  Assert(size() == v.size(), ExcDimensionMismatch(size(), v.size()));
 
  size_type index_v = 0;
  for (size_type b = 0; b < n_blocks(); ++b)
    for (size_type i = 0; i < block(b).size(); ++i, ++index_v)
      block(b)(i) = v(index_v);
 
  return *this;
}
 
 
 
template <typename VectorType>
template <typename VectorType2>
inline bool
BlockVectorBase<VectorType>::operator==(
  const BlockVectorBase<VectorType2> &v) const
{
  Assert(block_indices == v.block_indices, ExcDifferentBlockIndices());
 
  for (size_type i = 0; i < n_blocks(); ++i)
    if (!(components[i] == v.components[i]))
      return false;
 
  return true;
}
 
 
 
template <typename VectorType>
inline BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator*=(const value_type factor)
{
  AssertIsFinite(factor);
 
  for (size_type i = 0; i < n_blocks(); ++i)
    components[i] *= factor;
 
  return *this;
}
 
 
 
template <typename VectorType>
inline BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator/=(const value_type factor)
{
  AssertIsFinite(factor);
  Assert(factor != 0., ExcDivideByZero());
 
  const value_type inverse_factor = value_type(1.) / factor;
 
  for (size_type i = 0; i < n_blocks(); ++i)
    components[i] *= inverse_factor;
 
  return *this;
}
 
 
template <typename VectorType>
inline typename BlockVectorBase<VectorType>::value_type
BlockVectorBase<VectorType>::operator()(const size_type i) const
{
  const std::pair<unsigned int, size_type> local_index =
    block_indices.global_to_local(i);
  return components[local_index.first](local_index.second);
}
 
 
 
template <typename VectorType>
inline typename BlockVectorBase<VectorType>::reference
BlockVectorBase<VectorType>::operator()(const size_type i)
{
  const std::pair<unsigned int, size_type> local_index =
    block_indices.global_to_local(i);
  return components[local_index.first](local_index.second);
}
 
 
 
template <typename VectorType>
inline typename BlockVectorBase<VectorType>::value_type
BlockVectorBase<VectorType>::operator[](const size_type i) const
{
  return operator()(i);
}
 
 
 
template <typename VectorType>
inline typename BlockVectorBase<VectorType>::reference
BlockVectorBase<VectorType>::operator[](const size_type i)
{
  return operator()(i);
}
 
 
 
template <typename VectorType>
template <typename OtherNumber>
inline void
BlockVectorBase<VectorType>::extract_subvector_to(
  const std::vector<size_type> &indices,
  std::vector<OtherNumber>     &values) const
{
  for (size_type i = 0; i < indices.size(); ++i)
    values[i] = operator()(indices[i]);
}
 
 
 
template <typename VectorType>
inline void
BlockVectorBase<VectorType>::extract_subvector_to(
  const ArrayView<const types::global_dof_index> &indices,
  const ArrayView<value_type>                    &entries) const
{
  AssertDimension(indices.size(), entries.size());
  for (unsigned int i = 0; i < indices.size(); ++i)
    {
      AssertIndexRange(indices[i], size());
      entries[i] = (*this)[indices[i]];
    }
}
 
 
 
template <typename VectorType>
template <typename ForwardIterator, typename OutputIterator>
inline void
BlockVectorBase<VectorType>::extract_subvector_to(
  ForwardIterator       indices_begin,
  const ForwardIterator indices_end,
  OutputIterator        values_begin) const
{
  while (indices_begin != indices_end)
    {
      *values_begin = operator()(*indices_begin);
      ++indices_begin;
      ++values_begin;
    }
}
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif