block_matrix_base.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2020 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_block_matrix_base_h
#define dealii_block_matrix_base_h


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

#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/table.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/utilities.h>

#include <deal.II/lac/block_indices.h>
#include <deal.II/lac/exceptions.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/matrix_iterator.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_operation.h>

#include <cmath>

DEAL_II_NAMESPACE_OPEN


// Forward declaration
#ifndef DOXYGEN
template <typename>
class MatrixIterator;
#endif


namespace BlockMatrixIterators
{
  template <class BlockMatrixType>
  class AccessorBase
  {
  public:
    using size_type = types::global_dof_index;

    using value_type = typename BlockMatrixType::value_type;

    AccessorBase();

    unsigned int
    block_row() const;

    unsigned int
    block_column() const;

  protected:
    unsigned int row_block;

    unsigned int col_block;

    // Let the iterator class be a friend.
    template <typename>
    friend class MatrixIterator;
  };



  template <class BlockMatrixType, bool Constness>
  class Accessor;


  template <class BlockMatrixType>
  class Accessor<BlockMatrixType, false> : public AccessorBase<BlockMatrixType>
  {
  public:
    using size_type = types::global_dof_index;

    using MatrixType = BlockMatrixType;

    using value_type = typename BlockMatrixType::value_type;

    Accessor(BlockMatrixType *m, const size_type row, const size_type col);

    size_type
    row() const;

    size_type
    column() const;

    value_type
    value() const;

    void
    set_value(value_type newval) const;

  protected:
    BlockMatrixType *matrix;

    typename BlockMatrixType::BlockType::iterator base_iterator;

    void
    advance();

    bool
    operator==(const Accessor &a) const;

    template <typename>
    friend class MatrixIterator;
    friend class Accessor<BlockMatrixType, true>;
  };

  template <class BlockMatrixType>
  class Accessor<BlockMatrixType, true> : public AccessorBase<BlockMatrixType>
  {
  public:
    using size_type = types::global_dof_index;

    using MatrixType = const BlockMatrixType;

    using value_type = typename BlockMatrixType::value_type;

    Accessor(const BlockMatrixType *m,
             const size_type        row,
             const size_type        col);

    Accessor(const Accessor<BlockMatrixType, false> &);

    size_type
    row() const;

    size_type
    column() const;

    value_type
    value() const;

  protected:
    const BlockMatrixType *matrix;

    typename BlockMatrixType::BlockType::const_iterator base_iterator;

    void
    advance();

    bool
    operator==(const Accessor &a) const;

    // Let the iterator class be a friend.
    template <typename>
    friend class ::MatrixIterator;
  };
} // namespace BlockMatrixIterators



template <typename MatrixType>
class BlockMatrixBase : public Subscriptor
{
public:
  using BlockType = MatrixType;

  using value_type      = typename BlockType::value_type;
  using real_type       = typename numbers::NumberTraits<value_type>::real_type;
  using pointer         = value_type *;
  using const_pointer   = const value_type *;
  using reference       = value_type &;
  using const_reference = const value_type &;
  using size_type       = types::global_dof_index;

  using iterator =
    MatrixIterator<BlockMatrixIterators::Accessor<BlockMatrixBase, false>>;

  using const_iterator =
    MatrixIterator<BlockMatrixIterators::Accessor<BlockMatrixBase, true>>;


  BlockMatrixBase() = default;

  ~BlockMatrixBase() override;

  template <class BlockMatrixType>
  BlockMatrixBase &
  copy_from(const BlockMatrixType &source);

  BlockType &
  block(const unsigned int row, const unsigned int column);


  const BlockType &
  block(const unsigned int row, const unsigned int column) const;

  size_type
  m() const;

  size_type
  n() const;


  unsigned int
  n_block_rows() const;

  unsigned int
  n_block_cols() const;

  void
  set(const size_type i, const size_type j, const value_type value);

  template <typename number>
  void
  set(const std::vector<size_type> &indices,
      const FullMatrix<number> &    full_matrix,
      const bool                    elide_zero_values = false);

  template <typename number>
  void
  set(const std::vector<size_type> &row_indices,
      const std::vector<size_type> &col_indices,
      const FullMatrix<number> &    full_matrix,
      const bool                    elide_zero_values = false);

  template <typename number>
  void
  set(const size_type               row,
      const std::vector<size_type> &col_indices,
      const std::vector<number> &   values,
      const bool                    elide_zero_values = false);

  template <typename number>
  void
  set(const size_type  row,
      const size_type  n_cols,
      const size_type *col_indices,
      const number *   values,
      const bool       elide_zero_values = false);

  void
  add(const size_type i, const size_type j, const value_type value);

  template <typename number>
  void
  add(const std::vector<size_type> &indices,
      const FullMatrix<number> &    full_matrix,
      const bool                    elide_zero_values = true);

  template <typename number>
  void
  add(const std::vector<size_type> &row_indices,
      const std::vector<size_type> &col_indices,
      const FullMatrix<number> &    full_matrix,
      const bool                    elide_zero_values = true);

  template <typename number>
  void
  add(const size_type               row,
      const std::vector<size_type> &col_indices,
      const std::vector<number> &   values,
      const bool                    elide_zero_values = true);

  template <typename number>
  void
  add(const size_type  row,
      const size_type  n_cols,
      const size_type *col_indices,
      const number *   values,
      const bool       elide_zero_values      = true,
      const bool       col_indices_are_sorted = false);

  void
  add(const value_type factor, const BlockMatrixBase<MatrixType> &matrix);

  value_type
  operator()(const size_type i, const size_type j) const;

  value_type
  el(const size_type i, const size_type j) const;

  value_type
  diag_element(const size_type i) const;

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

  BlockMatrixBase &
  operator*=(const value_type factor);

  BlockMatrixBase &
  operator/=(const value_type factor);

  template <class BlockVectorType>
  void
  vmult_add(BlockVectorType &dst, const BlockVectorType &src) const;

  template <class BlockVectorType>
  void
  Tvmult_add(BlockVectorType &dst, const BlockVectorType &src) const;

  template <class BlockVectorType>
  value_type
  matrix_norm_square(const BlockVectorType &v) const;

  real_type
  frobenius_norm() const;

  template <class BlockVectorType>
  value_type
  matrix_scalar_product(const BlockVectorType &u,
                        const BlockVectorType &v) const;

  template <class BlockVectorType>
  value_type
  residual(BlockVectorType &      dst,
           const BlockVectorType &x,
           const BlockVectorType &b) const;

  void
  print(std::ostream &out, const bool alternative_output = false) const;

  iterator
  begin();

  iterator
  end();

  iterator
  begin(const size_type r);

  iterator
  end(const size_type r);
  const_iterator
  begin() const;

  const_iterator
  end() const;

  const_iterator
  begin(const size_type r) const;

  const_iterator
  end(const size_type r) const;

  const BlockIndices &
  get_row_indices() const;

  const BlockIndices &
  get_column_indices() const;

  std::size_t
  memory_consumption() const;

  DeclException4(ExcIncompatibleRowNumbers,
                 int,
                 int,
                 int,
                 int,
                 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
                 << ',' << arg4 << "] have differing row numbers.");
  DeclException4(ExcIncompatibleColNumbers,
                 int,
                 int,
                 int,
                 int,
                 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
                 << ',' << arg4 << "] have differing column numbers.");
protected:
  void
  clear();

  BlockIndices row_block_indices;
  BlockIndices column_block_indices;

  Table<2, SmartPointer<BlockType, BlockMatrixBase<MatrixType>>> sub_objects;

  void
  collect_sizes();

  template <class BlockVectorType>
  void
  vmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;

  template <class BlockVectorType, class VectorType>
  void
  vmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;

  template <class BlockVectorType, class VectorType>
  void
  vmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;

  template <class VectorType>
  void
  vmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;

  template <class BlockVectorType>
  void
  Tvmult_block_block(BlockVectorType &dst, const BlockVectorType &src) const;

  template <class BlockVectorType, class VectorType>
  void
  Tvmult_block_nonblock(BlockVectorType &dst, const VectorType &src) const;

  template <class BlockVectorType, class VectorType>
  void
  Tvmult_nonblock_block(VectorType &dst, const BlockVectorType &src) const;

  template <class VectorType>
  void
  Tvmult_nonblock_nonblock(VectorType &dst, const VectorType &src) const;


protected:
  void
  prepare_add_operation();

  void
  prepare_set_operation();


private:
  struct TemporaryData
  {
    std::vector<size_type> counter_within_block;

    std::vector<std::vector<size_type>> column_indices;

    std::vector<std::vector<value_type>> column_values;

    std::mutex mutex;

    TemporaryData &
    operator=(const TemporaryData &)
    {
      return *this;
    }
  };

  TemporaryData temporary_data;

  // Make the iterator class a friend. We have to work around a compiler bug
  // here again.
  template <typename, bool>
  friend class BlockMatrixIterators::Accessor;

  template <typename>
  friend class MatrixIterator;
};


#ifndef DOXYGEN
/* ------------------------- Template functions ---------------------- */


namespace BlockMatrixIterators
{
  template <class BlockMatrixType>
  inline AccessorBase<BlockMatrixType>::AccessorBase()
    : row_block(0)
    , col_block(0)
  {}


  template <class BlockMatrixType>
  inline unsigned int
  AccessorBase<BlockMatrixType>::block_row() const
  {
    Assert(row_block != numbers::invalid_unsigned_int, ExcIteratorPastEnd());

    return row_block;
  }


  template <class BlockMatrixType>
  inline unsigned int
  AccessorBase<BlockMatrixType>::block_column() const
  {
    Assert(col_block != numbers::invalid_unsigned_int, ExcIteratorPastEnd());

    return col_block;
  }


  template <class BlockMatrixType>
  inline Accessor<BlockMatrixType, true>::Accessor(
    const BlockMatrixType *matrix,
    const size_type        row,
    const size_type        col)
    : matrix(matrix)
    , base_iterator(matrix->block(0, 0).begin())
  {
    (void)col;
    Assert(col == 0, ExcNotImplemented());

    // check if this is a regular row or
    // the end of the matrix
    if (row < matrix->m())
      {
        const std::pair<unsigned int, size_type> indices =
          matrix->row_block_indices.global_to_local(row);

        // find the first block that does
        // have an entry in this row
        for (unsigned int bc = 0; bc < matrix->n_block_cols(); ++bc)
          {
            base_iterator =
              matrix->block(indices.first, bc).begin(indices.second);
            if (base_iterator !=
                matrix->block(indices.first, bc).end(indices.second))
              {
                this->row_block = indices.first;
                this->col_block = bc;
                return;
              }
          }

        // hm, there is no block that has
        // an entry in this column. we need
        // to take the next entry then,
        // which may be the first entry of
        // the next row, or recursively the
        // next row, or so on
        *this = Accessor(matrix, row + 1, 0);
      }
    else
      {
        // we were asked to create the end
        // iterator for this matrix
        this->row_block = numbers::invalid_unsigned_int;
        this->col_block = numbers::invalid_unsigned_int;
      }
  }


  //   template <class BlockMatrixType>
  //   inline
  //   Accessor<BlockMatrixType, true>::Accessor (const
  //   Accessor<BlockMatrixType, true>& other)
  //                :
  //                matrix(other.matrix),
  //                base_iterator(other.base_iterator)
  //   {
  //     this->row_block = other.row_block;
  //     this->col_block = other.col_block;
  //   }


  template <class BlockMatrixType>
  inline Accessor<BlockMatrixType, true>::Accessor(
    const Accessor<BlockMatrixType, false> &other)
    : matrix(other.matrix)
    , base_iterator(other.base_iterator)
  {
    this->row_block = other.row_block;
    this->col_block = other.col_block;
  }


  template <class BlockMatrixType>
  inline typename Accessor<BlockMatrixType, true>::size_type
  Accessor<BlockMatrixType, true>::row() const
  {
    Assert(this->row_block != numbers::invalid_unsigned_int,
           ExcIteratorPastEnd());

    return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
            base_iterator->row());
  }


  template <class BlockMatrixType>
  inline typename Accessor<BlockMatrixType, true>::size_type
  Accessor<BlockMatrixType, true>::column() const
  {
    Assert(this->col_block != numbers::invalid_unsigned_int,
           ExcIteratorPastEnd());

    return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
            base_iterator->column());
  }


  template <class BlockMatrixType>
  inline typename Accessor<BlockMatrixType, true>::value_type
  Accessor<BlockMatrixType, true>::value() const
  {
    Assert(this->row_block != numbers::invalid_unsigned_int,
           ExcIteratorPastEnd());
    Assert(this->col_block != numbers::invalid_unsigned_int,
           ExcIteratorPastEnd());

    return base_iterator->value();
  }



  template <class BlockMatrixType>
  inline void
  Accessor<BlockMatrixType, true>::advance()
  {
    Assert(this->row_block != numbers::invalid_unsigned_int,
           ExcIteratorPastEnd());
    Assert(this->col_block != numbers::invalid_unsigned_int,
           ExcIteratorPastEnd());

    // Remember current row inside block
    size_type local_row = base_iterator->row();

    // Advance one element inside the
    // current block
    ++base_iterator;

    // while we hit the end of the row of a
    // block (which may happen multiple
    // times if rows inside a block are
    // empty), we have to jump to the next
    // block and take the
    while (base_iterator ==
           matrix->block(this->row_block, this->col_block).end(local_row))
      {
        // jump to next block in this block
        // row, if possible, otherwise go
        // to next row
        if (this->col_block < matrix->n_block_cols() - 1)
          {
            ++this->col_block;
            base_iterator =
              matrix->block(this->row_block, this->col_block).begin(local_row);
          }
        else
          {
            // jump back to next row in
            // first block column
            this->col_block = 0;
            ++local_row;

            // see if this has brought us
            // past the number of rows in
            // this block. if so see
            // whether we've just fallen
            // off the end of the whole
            // matrix
            if (local_row ==
                matrix->block(this->row_block, this->col_block).m())
              {
                local_row = 0;
                ++this->row_block;
                if (this->row_block == matrix->n_block_rows())
                  {
                    this->row_block = numbers::invalid_unsigned_int;
                    this->col_block = numbers::invalid_unsigned_int;
                    return;
                  }
              }

            base_iterator =
              matrix->block(this->row_block, this->col_block).begin(local_row);
          }
      }
  }


  template <class BlockMatrixType>
  inline bool
  Accessor<BlockMatrixType, true>::operator==(const Accessor &a) const
  {
    if (matrix != a.matrix)
      return false;

    if (this->row_block == a.row_block && this->col_block == a.col_block)
      // end iterators do not necessarily
      // have to have the same
      // base_iterator representation, but
      // valid iterators have to
      return (((this->row_block == numbers::invalid_unsigned_int) &&
               (this->col_block == numbers::invalid_unsigned_int)) ||
              (base_iterator == a.base_iterator));

    return false;
  }

  //----------------------------------------------------------------------//


  template <class BlockMatrixType>
  inline Accessor<BlockMatrixType, false>::Accessor(BlockMatrixType *matrix,
                                                    const size_type  row,
                                                    const size_type  col)
    : matrix(matrix)
    , base_iterator(matrix->block(0, 0).begin())
  {
    (void)col;
    Assert(col == 0, ExcNotImplemented());
    // check if this is a regular row or
    // the end of the matrix
    if (row < matrix->m())
      {
        const std::pair<unsigned int, size_type> indices =
          matrix->row_block_indices.global_to_local(row);

        // find the first block that does
        // have an entry in this row
        for (size_type bc = 0; bc < matrix->n_block_cols(); ++bc)
          {
            base_iterator =
              matrix->block(indices.first, bc).begin(indices.second);
            if (base_iterator !=
                matrix->block(indices.first, bc).end(indices.second))
              {
                this->row_block = indices.first;
                this->col_block = bc;
                return;
              }
          }

        // hm, there is no block that has
        // an entry in this column. we need
        // to take the next entry then,
        // which may be the first entry of
        // the next row, or recursively the
        // next row, or so on
        *this = Accessor(matrix, row + 1, 0);
      }
    else
      {
        // we were asked to create the end
        // iterator for this matrix
        this->row_block = numbers::invalid_size_type;
        this->col_block = numbers::invalid_size_type;
      }
  }


  template <class BlockMatrixType>
  inline typename Accessor<BlockMatrixType, false>::size_type
  Accessor<BlockMatrixType, false>::row() const
  {
    Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());

    return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
            base_iterator->row());
  }


  template <class BlockMatrixType>
  inline typename Accessor<BlockMatrixType, false>::size_type
  Accessor<BlockMatrixType, false>::column() const
  {
    Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());

    return (matrix->column_block_indices.local_to_global(this->col_block, 0) +
            base_iterator->column());
  }


  template <class BlockMatrixType>
  inline typename Accessor<BlockMatrixType, false>::value_type
  Accessor<BlockMatrixType, false>::value() const
  {
    Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
    Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());

    return base_iterator->value();
  }



  template <class BlockMatrixType>
  inline void
  Accessor<BlockMatrixType, false>::set_value(
    typename Accessor<BlockMatrixType, false>::value_type newval) const
  {
    Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
    Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());

    base_iterator->value() = newval;
  }



  template <class BlockMatrixType>
  inline void
  Accessor<BlockMatrixType, false>::advance()
  {
    Assert(this->row_block != numbers::invalid_size_type, ExcIteratorPastEnd());
    Assert(this->col_block != numbers::invalid_size_type, ExcIteratorPastEnd());

    // Remember current row inside block
    size_type local_row = base_iterator->row();

    // Advance one element inside the
    // current block
    ++base_iterator;

    // while we hit the end of the row of a
    // block (which may happen multiple
    // times if rows inside a block are
    // empty), we have to jump to the next
    // block and take the
    while (base_iterator ==
           matrix->block(this->row_block, this->col_block).end(local_row))
      {
        // jump to next block in this block
        // row, if possible, otherwise go
        // to next row
        if (this->col_block < matrix->n_block_cols() - 1)
          {
            ++this->col_block;
            base_iterator =
              matrix->block(this->row_block, this->col_block).begin(local_row);
          }
        else
          {
            // jump back to next row in
            // first block column
            this->col_block = 0;
            ++local_row;

            // see if this has brought us
            // past the number of rows in
            // this block. if so see
            // whether we've just fallen
            // off the end of the whole
            // matrix
            if (local_row ==
                matrix->block(this->row_block, this->col_block).m())
              {
                local_row = 0;
                ++this->row_block;
                if (this->row_block == matrix->n_block_rows())
                  {
                    this->row_block = numbers::invalid_size_type;
                    this->col_block = numbers::invalid_size_type;
                    return;
                  }
              }

            base_iterator =
              matrix->block(this->row_block, this->col_block).begin(local_row);
          }
      }
  }



  template <class BlockMatrixType>
  inline bool
  Accessor<BlockMatrixType, false>::operator==(const Accessor &a) const
  {
    if (matrix != a.matrix)
      return false;

    if (this->row_block == a.row_block && this->col_block == a.col_block)
      // end iterators do not necessarily
      // have to have the same
      // base_iterator representation, but
      // valid iterators have to
      return (((this->row_block == numbers::invalid_size_type) &&
               (this->col_block == numbers::invalid_size_type)) ||
              (base_iterator == a.base_iterator));

    return false;
  }
} // namespace BlockMatrixIterators


//---------------------------------------------------------------------------

template <typename MatrixType>
inline BlockMatrixBase<MatrixType>::~BlockMatrixBase()
{
  try
    {
      clear();
    }
  catch (...)
    {}
}


template <class MatrixType>
template <class BlockMatrixType>
inline BlockMatrixBase<MatrixType> &
BlockMatrixBase<MatrixType>::copy_from(const BlockMatrixType &source)
{
  for (unsigned int r = 0; r < n_block_rows(); ++r)
    for (unsigned int c = 0; c < n_block_cols(); ++c)
      block(r, c).copy_from(source.block(r, c));

  return *this;
}


template <class MatrixType>
std::size_t
BlockMatrixBase<MatrixType>::memory_consumption() const
{
  std::size_t mem =
    MemoryConsumption::memory_consumption(row_block_indices) +
    MemoryConsumption::memory_consumption(column_block_indices) +
    MemoryConsumption::memory_consumption(sub_objects) +
    MemoryConsumption::memory_consumption(temporary_data.counter_within_block) +
    MemoryConsumption::memory_consumption(temporary_data.column_indices) +
    MemoryConsumption::memory_consumption(temporary_data.column_values) +
    sizeof(temporary_data.mutex);

  for (unsigned int r = 0; r < n_block_rows(); ++r)
    for (unsigned int c = 0; c < n_block_cols(); ++c)
      {
        MatrixType *p = this->sub_objects[r][c];
        mem += MemoryConsumption::memory_consumption(*p);
      }

  return mem;
}



template <class MatrixType>
inline void
BlockMatrixBase<MatrixType>::clear()
{
  for (unsigned int r = 0; r < n_block_rows(); ++r)
    for (unsigned int c = 0; c < n_block_cols(); ++c)
      {
        MatrixType *p           = this->sub_objects[r][c];
        this->sub_objects[r][c] = nullptr;
        delete p;
      }
  sub_objects.reinit(0, 0);

  // reset block indices to empty
  row_block_indices = column_block_indices = BlockIndices();
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::BlockType &
BlockMatrixBase<MatrixType>::block(const unsigned int row,
                                   const unsigned int column)
{
  AssertIndexRange(row, n_block_rows());
  AssertIndexRange(column, n_block_cols());

  return *sub_objects[row][column];
}



template <class MatrixType>
inline const typename BlockMatrixBase<MatrixType>::BlockType &
BlockMatrixBase<MatrixType>::block(const unsigned int row,
                                   const unsigned int column) const
{
  AssertIndexRange(row, n_block_rows());
  AssertIndexRange(column, n_block_cols());

  return *sub_objects[row][column];
}


template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::size_type
BlockMatrixBase<MatrixType>::m() const
{
  return row_block_indices.total_size();
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::size_type
BlockMatrixBase<MatrixType>::n() const
{
  return column_block_indices.total_size();
}



template <class MatrixType>
inline unsigned int
BlockMatrixBase<MatrixType>::n_block_cols() const
{
  return column_block_indices.size();
}



template <class MatrixType>
inline unsigned int
BlockMatrixBase<MatrixType>::n_block_rows() const
{
  return row_block_indices.size();
}



// Write the single set manually,
// since the other function has a lot
// of overhead in that case.
template <class MatrixType>
inline void
BlockMatrixBase<MatrixType>::set(const size_type  i,
                                 const size_type  j,
                                 const value_type value)
{
  prepare_set_operation();

  AssertIsFinite(value);

  const std::pair<unsigned int, size_type>
    row_index = row_block_indices.global_to_local(i),
    col_index = column_block_indices.global_to_local(j);
  block(row_index.first, col_index.first)
    .set(row_index.second, col_index.second, value);
}



template <class MatrixType>
template <typename number>
inline void
BlockMatrixBase<MatrixType>::set(const std::vector<size_type> &row_indices,
                                 const std::vector<size_type> &col_indices,
                                 const FullMatrix<number> &    values,
                                 const bool elide_zero_values)
{
  Assert(row_indices.size() == values.m(),
         ExcDimensionMismatch(row_indices.size(), values.m()));
  Assert(col_indices.size() == values.n(),
         ExcDimensionMismatch(col_indices.size(), values.n()));

  for (size_type i = 0; i < row_indices.size(); ++i)
    set(row_indices[i],
        col_indices.size(),
        col_indices.data(),
        &values(i, 0),
        elide_zero_values);
}



template <class MatrixType>
template <typename number>
inline void
BlockMatrixBase<MatrixType>::set(const std::vector<size_type> &indices,
                                 const FullMatrix<number> &    values,
                                 const bool elide_zero_values)
{
  Assert(indices.size() == values.m(),
         ExcDimensionMismatch(indices.size(), values.m()));
  Assert(values.n() == values.m(), ExcNotQuadratic());

  for (size_type i = 0; i < indices.size(); ++i)
    set(indices[i],
        indices.size(),
        indices.data(),
        &values(i, 0),
        elide_zero_values);
}



template <class MatrixType>
template <typename number>
inline void
BlockMatrixBase<MatrixType>::set(const size_type               row,
                                 const std::vector<size_type> &col_indices,
                                 const std::vector<number> &   values,
                                 const bool elide_zero_values)
{
  Assert(col_indices.size() == values.size(),
         ExcDimensionMismatch(col_indices.size(), values.size()));

  set(row,
      col_indices.size(),
      col_indices.data(),
      values.data(),
      elide_zero_values);
}



// This is a very messy function, since
// we need to calculate to each position
// the location in the global array.
template <class MatrixType>
template <typename number>
inline void
BlockMatrixBase<MatrixType>::set(const size_type  row,
                                 const size_type  n_cols,
                                 const size_type *col_indices,
                                 const number *   values,
                                 const bool       elide_zero_values)
{
  prepare_set_operation();

  // lock access to the temporary data structure to
  // allow multiple threads to call this function concurrently
  std::lock_guard<std::mutex> lock(temporary_data.mutex);

  // Resize scratch arrays
  if (temporary_data.column_indices.size() < this->n_block_cols())
    {
      temporary_data.column_indices.resize(this->n_block_cols());
      temporary_data.column_values.resize(this->n_block_cols());
      temporary_data.counter_within_block.resize(this->n_block_cols());
    }

  // Resize sub-arrays to n_cols. This
  // is a bit wasteful, but we resize
  // only a few times (then the maximum
  // row length won't increase that
  // much any more). At least we know
  // that all arrays are going to be of
  // the same size, so we can check
  // whether the size of one is large
  // enough before actually going
  // through all of them.
  if (temporary_data.column_indices[0].size() < n_cols)
    {
      for (unsigned int i = 0; i < this->n_block_cols(); ++i)
        {
          temporary_data.column_indices[i].resize(n_cols);
          temporary_data.column_values[i].resize(n_cols);
        }
    }

  // Reset the number of added elements
  // in each block to zero.
  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
    temporary_data.counter_within_block[i] = 0;

  // Go through the column indices to
  // find out which portions of the
  // values should be set in which
  // block of the matrix. We need to
  // touch all the data, since we can't
  // be sure that the data of one block
  // is stored contiguously (in fact,
  // indices will be intermixed when it
  // comes from an element matrix).
  for (size_type j = 0; j < n_cols; ++j)
    {
      number value = values[j];

      if (value == number() && elide_zero_values == true)
        continue;

      const std::pair<unsigned int, size_type> col_index =
        this->column_block_indices.global_to_local(col_indices[j]);

      const size_type local_index =
        temporary_data.counter_within_block[col_index.first]++;

      temporary_data.column_indices[col_index.first][local_index] =
        col_index.second;
      temporary_data.column_values[col_index.first][local_index] = value;
    }

#  ifdef DEBUG
  // If in debug mode, do a check whether
  // the right length has been obtained.
  size_type length = 0;
  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
    length += temporary_data.counter_within_block[i];
  Assert(length <= n_cols, ExcInternalError());
#  endif

  // Now we found out about where the
  // individual columns should start and
  // where we should start reading out
  // data. Now let's write the data into
  // the individual blocks!
  const std::pair<unsigned int, size_type> row_index =
    this->row_block_indices.global_to_local(row);
  for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
    {
      if (temporary_data.counter_within_block[block_col] == 0)
        continue;

      block(row_index.first, block_col)
        .set(row_index.second,
             temporary_data.counter_within_block[block_col],
             temporary_data.column_indices[block_col].data(),
             temporary_data.column_values[block_col].data(),
             false);
    }
}



template <class MatrixType>
inline void
BlockMatrixBase<MatrixType>::add(const size_type  i,
                                 const size_type  j,
                                 const value_type value)
{
  AssertIsFinite(value);

  prepare_add_operation();

  // save some cycles for zero additions, but
  // only if it is safe for the matrix we are
  // working with
  using MatrixTraits = typename MatrixType::Traits;
  if ((MatrixTraits::zero_addition_can_be_elided == true) &&
      (value == value_type()))
    return;

  const std::pair<unsigned int, size_type>
    row_index = row_block_indices.global_to_local(i),
    col_index = column_block_indices.global_to_local(j);
  block(row_index.first, col_index.first)
    .add(row_index.second, col_index.second, value);
}



template <class MatrixType>
template <typename number>
inline void
BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &row_indices,
                                 const std::vector<size_type> &col_indices,
                                 const FullMatrix<number> &    values,
                                 const bool elide_zero_values)
{
  Assert(row_indices.size() == values.m(),
         ExcDimensionMismatch(row_indices.size(), values.m()));
  Assert(col_indices.size() == values.n(),
         ExcDimensionMismatch(col_indices.size(), values.n()));

  for (size_type i = 0; i < row_indices.size(); ++i)
    add(row_indices[i],
        col_indices.size(),
        col_indices.data(),
        &values(i, 0),
        elide_zero_values);
}



template <class MatrixType>
template <typename number>
inline void
BlockMatrixBase<MatrixType>::add(const std::vector<size_type> &indices,
                                 const FullMatrix<number> &    values,
                                 const bool elide_zero_values)
{
  Assert(indices.size() == values.m(),
         ExcDimensionMismatch(indices.size(), values.m()));
  Assert(values.n() == values.m(), ExcNotQuadratic());

  for (size_type i = 0; i < indices.size(); ++i)
    add(indices[i],
        indices.size(),
        indices.data(),
        &values(i, 0),
        elide_zero_values);
}



template <class MatrixType>
template <typename number>
inline void
BlockMatrixBase<MatrixType>::add(const size_type               row,
                                 const std::vector<size_type> &col_indices,
                                 const std::vector<number> &   values,
                                 const bool elide_zero_values)
{
  Assert(col_indices.size() == values.size(),
         ExcDimensionMismatch(col_indices.size(), values.size()));

  add(row,
      col_indices.size(),
      col_indices.data(),
      values.data(),
      elide_zero_values);
}



// This is a very messy function, since
// we need to calculate to each position
// the location in the global array.
template <class MatrixType>
template <typename number>
inline void
BlockMatrixBase<MatrixType>::add(const size_type  row,
                                 const size_type  n_cols,
                                 const size_type *col_indices,
                                 const number *   values,
                                 const bool       elide_zero_values,
                                 const bool       col_indices_are_sorted)
{
  prepare_add_operation();

  // TODO: Look over this to find out
  // whether we can do that more
  // efficiently.
  if (col_indices_are_sorted == true)
    {
#  ifdef DEBUG
      // check whether indices really are
      // sorted.
      size_type before = col_indices[0];
      for (size_type i = 1; i < n_cols; ++i)
        if (col_indices[i] <= before)
          Assert(false,
                 ExcMessage("Flag col_indices_are_sorted is set, but "
                            "indices appear to not be sorted.")) else before =
            col_indices[i];
#  endif
      const std::pair<unsigned int, size_type> row_index =
        this->row_block_indices.global_to_local(row);

      if (this->n_block_cols() > 1)
        {
          const size_type *first_block =
            Utilities::lower_bound(col_indices,
                                   col_indices + n_cols,
                                   this->column_block_indices.block_start(1));

          const size_type n_zero_block_indices = first_block - col_indices;
          block(row_index.first, 0)
            .add(row_index.second,
                 n_zero_block_indices,
                 col_indices,
                 values,
                 elide_zero_values,
                 col_indices_are_sorted);

          if (n_zero_block_indices < n_cols)
            this->add(row,
                      n_cols - n_zero_block_indices,
                      first_block,
                      values + n_zero_block_indices,
                      elide_zero_values,
                      false);
        }
      else
        {
          block(row_index.first, 0)
            .add(row_index.second,
                 n_cols,
                 col_indices,
                 values,
                 elide_zero_values,
                 col_indices_are_sorted);
        }

      return;
    }

  // Lock scratch arrays, then resize them
  std::lock_guard<std::mutex> lock(temporary_data.mutex);

  if (temporary_data.column_indices.size() < this->n_block_cols())
    {
      temporary_data.column_indices.resize(this->n_block_cols());
      temporary_data.column_values.resize(this->n_block_cols());
      temporary_data.counter_within_block.resize(this->n_block_cols());
    }

  // Resize sub-arrays to n_cols. This
  // is a bit wasteful, but we resize
  // only a few times (then the maximum
  // row length won't increase that
  // much any more). At least we know
  // that all arrays are going to be of
  // the same size, so we can check
  // whether the size of one is large
  // enough before actually going
  // through all of them.
  if (temporary_data.column_indices[0].size() < n_cols)
    {
      for (unsigned int i = 0; i < this->n_block_cols(); ++i)
        {
          temporary_data.column_indices[i].resize(n_cols);
          temporary_data.column_values[i].resize(n_cols);
        }
    }

  // Reset the number of added elements
  // in each block to zero.
  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
    temporary_data.counter_within_block[i] = 0;

  // Go through the column indices to
  // find out which portions of the
  // values should be written into
  // which block of the matrix. We need
  // to touch all the data, since we
  // can't be sure that the data of one
  // block is stored contiguously (in
  // fact, data will be intermixed when
  // it comes from an element matrix).
  for (size_type j = 0; j < n_cols; ++j)
    {
      number value = values[j];

      if (value == number() && elide_zero_values == true)
        continue;

      const std::pair<unsigned int, size_type> col_index =
        this->column_block_indices.global_to_local(col_indices[j]);

      const size_type local_index =
        temporary_data.counter_within_block[col_index.first]++;

      temporary_data.column_indices[col_index.first][local_index] =
        col_index.second;
      temporary_data.column_values[col_index.first][local_index] = value;
    }

#  ifdef DEBUG
  // If in debug mode, do a check whether
  // the right length has been obtained.
  size_type length = 0;
  for (unsigned int i = 0; i < this->n_block_cols(); ++i)
    length += temporary_data.counter_within_block[i];
  Assert(length <= n_cols, ExcInternalError());
#  endif

  // Now we found out about where the
  // individual columns should start and
  // where we should start reading out
  // data. Now let's write the data into
  // the individual blocks!
  const std::pair<unsigned int, size_type> row_index =
    this->row_block_indices.global_to_local(row);
  for (unsigned int block_col = 0; block_col < n_block_cols(); ++block_col)
    {
      if (temporary_data.counter_within_block[block_col] == 0)
        continue;

      block(row_index.first, block_col)
        .add(row_index.second,
             temporary_data.counter_within_block[block_col],
             temporary_data.column_indices[block_col].data(),
             temporary_data.column_values[block_col].data(),
             false,
             col_indices_are_sorted);
    }
}



template <class MatrixType>
inline void
BlockMatrixBase<MatrixType>::add(const value_type                   factor,
                                 const BlockMatrixBase<MatrixType> &matrix)
{
  AssertIsFinite(factor);

  prepare_add_operation();

  // save some cycles for zero additions, but
  // only if it is safe for the matrix we are
  // working with
  using MatrixTraits = typename MatrixType::Traits;
  if ((MatrixTraits::zero_addition_can_be_elided == true) && (factor == 0))
    return;

  for (unsigned int row = 0; row < n_block_rows(); ++row)
    for (unsigned int col = 0; col < n_block_cols(); ++col)
      // This function should throw if the sparsity
      // patterns of the two blocks differ
      block(row, col).add(factor, matrix.block(row, col));
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::operator()(const size_type i,
                                        const size_type j) const
{
  const std::pair<unsigned int, size_type>
    row_index = row_block_indices.global_to_local(i),
    col_index = column_block_indices.global_to_local(j);
  return block(row_index.first, col_index.first)(row_index.second,
                                                 col_index.second);
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::el(const size_type i, const size_type j) const
{
  const std::pair<unsigned int, size_type>
    row_index = row_block_indices.global_to_local(i),
    col_index = column_block_indices.global_to_local(j);
  return block(row_index.first, col_index.first)
    .el(row_index.second, col_index.second);
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::diag_element(const size_type i) const
{
  Assert(n_block_rows() == n_block_cols(), ExcNotQuadratic());

  const std::pair<unsigned int, size_type> index =
    row_block_indices.global_to_local(i);
  return block(index.first, index.first).diag_element(index.second);
}



template <class MatrixType>
inline void
BlockMatrixBase<MatrixType>::compress(
  ::VectorOperation::values operation)
{
  for (unsigned int r = 0; r < n_block_rows(); ++r)
    for (unsigned int c = 0; c < n_block_cols(); ++c)
      block(r, c).compress(operation);
}



template <class MatrixType>
inline BlockMatrixBase<MatrixType> &
BlockMatrixBase<MatrixType>::operator*=(const value_type factor)
{
  Assert(n_block_cols() != 0, ExcNotInitialized());
  Assert(n_block_rows() != 0, ExcNotInitialized());

  for (unsigned int r = 0; r < n_block_rows(); ++r)
    for (unsigned int c = 0; c < n_block_cols(); ++c)
      block(r, c) *= factor;

  return *this;
}



template <class MatrixType>
inline BlockMatrixBase<MatrixType> &
BlockMatrixBase<MatrixType>::operator/=(const value_type factor)
{
  Assert(n_block_cols() != 0, ExcNotInitialized());
  Assert(n_block_rows() != 0, ExcNotInitialized());
  Assert(factor != 0, ExcDivideByZero());

  const value_type factor_inv = 1. / factor;

  for (unsigned int r = 0; r < n_block_rows(); ++r)
    for (unsigned int c = 0; c < n_block_cols(); ++c)
      block(r, c) *= factor_inv;

  return *this;
}



template <class MatrixType>
const BlockIndices &
BlockMatrixBase<MatrixType>::get_row_indices() const
{
  return this->row_block_indices;
}



template <class MatrixType>
const BlockIndices &
BlockMatrixBase<MatrixType>::get_column_indices() const
{
  return this->column_block_indices;
}



template <class MatrixType>
template <class BlockVectorType>
void
BlockMatrixBase<MatrixType>::vmult_block_block(BlockVectorType &      dst,
                                               const BlockVectorType &src) const
{
  Assert(dst.n_blocks() == n_block_rows(),
         ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
  Assert(src.n_blocks() == n_block_cols(),
         ExcDimensionMismatch(src.n_blocks(), n_block_cols()));

  for (size_type row = 0; row < n_block_rows(); ++row)
    {
      block(row, 0).vmult(dst.block(row), src.block(0));
      for (size_type col = 1; col < n_block_cols(); ++col)
        block(row, col).vmult_add(dst.block(row), src.block(col));
    };
}



template <class MatrixType>
template <class BlockVectorType, class VectorType>
void
BlockMatrixBase<MatrixType>::vmult_nonblock_block(
  VectorType &           dst,
  const BlockVectorType &src) const
{
  Assert(n_block_rows() == 1, ExcDimensionMismatch(1, n_block_rows()));
  Assert(src.n_blocks() == n_block_cols(),
         ExcDimensionMismatch(src.n_blocks(), n_block_cols()));

  block(0, 0).vmult(dst, src.block(0));
  for (size_type col = 1; col < n_block_cols(); ++col)
    block(0, col).vmult_add(dst, src.block(col));
}



template <class MatrixType>
template <class BlockVectorType, class VectorType>
void
BlockMatrixBase<MatrixType>::vmult_block_nonblock(BlockVectorType & dst,
                                                  const VectorType &src) const
{
  Assert(dst.n_blocks() == n_block_rows(),
         ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
  Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));

  for (size_type row = 0; row < n_block_rows(); ++row)
    block(row, 0).vmult(dst.block(row), src);
}



template <class MatrixType>
template <class VectorType>
void
BlockMatrixBase<MatrixType>::vmult_nonblock_nonblock(
  VectorType &      dst,
  const VectorType &src) const
{
  Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));
  Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));

  block(0, 0).vmult(dst, src);
}



template <class MatrixType>
template <class BlockVectorType>
void
BlockMatrixBase<MatrixType>::vmult_add(BlockVectorType &      dst,
                                       const BlockVectorType &src) const
{
  Assert(dst.n_blocks() == n_block_rows(),
         ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
  Assert(src.n_blocks() == n_block_cols(),
         ExcDimensionMismatch(src.n_blocks(), n_block_cols()));

  for (unsigned int row = 0; row < n_block_rows(); ++row)
    for (unsigned int col = 0; col < n_block_cols(); ++col)
      block(row, col).vmult_add(dst.block(row), src.block(col));
}



template <class MatrixType>
template <class BlockVectorType>
void
BlockMatrixBase<MatrixType>::Tvmult_block_block(
  BlockVectorType &      dst,
  const BlockVectorType &src) const
{
  Assert(dst.n_blocks() == n_block_cols(),
         ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
  Assert(src.n_blocks() == n_block_rows(),
         ExcDimensionMismatch(src.n_blocks(), n_block_rows()));

  dst = 0.;

  for (unsigned int row = 0; row < n_block_rows(); ++row)
    {
      for (unsigned int col = 0; col < n_block_cols(); ++col)
        block(row, col).Tvmult_add(dst.block(col), src.block(row));
    };
}



template <class MatrixType>
template <class BlockVectorType, class VectorType>
void
BlockMatrixBase<MatrixType>::Tvmult_block_nonblock(BlockVectorType & dst,
                                                   const VectorType &src) const
{
  Assert(dst.n_blocks() == n_block_cols(),
         ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
  Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));

  dst = 0.;

  for (unsigned int col = 0; col < n_block_cols(); ++col)
    block(0, col).Tvmult_add(dst.block(col), src);
}



template <class MatrixType>
template <class BlockVectorType, class VectorType>
void
BlockMatrixBase<MatrixType>::Tvmult_nonblock_block(
  VectorType &           dst,
  const BlockVectorType &src) const
{
  Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
  Assert(src.n_blocks() == n_block_rows(),
         ExcDimensionMismatch(src.n_blocks(), n_block_rows()));

  block(0, 0).Tvmult(dst, src.block(0));

  for (size_type row = 1; row < n_block_rows(); ++row)
    block(row, 0).Tvmult_add(dst, src.block(row));
}



template <class MatrixType>
template <class VectorType>
void
BlockMatrixBase<MatrixType>::Tvmult_nonblock_nonblock(
  VectorType &      dst,
  const VectorType &src) const
{
  Assert(1 == n_block_cols(), ExcDimensionMismatch(1, n_block_cols()));
  Assert(1 == n_block_rows(), ExcDimensionMismatch(1, n_block_rows()));

  block(0, 0).Tvmult(dst, src);
}



template <class MatrixType>
template <class BlockVectorType>
void
BlockMatrixBase<MatrixType>::Tvmult_add(BlockVectorType &      dst,
                                        const BlockVectorType &src) const
{
  Assert(dst.n_blocks() == n_block_cols(),
         ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
  Assert(src.n_blocks() == n_block_rows(),
         ExcDimensionMismatch(src.n_blocks(), n_block_rows()));

  for (unsigned int row = 0; row < n_block_rows(); ++row)
    for (unsigned int col = 0; col < n_block_cols(); ++col)
      block(row, col).Tvmult_add(dst.block(col), src.block(row));
}



template <class MatrixType>
template <class BlockVectorType>
typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::matrix_norm_square(const BlockVectorType &v) const
{
  Assert(n_block_rows() == n_block_cols(), ExcNotQuadratic());
  Assert(v.n_blocks() == n_block_rows(),
         ExcDimensionMismatch(v.n_blocks(), n_block_rows()));

  value_type norm_sqr = 0;
  for (unsigned int row = 0; row < n_block_rows(); ++row)
    for (unsigned int col = 0; col < n_block_cols(); ++col)
      if (row == col)
        norm_sqr += block(row, col).matrix_norm_square(v.block(row));
      else
        norm_sqr +=
          block(row, col).matrix_scalar_product(v.block(row), v.block(col));
  return norm_sqr;
}



template <class MatrixType>
typename BlockMatrixBase<MatrixType>::real_type
BlockMatrixBase<MatrixType>::frobenius_norm() const
{
  value_type norm_sqr = 0;

  // For each block, get the Frobenius norm, and add the square to the
  // accumulator for the full matrix
  for (unsigned int row = 0; row < n_block_rows(); ++row)
    {
      for (unsigned int col = 0; col < n_block_cols(); ++col)
        {
          const value_type block_norm = block(row, col).frobenius_norm();
          norm_sqr += block_norm * block_norm;
        }
    }

  return std::sqrt(norm_sqr);
}



template <class MatrixType>
template <class BlockVectorType>
typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::matrix_scalar_product(
  const BlockVectorType &u,
  const BlockVectorType &v) const
{
  Assert(u.n_blocks() == n_block_rows(),
         ExcDimensionMismatch(u.n_blocks(), n_block_rows()));
  Assert(v.n_blocks() == n_block_cols(),
         ExcDimensionMismatch(v.n_blocks(), n_block_cols()));

  value_type result = 0;
  for (unsigned int row = 0; row < n_block_rows(); ++row)
    for (unsigned int col = 0; col < n_block_cols(); ++col)
      result +=
        block(row, col).matrix_scalar_product(u.block(row), v.block(col));
  return result;
}



template <class MatrixType>
template <class BlockVectorType>
typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::residual(BlockVectorType &      dst,
                                      const BlockVectorType &x,
                                      const BlockVectorType &b) const
{
  Assert(dst.n_blocks() == n_block_rows(),
         ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
  Assert(b.n_blocks() == n_block_rows(),
         ExcDimensionMismatch(b.n_blocks(), n_block_rows()));
  Assert(x.n_blocks() == n_block_cols(),
         ExcDimensionMismatch(x.n_blocks(), n_block_cols()));
  // in block notation, the residual is
  // r_i = b_i - \sum_j A_ij x_j.
  // this can be written as
  // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j.
  //
  // for the first two terms, we can
  // call the residual function of
  // A_i0. for the other terms, we
  // use vmult_add. however, we want
  // to subtract, so in order to
  // avoid a temporary vector, we
  // perform a sign change of the
  // first two term before, and after
  // adding up
  for (unsigned int row = 0; row < n_block_rows(); ++row)
    {
      block(row, 0).residual(dst.block(row), x.block(0), b.block(row));

      for (size_type i = 0; i < dst.block(row).size(); ++i)
        dst.block(row)(i) = -dst.block(row)(i);

      for (unsigned int col = 1; col < n_block_cols(); ++col)
        block(row, col).vmult_add(dst.block(row), x.block(col));

      for (size_type i = 0; i < dst.block(row).size(); ++i)
        dst.block(row)(i) = -dst.block(row)(i);
    };

  value_type res = 0;
  for (size_type row = 0; row < n_block_rows(); ++row)
    res += dst.block(row).norm_sqr();
  return std::sqrt(res);
}



template <class MatrixType>
inline void
BlockMatrixBase<MatrixType>::print(std::ostream &out,
                                   const bool    alternative_output) const
{
  for (unsigned int row = 0; row < n_block_rows(); ++row)
    for (unsigned int col = 0; col < n_block_cols(); ++col)
      {
        if (!alternative_output)
          out << "Block (" << row << ", " << col << ")" << std::endl;

        block(row, col).print(out, alternative_output);
      }
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::const_iterator
BlockMatrixBase<MatrixType>::begin() const
{
  return const_iterator(this, 0);
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::const_iterator
BlockMatrixBase<MatrixType>::end() const
{
  return const_iterator(this, m());
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::const_iterator
BlockMatrixBase<MatrixType>::begin(const size_type r) const
{
  AssertIndexRange(r, m());
  return const_iterator(this, r);
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::const_iterator
BlockMatrixBase<MatrixType>::end(const size_type r) const
{
  AssertIndexRange(r, m());
  return const_iterator(this, r + 1);
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::iterator
BlockMatrixBase<MatrixType>::begin()
{
  return iterator(this, 0);
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::iterator
BlockMatrixBase<MatrixType>::end()
{
  return iterator(this, m());
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::iterator
BlockMatrixBase<MatrixType>::begin(const size_type r)
{
  AssertIndexRange(r, m());
  return iterator(this, r);
}



template <class MatrixType>
inline typename BlockMatrixBase<MatrixType>::iterator
BlockMatrixBase<MatrixType>::end(const size_type r)
{
  AssertIndexRange(r, m());
  return iterator(this, r + 1);
}



template <class MatrixType>
void
BlockMatrixBase<MatrixType>::collect_sizes()
{
  std::vector<size_type> row_sizes(this->n_block_rows());
  std::vector<size_type> col_sizes(this->n_block_cols());

  // first find out the row sizes
  // from the first block column
  for (unsigned int r = 0; r < this->n_block_rows(); ++r)
    row_sizes[r] = sub_objects[r][0]->m();
  // then check that the following
  // block columns have the same
  // sizes
  for (unsigned int c = 1; c < this->n_block_cols(); ++c)
    for (unsigned int r = 0; r < this->n_block_rows(); ++r)
      Assert(row_sizes[r] == sub_objects[r][c]->m(),
             ExcIncompatibleRowNumbers(r, 0, r, c));

  // finally initialize the row
  // indices with this array
  this->row_block_indices.reinit(row_sizes);


  // then do the same with the columns
  for (unsigned int c = 0; c < this->n_block_cols(); ++c)
    col_sizes[c] = sub_objects[0][c]->n();
  for (unsigned int r = 1; r < this->n_block_rows(); ++r)
    for (unsigned int c = 0; c < this->n_block_cols(); ++c)
      Assert(col_sizes[c] == sub_objects[r][c]->n(),
             ExcIncompatibleRowNumbers(0, c, r, c));

  // finally initialize the row
  // indices with this array
  this->column_block_indices.reinit(col_sizes);
}



template <class MatrixType>
void
BlockMatrixBase<MatrixType>::prepare_add_operation()
{
  for (unsigned int row = 0; row < n_block_rows(); ++row)
    for (unsigned int col = 0; col < n_block_cols(); ++col)
      block(row, col).prepare_add();
}



template <class MatrixType>
void
BlockMatrixBase<MatrixType>::prepare_set_operation()
{
  for (unsigned int row = 0; row < n_block_rows(); ++row)
    for (unsigned int col = 0; col < n_block_cols(); ++col)
      block(row, col).prepare_set();
}

#endif // DOXYGEN


DEAL_II_NAMESPACE_CLOSE

#endif // dealii_block_matrix_base_h