block_linear_operator.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2010 - 2019 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_linear_operator_h
#define dealii_block_linear_operator_h

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

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

#include <deal.II/lac/linear_operator.h>


DEAL_II_NAMESPACE_OPEN

// Forward declarations:

namespace internal
{
  namespace BlockLinearOperatorImplementation
  {
    template <typename PayloadBlockType =
                internal::LinearOperatorImplementation::EmptyPayload>
    class EmptyBlockPayload;
  }
} // namespace internal

template <typename Number>
class BlockVector;

template <typename Range  = BlockVector<double>,
          typename Domain = Range,
          typename BlockPayload =
            internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
class BlockLinearOperator;

template <typename Range  = BlockVector<double>,
          typename Domain = Range,
          typename BlockPayload =
            internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
          typename BlockMatrixType>
BlockLinearOperator<Range, Domain, BlockPayload>
block_operator(const BlockMatrixType &matrix);

template <std::size_t m,
          std::size_t n,
          typename Range  = BlockVector<double>,
          typename Domain = Range,
          typename BlockPayload =
            internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
BlockLinearOperator<Range, Domain, BlockPayload>
block_operator(
  const std::array<std::array<LinearOperator<typename Range::BlockType,
                                             typename Domain::BlockType,
                                             typename BlockPayload::BlockType>,
                              n>,
                   m> &);

template <std::size_t m,
          typename Range  = BlockVector<double>,
          typename Domain = Range,
          typename BlockPayload =
            internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
BlockLinearOperator<Range, Domain, BlockPayload>
block_diagonal_operator(
  const std::array<LinearOperator<typename Range::BlockType,
                                  typename Domain::BlockType,
                                  typename BlockPayload::BlockType>,
                   m> &);

template <std::size_t m,
          typename Range  = BlockVector<double>,
          typename Domain = Range,
          typename BlockPayload =
            internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
BlockLinearOperator<Range, Domain, BlockPayload>
block_diagonal_operator(
  const LinearOperator<typename Range::BlockType,
                       typename Domain::BlockType,
                       typename BlockPayload::BlockType> &op);

// This is a workaround for a bug in <=gcc-4.7 that does not like partial
// template default values in combination with local lambda expressions [1]
//
// [1] https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53624
//
// Forward declare functions with partial template defaults:

template <typename Range  = BlockVector<double>,
          typename Domain = Range,
          typename BlockPayload =
            internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
          typename BlockMatrixType>
BlockLinearOperator<Range, Domain, BlockPayload>
block_diagonal_operator(const BlockMatrixType &block_matrix);

template <typename Range  = BlockVector<double>,
          typename Domain = Range,
          typename BlockPayload =
            internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
LinearOperator<Domain, Range, typename BlockPayload::BlockType>
block_forward_substitution(
  const BlockLinearOperator<Range, Domain, BlockPayload> &,
  const BlockLinearOperator<Domain, Range, BlockPayload> &);

template <typename Range  = BlockVector<double>,
          typename Domain = Range,
          typename BlockPayload =
            internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
LinearOperator<Domain, Range, typename BlockPayload::BlockType>
block_back_substitution(
  const BlockLinearOperator<Range, Domain, BlockPayload> &,
  const BlockLinearOperator<Domain, Range, BlockPayload> &);

// end of workaround



template <typename Range, typename Domain, typename BlockPayload>
class BlockLinearOperator
  : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
{
public:
  using BlockType = LinearOperator<typename Range::BlockType,
                                   typename Domain::BlockType,
                                   typename BlockPayload::BlockType>;

  BlockLinearOperator(const BlockPayload &payload)
    : LinearOperator<Range, Domain, typename BlockPayload::BlockType>(
        typename BlockPayload::BlockType(payload, payload))
  {
    n_block_rows = []() -> unsigned int {
      Assert(
        false,
        ExcMessage(
          "Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
      return 0;
    };

    n_block_cols = []() -> unsigned int {
      Assert(
        false,
        ExcMessage(
          "Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
      return 0;
    };

    block = [](unsigned int, unsigned int) -> BlockType {
      Assert(
        false,
        ExcMessage(
          "Uninitialized BlockLinearOperator<Range, Domain>::block called"));
      return BlockType();
    };
  }

  BlockLinearOperator(
    const BlockLinearOperator<Range, Domain, BlockPayload> &) = default;

  template <typename Op>
  BlockLinearOperator(const Op &op)
  {
    *this = block_operator<Range, Domain, BlockPayload, Op>(op);
  }

  template <std::size_t m, std::size_t n>
  BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
  {
    *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
  }

  template <std::size_t m>
  BlockLinearOperator(const std::array<BlockType, m> &ops)
  {
    *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
  }

  BlockLinearOperator<Range, Domain, BlockPayload> &
  operator=(const BlockLinearOperator<Range, Domain, BlockPayload> &) = default;

  template <typename Op>
  BlockLinearOperator<Range, Domain, BlockPayload> &
  operator=(const Op &op)
  {
    *this = block_operator<Range, Domain, BlockPayload, Op>(op);
    return *this;
  }

  template <std::size_t m, std::size_t n>
  BlockLinearOperator<Range, Domain, BlockPayload> &
  operator=(const std::array<std::array<BlockType, n>, m> &ops)
  {
    *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
    return *this;
  }

  template <std::size_t m>
  BlockLinearOperator<Range, Domain, BlockPayload> &
  operator=(const std::array<BlockType, m> &ops)
  {
    *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
    return *this;
  }

  std::function<unsigned int()> n_block_rows;

  std::function<unsigned int()> n_block_cols;

  std::function<BlockType(unsigned int, unsigned int)> block;
};


namespace internal
{
  namespace BlockLinearOperatorImplementation
  {
    // A helper function to apply a given vmult, or Tvmult to a vector with
    // intermediate storage, similar to the corresponding helper
    // function for LinearOperator. Here, two operators are used.
    // The first one takes care of the first "column" and typically doesn't add.
    // On the other hand, the second operator is normally an adding one.
    template <typename Function1,
              typename Function2,
              typename Range,
              typename Domain>
    void
    apply_with_intermediate_storage(const Function1 &first_op,
                                    const Function2 &loop_op,
                                    Range &          v,
                                    const Domain &   u,
                                    bool             add)
    {
      GrowingVectorMemory<Range> vector_memory;

      typename VectorMemory<Range>::Pointer tmp(vector_memory);
      tmp->reinit(v, /*bool omit_zeroing_entries =*/true);

      const unsigned int n = u.n_blocks();
      const unsigned int m = v.n_blocks();

      for (unsigned int i = 0; i < m; ++i)
        {
          first_op(*tmp, u, i, 0);
          for (unsigned int j = 1; j < n; ++j)
            loop_op(*tmp, u, i, j);
        }

      if (add)
        v += *tmp;
      else
        v = *tmp;
    }

    // Populate the LinearOperator interfaces with the help of the
    // BlockLinearOperator functions
    template <typename Range, typename Domain, typename BlockPayload>
    inline void
    populate_linear_operator_functions(
      ::BlockLinearOperator<Range, Domain, BlockPayload> &op)
    {
      op.reinit_range_vector = [=](Range &v, bool omit_zeroing_entries) {
        const unsigned int m = op.n_block_rows();

        // Reinitialize the block vector to m blocks:
        v.reinit(m);

        // And reinitialize every individual block with reinit_range_vectors:
        for (unsigned int i = 0; i < m; ++i)
          op.block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);

        v.collect_sizes();
      };

      op.reinit_domain_vector = [=](Domain &v, bool omit_zeroing_entries) {
        const unsigned int n = op.n_block_cols();

        // Reinitialize the block vector to n blocks:
        v.reinit(n);

        // And reinitialize every individual block with reinit_domain_vectors:
        for (unsigned int i = 0; i < n; ++i)
          op.block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);

        v.collect_sizes();
      };

      op.vmult = [&op](Range &v, const Domain &u) {
        const unsigned int m = op.n_block_rows();
        const unsigned int n = op.n_block_cols();
        Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
        Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));

        if (PointerComparison::equal(&v, &u))
          {
            const auto first_op = [&op](Range &            v,
                                        const Domain &     u,
                                        const unsigned int i,
                                        const unsigned int j) {
              op.block(i, j).vmult(v.block(i), u.block(j));
            };

            const auto loop_op = [&op](Range &            v,
                                       const Domain &     u,
                                       const unsigned int i,
                                       const unsigned int j) {
              op.block(i, j).vmult_add(v.block(i), u.block(j));
            };

            apply_with_intermediate_storage(first_op, loop_op, v, u, false);
          }
        else
          {
            for (unsigned int i = 0; i < m; ++i)
              {
                op.block(i, 0).vmult(v.block(i), u.block(0));
                for (unsigned int j = 1; j < n; ++j)
                  op.block(i, j).vmult_add(v.block(i), u.block(j));
              }
          }
      };

      op.vmult_add = [&op](Range &v, const Domain &u) {
        const unsigned int m = op.n_block_rows();
        const unsigned int n = op.n_block_cols();
        Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
        Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));

        if (PointerComparison::equal(&v, &u))
          {
            const auto first_op = [&op](Range &            v,
                                        const Domain &     u,
                                        const unsigned int i,
                                        const unsigned int j) {
              op.block(i, j).vmult(v.block(i), u.block(j));
            };

            const auto loop_op = [&op](Range &            v,
                                       const Domain &     u,
                                       const unsigned int i,
                                       const unsigned int j) {
              op.block(i, j).vmult_add(v.block(i), u.block(j));
            };

            apply_with_intermediate_storage(first_op, loop_op, v, u, true);
          }
        else
          {
            for (unsigned int i = 0; i < m; ++i)
              for (unsigned int j = 0; j < n; ++j)
                op.block(i, j).vmult_add(v.block(i), u.block(j));
          }
      };

      op.Tvmult = [&op](Domain &v, const Range &u) {
        const unsigned int n = op.n_block_cols();
        const unsigned int m = op.n_block_rows();
        Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
        Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));

        if (PointerComparison::equal(&v, &u))
          {
            const auto first_op = [&op](Range &            v,
                                        const Domain &     u,
                                        const unsigned int i,
                                        const unsigned int j) {
              op.block(j, i).Tvmult(v.block(i), u.block(j));
            };

            const auto loop_op = [&op](Range &            v,
                                       const Domain &     u,
                                       const unsigned int i,
                                       const unsigned int j) {
              op.block(j, i).Tvmult_add(v.block(i), u.block(j));
            };

            apply_with_intermediate_storage(first_op, loop_op, v, u, false);
          }
        else
          {
            for (unsigned int i = 0; i < n; ++i)
              {
                op.block(0, i).Tvmult(v.block(i), u.block(0));
                for (unsigned int j = 1; j < m; ++j)
                  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
              }
          }
      };

      op.Tvmult_add = [&op](Domain &v, const Range &u) {
        const unsigned int n = op.n_block_cols();
        const unsigned int m = op.n_block_rows();
        Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
        Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));

        if (PointerComparison::equal(&v, &u))
          {
            const auto first_op = [&op](Range &            v,
                                        const Domain &     u,
                                        const unsigned int i,
                                        const unsigned int j) {
              op.block(j, i).Tvmult(v.block(i), u.block(j));
            };

            const auto loop_op = [&op](Range &            v,
                                       const Domain &     u,
                                       const unsigned int i,
                                       const unsigned int j) {
              op.block(j, i).Tvmult_add(v.block(i), u.block(j));
            };

            apply_with_intermediate_storage(first_op, loop_op, v, u, true);
          }
        else
          {
            for (unsigned int i = 0; i < n; ++i)
              for (unsigned int j = 0; j < m; ++j)
                op.block(j, i).Tvmult_add(v.block(i), u.block(j));
          }
      };
    }



    template <typename PayloadBlockType>
    class EmptyBlockPayload
    {
    public:
      using BlockType = PayloadBlockType;

      template <typename... Args>
      EmptyBlockPayload(const Args &...)
      {}
    };

  } // namespace BlockLinearOperatorImplementation
} // namespace internal




template <typename Range,
          typename Domain,
          typename BlockPayload,
          typename BlockMatrixType>
BlockLinearOperator<Range, Domain, BlockPayload>
block_operator(const BlockMatrixType &block_matrix)
{
  using BlockType =
    typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType;

  BlockLinearOperator<Range, Domain, BlockPayload> return_op{
    BlockPayload(block_matrix, block_matrix)};

  return_op.n_block_rows = [&block_matrix]() -> unsigned int {
    return block_matrix.n_block_rows();
  };

  return_op.n_block_cols = [&block_matrix]() -> unsigned int {
    return block_matrix.n_block_cols();
  };

  return_op.block = [&block_matrix](unsigned int i,
                                    unsigned int j) -> BlockType {
#ifdef DEBUG
    const unsigned int m = block_matrix.n_block_rows();
    const unsigned int n = block_matrix.n_block_cols();
    Assert(i < m, ExcIndexRange(i, 0, m));
    Assert(j < n, ExcIndexRange(j, 0, n));
#endif

    return BlockType(block_matrix.block(i, j));
  };

  populate_linear_operator_functions(return_op);
  return return_op;
}



template <std::size_t m,
          std::size_t n,
          typename Range,
          typename Domain,
          typename BlockPayload>
BlockLinearOperator<Range, Domain, BlockPayload>
block_operator(
  const std::array<std::array<LinearOperator<typename Range::BlockType,
                                             typename Domain::BlockType,
                                             typename BlockPayload::BlockType>,
                              n>,
                   m> &ops)
{
  static_assert(m > 0 && n > 0,
                "a blocked LinearOperator must consist of at least one block");

  using BlockType =
    typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType;

  // TODO: Create block payload so that this can be initialized correctly
  BlockLinearOperator<Range, Domain, BlockPayload> return_op{BlockPayload()};

  return_op.n_block_rows = []() -> unsigned int { return m; };

  return_op.n_block_cols = []() -> unsigned int { return n; };

  return_op.block = [ops](unsigned int i, unsigned int j) -> BlockType {
    Assert(i < m, ExcIndexRange(i, 0, m));
    Assert(j < n, ExcIndexRange(j, 0, n));

    return ops[i][j];
  };

  populate_linear_operator_functions(return_op);
  return return_op;
}



template <typename Range,
          typename Domain,
          typename BlockPayload,
          typename BlockMatrixType>
BlockLinearOperator<Range, Domain, BlockPayload>
block_diagonal_operator(const BlockMatrixType &block_matrix)
{
  using BlockType =
    typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType;

  BlockLinearOperator<Range, Domain, BlockPayload> return_op{
    BlockPayload(block_matrix, block_matrix)};

  return_op.n_block_rows = [&block_matrix]() -> unsigned int {
    return block_matrix.n_block_rows();
  };

  return_op.n_block_cols = [&block_matrix]() -> unsigned int {
    return block_matrix.n_block_cols();
  };

  return_op.block = [&block_matrix](unsigned int i,
                                    unsigned int j) -> BlockType {
#ifdef DEBUG
    const unsigned int m = block_matrix.n_block_rows();
    const unsigned int n = block_matrix.n_block_cols();
    Assert(m == n, ExcDimensionMismatch(m, n));
    Assert(i < m, ExcIndexRange(i, 0, m));
    Assert(j < n, ExcIndexRange(j, 0, n));
#endif
    if (i == j)
      return BlockType(block_matrix.block(i, j));
    else
      return null_operator(BlockType(block_matrix.block(i, j)));
  };

  populate_linear_operator_functions(return_op);
  return return_op;
}



template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
BlockLinearOperator<Range, Domain, BlockPayload>
block_diagonal_operator(
  const std::array<LinearOperator<typename Range::BlockType,
                                  typename Domain::BlockType,
                                  typename BlockPayload::BlockType>,
                   m> &ops)
{
  static_assert(
    m > 0, "a blockdiagonal LinearOperator must consist of at least one block");

  using BlockType =
    typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType;

  std::array<std::array<BlockType, m>, m> new_ops;

  // This is a bit tricky. We have to make sure that the off-diagonal
  // elements of return_op.ops are populated correctly. They must be
  // null_operators, but with correct reinit_domain_vector and
  // reinit_range_vector functions.
  for (unsigned int i = 0; i < m; ++i)
    for (unsigned int j = 0; j < m; ++j)
      if (i == j)
        {
          // diagonal elements are easy:
          new_ops[i][j] = ops[i];
        }
      else
        {
          // create a null-operator...
          new_ops[i][j] = null_operator(ops[i]);
          // ... and fix up reinit_domain_vector:
          new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
        }

  return block_operator<m, m, Range, Domain>(new_ops);
}



template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
BlockLinearOperator<Range, Domain, BlockPayload>
block_diagonal_operator(
  const LinearOperator<typename Range::BlockType,
                       typename Domain::BlockType,
                       typename BlockPayload::BlockType> &op)
{
  static_assert(m > 0,
                "a blockdiagonal LinearOperator must consist of at least "
                "one block");

  using BlockType =
    typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType;
  std::array<BlockType, m> new_ops;
  new_ops.fill(op);

  return block_diagonal_operator(new_ops);
}





template <typename Range, typename Domain, typename BlockPayload>
LinearOperator<Domain, Range, typename BlockPayload::BlockType>
block_forward_substitution(
  const BlockLinearOperator<Range, Domain, BlockPayload> &block_operator,
  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
{
  LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op{
    typename BlockPayload::BlockType(diagonal_inverse)};

  return_op.reinit_range_vector  = diagonal_inverse.reinit_range_vector;
  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;

  return_op.vmult = [block_operator, diagonal_inverse](Range &      v,
                                                       const Range &u) {
    const unsigned int m = block_operator.n_block_rows();
    Assert(block_operator.n_block_cols() == m,
           ExcDimensionMismatch(block_operator.n_block_cols(), m));
    Assert(diagonal_inverse.n_block_rows() == m,
           ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
    Assert(diagonal_inverse.n_block_cols() == m,
           ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));

    if (m == 0)
      return;

    diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
    for (unsigned int i = 1; i < m; ++i)
      {
        auto &dst = v.block(i);
        dst       = u.block(i);
        dst *= -1.;
        for (unsigned int j = 0; j < i; ++j)
          block_operator.block(i, j).vmult_add(dst, v.block(j));
        dst *= -1.;
        diagonal_inverse.block(i, i).vmult(dst,
                                           dst); // uses intermediate storage
      }
  };

  return_op.vmult_add = [block_operator, diagonal_inverse](Range &      v,
                                                           const Range &u) {
    const unsigned int m = block_operator.n_block_rows();
    Assert(block_operator.n_block_cols() == m,
           ExcDimensionMismatch(block_operator.n_block_cols(), m));
    Assert(diagonal_inverse.n_block_rows() == m,
           ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
    Assert(diagonal_inverse.n_block_cols() == m,
           ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));

    if (m == 0)
      return;

    GrowingVectorMemory<typename Range::BlockType>            vector_memory;
    typename VectorMemory<typename Range::BlockType>::Pointer tmp(
      vector_memory);

    diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0));

    for (unsigned int i = 1; i < m; ++i)
      {
        diagonal_inverse.block(i, i).reinit_range_vector(
          *tmp, /*bool omit_zeroing_entries=*/true);
        *tmp = u.block(i);
        *tmp *= -1.;
        for (unsigned int j = 0; j < i; ++j)
          block_operator.block(i, j).vmult_add(*tmp, v.block(j));
        *tmp *= -1.;
        diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
      }
  };

  return return_op;
}



template <typename Range, typename Domain, typename BlockPayload>
LinearOperator<Domain, Range, typename BlockPayload::BlockType>
block_back_substitution(
  const BlockLinearOperator<Range, Domain, BlockPayload> &block_operator,
  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
{
  LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op{
    typename BlockPayload::BlockType(diagonal_inverse)};

  return_op.reinit_range_vector  = diagonal_inverse.reinit_range_vector;
  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;

  return_op.vmult = [block_operator, diagonal_inverse](Range &      v,
                                                       const Range &u) {
    const unsigned int m = block_operator.n_block_rows();
    Assert(block_operator.n_block_cols() == m,
           ExcDimensionMismatch(block_operator.n_block_cols(), m));
    Assert(diagonal_inverse.n_block_rows() == m,
           ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
    Assert(diagonal_inverse.n_block_cols() == m,
           ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));

    if (m == 0)
      return;

    diagonal_inverse.block(m - 1, m - 1).vmult(v.block(m - 1), u.block(m - 1));

    for (int i = m - 2; i >= 0; --i)
      {
        auto &dst = v.block(i);
        dst       = u.block(i);
        dst *= -1.;
        for (unsigned int j = i + 1; j < m; ++j)
          block_operator.block(i, j).vmult_add(dst, v.block(j));
        dst *= -1.;
        diagonal_inverse.block(i, i).vmult(dst,
                                           dst); // uses intermediate storage
      }
  };

  return_op.vmult_add = [block_operator, diagonal_inverse](Range &      v,
                                                           const Range &u) {
    const unsigned int m = block_operator.n_block_rows();
    Assert(block_operator.n_block_cols() == m,
           ExcDimensionMismatch(block_operator.n_block_cols(), m));
    Assert(diagonal_inverse.n_block_rows() == m,
           ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
    Assert(diagonal_inverse.n_block_cols() == m,
           ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
    GrowingVectorMemory<typename Range::BlockType>            vector_memory;
    typename VectorMemory<typename Range::BlockType>::Pointer tmp(
      vector_memory);

    if (m == 0)
      return;

    diagonal_inverse.block(m - 1, m - 1)
      .vmult_add(v.block(m - 1), u.block(m - 1));

    for (int i = m - 2; i >= 0; --i)
      {
        diagonal_inverse.block(i, i).reinit_range_vector(
          *tmp, /*bool omit_zeroing_entries=*/true);
        *tmp = u.block(i);
        *tmp *= -1.;
        for (unsigned int j = i + 1; j < m; ++j)
          block_operator.block(i, j).vmult_add(*tmp, v.block(j));
        *tmp *= -1.;
        diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
      }
  };

  return return_op;
}


DEAL_II_NAMESPACE_CLOSE

#endif