block_matrix_base.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2023 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/mutex.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/table.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>
#include <mutex>
 
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;
 
    Threads::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