block_linear_operator.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2010 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_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;
  }
}

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 <size_t m, 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 <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 <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:

  typedef LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType> 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 <size_t m, size_t n>
  BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
  {
    *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
  }

  template <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 <size_t m, 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 <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
  {
    // 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 = [=](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));

        // Bug: We need a mechanism similar to
        // "apply_with_intermediate_storage" for LinearOperator to do the
        // matrix vector multiplication correctly. Currently, if u and v
        // are equal, the first vmult will garble up the ith block and
        // subsequent multiplications are wrong.
        Assert(!PointerComparison::equal(&v, &u),
               ExcMessage("BlockLinearOperator::vmult currently requires that "
                          "source and destination vectors are different memory "
                          "locations"));

        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 = [=](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));

        // Bug: We need a mechanism similar to
        // "apply_with_intermediate_storage" for LinearOperator to do the
        // matrix vector multiplication correctly. Currently, if u and v
        // are equal, the first vmult will garble up the ith block and
        // subsequent multiplications are wrong.
        Assert(!PointerComparison::equal(&v, &u),
               ExcMessage("BlockLinearOperator::vmult_add currently requires that "
                          "source and destination vectors are different memory "
                          "locations"));

        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 = [=](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));

        // Bug: We need a mechanism similar to
        // "apply_with_intermediate_storage" for LinearOperator to do the
        // matrix vector multiplication correctly. Currently, if u and v
        // are equal, the first vmult will garble up the ith block and
        // subsequent multiplications are wrong.
        Assert(!PointerComparison::equal(&v, &u),
               ExcMessage("BlockLinearOperator::Tvmult currently requires that "
                          "source and destination vectors are different memory "
                          "locations"));

        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 = [=](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));

        // Bug: We need a mechanism similar to
        // "apply_with_intermediate_storage" for LinearOperator to do the
        // matrix vector multiplication correctly. Currently, if u and v
        // are equal, the first vmult will garble up the ith block and
        // subsequent multiplications are wrong.
        Assert(!PointerComparison::equal(&v, &u),
               ExcMessage("BlockLinearOperator::Tvmult_add currently requires that "
                          "source and destination vectors are different memory "
                          "locations"));

        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:
      typedef PayloadBlockType BlockType;

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

  } /*namespace BlockLinearOperator*/
} /*namespace internal*/




template <typename Range,
          typename Domain,
          typename BlockPayload,
          typename BlockMatrixType>
BlockLinearOperator<Range, Domain, BlockPayload>
block_operator(const BlockMatrixType &block_matrix)
{
  typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType 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));
  };

  internal::BlockLinearOperatorImplementation::populate_linear_operator_functions(return_op);
  return return_op;
}



template <size_t m, 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");

  typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;

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

  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];
  };

  internal::BlockLinearOperatorImplementation::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)
{
  typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType 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)));
  };

  internal::BlockLinearOperatorImplementation::populate_linear_operator_functions(return_op);
  return return_op;
}



template <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");

  typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType 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 <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");

  typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType 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