From: Matthias Maier <tamiko@43-1.org>
Date: Sat, 8 Aug 2015 02:19:08 +0000 (-0500)
Subject: A BlockLinearOperator class for storing linear operator block structures
X-Git-Tag: v8.4.0-rc2~536^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b8fa7b12ef6ec018a11320b13dc20fbf776caaf2;p=dealii.git

A BlockLinearOperator class for storing linear operator block structures
---

diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h
index 61c84c22bd..690b132a0e 100644
--- a/include/deal.II/lac/block_linear_operator.h
+++ b/include/deal.II/lac/block_linear_operator.h
@@ -26,23 +26,371 @@
 DEAL_II_NAMESPACE_OPEN
 
 // Forward declarations:
+
 template <typename Number> class BlockVector;
 
+template <typename Range = BlockVector<double>,
+          typename Domain = Range>
+class BlockLinearOperator;
+
+template <typename Range = BlockVector<double>,
+          typename Domain = Range,
+          typename BlockMatrix>
+BlockLinearOperator<Range, Domain>
+block_operator(const BlockMatrix &matrix);
+
+template <size_t m, size_t n,
+          typename Range = BlockVector<double>,
+          typename Domain = Range>
+BlockLinearOperator<Range, Domain>
+block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &);
+
+template <size_t m,
+          typename Range = BlockVector<double>,
+          typename Domain = Range>
+BlockLinearOperator<Range, Domain>
+block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &);
+
+template <size_t m,
+          typename Range = BlockVector<double>,
+          typename Domain = Range>
+BlockLinearOperator<Range, Domain>
+block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::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 BlockMatrix>
+BlockLinearOperator<Range, Domain>
+block_diagonal_operator(const BlockMatrix &block_matrix);
+
+template <typename Range = BlockVector<double>,
+          typename Domain = Range>
+LinearOperator<Domain, Range>
+block_forward_substitution(const BlockLinearOperator<Range, Domain> &,
+                           const BlockLinearOperator<Domain, Range> &);
+
+template <typename Range = BlockVector<double>,
+          typename Domain = Range>
+LinearOperator<Domain, Range>
+block_back_substitution(const BlockLinearOperator<Range, Domain> &,
+                        const BlockLinearOperator<Domain, Range> &);
+
+// end of workaround
+
+
+
+/**
+ *  TODO: Description
+ */
+template <typename Range, typename Domain>
+class BlockLinearOperator : public LinearOperator<Range, Domain>
+{
+public:
+
+  typedef LinearOperator<typename Range::BlockType, typename Domain::BlockType> BlockType;
+
+  /**
+   * Create an empty BlockLinearOperator object.
+   *
+   * All<code>std::function</code> member objects of this class and its
+   * base class LinearOperator are initialized with default variants that
+   * throw an exception upon invocation.
+   */
+  BlockLinearOperator()
+    : LinearOperator<Range, Domain>()
+  {
+
+    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();
+    };
+  }
+
+  /**
+   * Default copy constructor.
+   */
+  BlockLinearOperator(const BlockLinearOperator<Range, Domain> &) =
+    default;
+
+  /**
+   * Templated copy constructor that creates a BlockLinearOperator object
+   * from an object @p op for which the conversion function
+   * <code>block_operator</code> is defined.
+   */
+  template<typename Op>
+  BlockLinearOperator(const Op &op)
+  {
+    *this = block_operator<Range, Domain, Op>(op);
+  }
+
+  /**
+   * Create a BlockLinearOperator from a two-dimensional array @p ops of
+   * LinearOperator. This constructor calls the corresponding
+   * block_operator() specialization.
+   */
+  template<size_t m, size_t n>
+  BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
+  {
+    *this = block_operator<m, n, Range, Domain>(ops);
+  }
+
+  /**
+   * Create a block-diagonal BlockLinearOperator from a one-dimensional
+   * array @p ops of LinearOperator. This constructor calls the
+   * corresponding block_operator() specialization.
+   */
+  template<size_t m>
+  BlockLinearOperator(const std::array<BlockType, m> &ops)
+  {
+    *this = block_diagonal_operator<m, Range, Domain>(ops);
+  }
+
+  /**
+   * Default copy assignment operator.
+   */
+  BlockLinearOperator<Range, Domain> &
+  operator=(const BlockLinearOperator<Range, Domain> &) = default;
+
+  /**
+   * Templated copy assignment operator for an object @p op for which the
+   * conversion function <code>block_operator</code> is defined.
+   */
+  template <typename Op>
+  BlockLinearOperator<Range, Domain> &operator=(const Op &op)
+  {
+    *this = block_operator<Range, Domain, Op>(op);
+    return *this;
+  }
+
+  /**
+   * Copy assignment from a two-dimensional array @p ops of LinearOperator.
+   * This assignment operator calls the corresponding block_operator()
+   * specialization.
+   */
+  template <size_t m, size_t n>
+  BlockLinearOperator<Range, Domain> &
+  operator=(const std::array<std::array<BlockType, n>, m> &ops)
+  {
+    *this = block_operator<m, n, Range, Domain>(ops);
+    return *this;
+  }
+
+  /**
+   * Copy assignment from a one-dimensional array @p ops of LinearOperator
+   * that creates a block-diagonal BlockLinearOperator.
+   * This assignment operator calls the corresponding block_operator()
+   * specialization.
+   */
+  template <size_t m>
+  BlockLinearOperator<Range, Domain> &
+  operator=(const std::array<BlockType, m> &ops)
+  {
+    *this = block_diagonal_operator<m, Range, Domain>(ops);
+    return *this;
+  }
+
+  /**
+   * TODO: Description
+   */
+  std::function<unsigned int()> n_block_rows;
+
+  /**
+   * TODO: Description
+   */
+  std::function<unsigned int()> n_block_cols;
+
+  /**
+   * TODO: Description
+   */
+  std::function<BlockType(unsigned int, unsigned int)> block;
+
+  //@}
+};
+
+
+
+namespace internal
+{
+  namespace BlockLinearOperator
+  {
+    // Populate the LinearOperator interfaces with the help of the
+    // BlockLinearOperator functions
+    template <typename Range, typename Domain>
+    inline void
+    populate_linear_operator_functions(
+      dealii::BlockLinearOperator<Range, Domain> &op)
+    {
+      op.reinit_range_vector = [=](Range &v, bool fast)
+      {
+        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 (size_t i = 0; i < m; ++i)
+          op.block(i, 0).reinit_range_vector(v.block(i), fast);
+
+        v.collect_sizes();
+      };
+
+      op.reinit_domain_vector = [=](Domain &v, bool fast)
+      {
+        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 (size_t i = 0; i < n; ++i)
+          op.block(0, i).reinit_domain_vector(v.block(i), fast);
+
+        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));
+
+        for (size_t i = 0; i < m; ++i)
+          {
+            op.block(i, 0).vmult(v.block(i), u.block(0));
+            for (size_t 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));
+
+        for (size_t i = 0; i < m; ++i)
+          for (size_t 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));
+
+        for (size_t i = 0; i < n; ++i)
+          {
+            op.block(0, i).Tvmult(v.block(i), u.block(0));
+            for (size_t 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));
+
+        for (size_t i = 0; i < n; ++i)
+          for (size_t j = 0; j < m; ++j)
+            op.block(j, i).Tvmult_add(v.block(i), u.block(j));
+      };
+    }
+  } /*namespace BlockLinearOperator*/
+} /*namespace internal*/
+
+
+
 /**
- * @name Creation of LinearOperator block structures
+ * @name Creation of a BlockLinearOperator
  */
 //@{
 
 /**
  * @relates LinearOperator
  *
- * A function that encapsulates a given collection @p ops of LinearOperators
- * into a block structure. Hereby, it is assumed that Range and Domain are
- * blockvectors, i.e., derived from
- * @ref BlockVectorBase.
- * The individual linear operators in @p ops must act on a the underlying
- * vector type of the block vectors, i.e., on Domain::BlockType yielding a
- * result in Range::BlockType.
+ * A function that encapsulates a @p block_matrix into a
+ * BlockLinearOperator.
+ *
+ * All changes made on the block structure and individual blocks of @p
+ * block_matrix after the creation of the BlockLinearOperator object are
+ * reflected by the operator object.
+ *
+ * @ingroup LAOperators
+ */
+template <typename Range,
+          typename Domain,
+          typename BlockMatrix>
+BlockLinearOperator<Range, Domain>
+block_operator(const BlockMatrix &block_matrix)
+{
+  typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+
+  BlockLinearOperator<Range, Domain> return_op;
+
+  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::BlockLinearOperator::populate_linear_operator_functions(return_op);
+  return return_op;
+}
+
+
+
+/**
+ * @relates LinearOperator
+ *
+ * A variant of above function that encapsulates a given collection @p ops
+ * of LinearOperators into a block structure. Here, it is assumed that
+ * Range and Domain are blockvectors, i.e., derived from @ref
+ * BlockVectorBase. The individual linear operators in @p ops must act on
+ * the underlying vector type of the block vectors, i.e., on
+ * Domain::BlockType yielding a result in Range::BlockType.
  *
  * The list @p ops is best passed as an initializer list. Consider for example
  * a linear operator block (acting on Vector<double>)
@@ -60,91 +408,97 @@ template <typename Number> class BlockVector;
  *
  * @ingroup LAOperators
  */
-template <unsigned int m, unsigned int n,
-          typename Range = BlockVector<double>,
-          typename Domain = Range>
-LinearOperator<Range, Domain>
+template <size_t m, size_t n, typename Range, typename Domain>
+BlockLinearOperator<Range, Domain>
 block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &ops)
 {
   static_assert(m > 0 && n > 0,
                 "a blocked LinearOperator must consist of at least one block");
 
-  LinearOperator<Range, Domain> return_op;
-
-  return_op.reinit_range_vector = [ops](Range &v, bool fast)
-  {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(m);
+  typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
 
-    // And reinitialize every individual block with reinit_range_vectors:
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i][0].reinit_range_vector(v.block(i), fast);
+  BlockLinearOperator<Range, Domain> return_op;
 
-    v.collect_sizes();
+  return_op.n_block_rows = []() -> unsigned int
+  {
+    return m;
   };
 
-  return_op.reinit_domain_vector = [ops](Domain &v, bool fast)
+  return_op.n_block_cols = []() -> unsigned int
   {
-    // 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)
-      ops[0][i].reinit_domain_vector(v.block(i), fast);
-
-    v.collect_sizes();
+    return n;
   };
 
-  return_op.vmult = [ops](Range &v, const Domain &u)
+  return_op.block = [ops](unsigned int i, unsigned int j) -> BlockType
   {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
+    Assert(i < m, ExcIndexRange (i, 0, m));
+    Assert(j < n, ExcIndexRange (j, 0, n));
 
-    for (unsigned int i = 0; i < m; ++i)
-      {
-        ops[i][0].vmult(v.block(i), u.block(0));
-        for (unsigned int j = 1; j < n; ++j)
-          ops[i][j].vmult_add(v.block(i), u.block(j));
-      }
+    return ops[i][j];
   };
 
-  return_op.vmult_add = [ops](Range &v, const Domain &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
+  internal::BlockLinearOperator::populate_linear_operator_functions(return_op);
+  return return_op;
+}
 
-    for (unsigned int i = 0; i < m; ++i)
-      for (unsigned int j = 0; j < n; ++j)
-        ops[i][j].vmult_add(v.block(i), u.block(j));
-  };
 
-  return_op.Tvmult = [ops](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
 
-    for (unsigned int i = 0; i < n; ++i)
-      {
-        ops[0][i].Tvmult(v.block(i), u.block(0));
-        for (unsigned int j = 1; j < m; ++j)
-          ops[j][i].Tvmult_add(v.block(i), u.block(j));
-      }
+/**
+ * @relates LinearOperator
+ *
+ * This function extracts the diagonal blocks of @p block_matrix (either a
+ * block matrix type or a BlockLinearOperator) and creates a
+ * BlockLinearOperator with the diagonal. Off-diagonal elements are
+ * initialized as null_operator (with correct reinit_range_vector and
+ * reinit_domain_vector methods).
+ *
+ * All changes made on the individual diagonal blocks of @p block_matrix
+ * after the creation of the BlockLinearOperator object are reflected by
+ * the operator object.
+ *
+ * @ingroup LAOperators
+ */
+template <typename Range,
+          typename Domain,
+          typename BlockMatrix>
+BlockLinearOperator<Range, Domain>
+block_diagonal_operator(const BlockMatrix &block_matrix)
+{
+  typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+
+  BlockLinearOperator<Range, Domain> return_op;
+
+  return_op.n_block_rows = [&block_matrix]() -> unsigned int
+  {
+    return block_matrix.n_block_rows();
   };
 
-  return_op.Tvmult_add = [ops](Domain &v, const Range &u)
+  return_op.n_block_cols = [&block_matrix]() -> unsigned int
   {
-    Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+    return block_matrix.n_block_cols();
+  };
 
-    for (unsigned int i = 0; i < n; ++i)
-      for (unsigned int j = 0; j < m; ++j)
-        ops[j][i].Tvmult_add(v.block(i), u.block(j));
+  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::BlockLinearOperator::populate_linear_operator_functions(return_op);
   return return_op;
 }
 
 
+
 /**
  * @relates LinearOperator
  *
@@ -162,81 +516,41 @@ block_operator(const std::array<std::array<LinearOperator<typename Range::BlockT
  *
  * @ingroup LAOperators
  */
-template <unsigned int m,
-          typename Range = BlockVector<double>,
-          typename Domain = Range>
-LinearOperator<Range, Domain>
+template <size_t m, typename Range, typename Domain>
+BlockLinearOperator<Range, Domain>
 block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &ops)
 {
   static_assert(m > 0,
                 "a blockdiagonal LinearOperator must consist of at least one block");
 
-  LinearOperator<Range, Domain> return_op;
-
-  return_op.reinit_range_vector = [ops](Range &v, bool fast)
-  {
-    // 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)
-      ops[i].reinit_range_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.reinit_domain_vector = [ops](Domain &v, bool fast)
-  {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(m);
-
-    // And reinitialize every individual block with reinit_domain_vectors:
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].reinit_domain_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.vmult = [ops](Range &v, const Domain &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].vmult(v.block(i), u.block(i));
-  };
-
-  return_op.vmult_add = [ops](Range &v, const Domain &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].vmult_add(v.block(i), u.block(i));
-  };
-
-  return_op.Tvmult = [ops](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].Tvmult(v.block(i), u.block(i));
-  };
-
-  return_op.Tvmult_add = [ops](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      ops[i].Tvmult_add(v.block(i), u.block(i));
-  };
-
-  return return_op;
+  typedef typename BlockLinearOperator<Range, Domain>::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 (size_t i = 0; i < m; ++i)
+    for (size_t 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(new_ops);
 }
 
 
+
 /**
  * @relates LinearOperator
  *
@@ -246,91 +560,29 @@ block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockTyp
  *
  * @ingroup LAOperators
  */
-template <unsigned int m,
-          typename Range = BlockVector<double>,
-          typename Domain = Range>
-LinearOperator<Range, Domain>
+template <size_t m, typename Range, typename Domain>
+BlockLinearOperator<Range, Domain>
 block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType> &op)
 {
-
   static_assert(m > 0,
                 "a blockdiagonal LinearOperator must consist of at least "
                 "one block");
 
-  LinearOperator<Range, Domain> return_op;
-
-  return_op.reinit_range_vector = [op](Range &v, bool fast)
-  {
-    // 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.reinit_range_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.reinit_domain_vector = [op](Domain &v, bool fast)
-  {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(m);
-
-    // And reinitialize every individual block with reinit_domain_vectors:
-    for (unsigned int i = 0; i < m; ++i)
-      op.reinit_domain_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
-
-  return_op.vmult = [op](Range &v, const Domain &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      op.vmult(v.block(i), u.block(i));
-  };
-
-  return_op.vmult_add = [op](Range &v, const Domain &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      op.vmult_add(v.block(i), u.block(i));
-  };
-
-  return_op.Tvmult = [op](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
+  typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+  std::array<BlockType, m> new_ops;
+  new_ops.fill(op);
 
-    for (unsigned int i = 0; i < m; ++i)
-      op.Tvmult(v.block(i), u.block(i));
-  };
-
-  return_op.Tvmult_add = [op](Domain &v, const Range &u)
-  {
-    Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
-    Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
-
-    for (unsigned int i = 0; i < m; ++i)
-      op.Tvmult_add(v.block(i), u.block(i));
-  };
-
-  return return_op;
+  return block_diagonal_operator(new_ops);
 }
 
-
 /**
  * @relates LinearOperator
  *
- * This function implements a forward substitution argument to invert a lower
- * block triangular matrix.
- * As argument, it takes an array of array of LinearOperators @p block_matrix
- * representing a lower block triangular matrix and an array of LinearOperators
- * @p inverse_diagonal representing inverses of diagonal blocks of @p block_matrix.
+ * This function implements forward substitution to invert a lower block
+ * triangular matrix. As arguments, it takes a BlockLinearOperator
+ * @p block_operator representing a block lower triangular matrix, as well as
+ * a BlockLinearOperator @p diagonal_inverse representing inverses of
+ * diagonal blocks of @p block_operator.
  *
  * Let us assume we have a linear system with the following block structure:
  *
@@ -345,85 +597,79 @@ block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename
  * and therefore:
  *    xn = Ann^-1 ( yn - A0n x0 - ... - A(n-1)n x(n-1) )
  *
- * @note Notice that we are not using the whole matrix:
- *       just the lower triangular block matrix obtained
- *       from @p block_matrix it is used.
+ * @note We are not using all blocks of the BlockLinearOperator arguments:
+ * Just the lower triangular block matrix of @p block_operator is used as
+ * well as the diagonal of @p diagonal_inverse.
  *
  * @ingroup LAOperators
  */
-template < unsigned int n,
-           typename Range = BlockVector<double>>
-LinearOperator<Range, Range>
-block_forward_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &,
-                           const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &);
-
-template < unsigned int n,
-           typename Range>
-LinearOperator<Range, Range>
-block_forward_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &block_matrix,
-                           const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &inverse_diagonal)
+template <typename Range, typename Domain>
+LinearOperator<Domain, Range>
+block_forward_substitution(const BlockLinearOperator<Range, Domain> &block_operator,
+                           const BlockLinearOperator<Domain, Range> &diagonal_inverse)
 {
   LinearOperator<Range, Range> return_op;
 
-  return_op.reinit_range_vector = [inverse_diagonal](Range &v, bool fast)
-  {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(n);
-
-    // And reinitialize every individual block with reinit_range_vectors:
-    for (unsigned int i = 0; i < n; ++i)
-      inverse_diagonal[i].reinit_range_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
+  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
+  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
 
-  return_op.reinit_domain_vector = [inverse_diagonal](Range &v, bool fast)
+  return_op.vmult = [block_operator, diagonal_inverse](Range &v, const Range &u)
   {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(n);
-
-    // And reinitialize every individual block with reinit_domain_vectors:
-    for (unsigned int i = 0; i < n; ++i)
-      inverse_diagonal[i].reinit_domain_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
+    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));
 
-  return_op.vmult = [block_matrix, inverse_diagonal](Range &v, const Range &u)
-  {
-    static GrowingVectorMemory<typename  Range::BlockType> vector_memory;
-    typename  Range::BlockType *tmp = vector_memory.alloc();
+    if (m == 0)
+      return;
 
-    inverse_diagonal[0].vmult(v.block(0), u.block(0));
-    for (unsigned int i=1; i<n; ++i)
+    diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
+    for (unsigned int i = 1; i < m; ++i)
       {
-        inverse_diagonal[i].reinit_range_vector(*tmp, /*bool fast=*/ true);
-        *tmp = u.block(i);
-        *tmp *= -1.;
-        for (unsigned int j=0; j<i; ++j)
-          block_matrix[i][j].vmult_add(*tmp, v.block(j));
-        *tmp *= -1.;
-        inverse_diagonal[i].vmult(v.block(i),*tmp);
+        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
       }
-
-    vector_memory.free(tmp);
   };
 
-  return_op.vmult_add = [block_matrix, inverse_diagonal](Range &v, const Range &u)
+  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;
+
     static GrowingVectorMemory<typename  Range::BlockType> vector_memory;
-    typename  Range::BlockType *tmp = vector_memory.alloc();
+    typename Range::BlockType *tmp = vector_memory.alloc();
+
+    diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0));
 
-    inverse_diagonal[0].vmult_add(v.block(0), u.block(0));
-    for (unsigned int i=1; i<n; ++i)
+    for (unsigned int i = 1; i < m; ++i)
       {
-        inverse_diagonal[i].reinit_range_vector(*tmp, /*bool fast=*/ true);
+        diagonal_inverse.block(i, i).reinit_range_vector(*tmp, /*bool fast=*/ true);
         *tmp = u.block(i);
         *tmp *= -1.;
-        for (unsigned int j=0; j<i; ++j)
-          block_matrix[i][j].vmult_add(*tmp, v.block(j));
+        for (unsigned int j = 0; j < i; ++j)
+          block_operator.block(i, j).vmult_add(*tmp, v.block(j));
         *tmp *= -1.;
-        inverse_diagonal[i].vmult_add(v.block(i),*tmp);
+        diagonal_inverse.block(i, i).vmult_add(v.block(i),*tmp);
       }
 
     vector_memory.free(tmp);
@@ -435,11 +681,11 @@ block_forward_substitution(const std::array<std::array<LinearOperator<typename R
 /**
  * @relates LinearOperator
  *
- * This function implements a back substitution argument to invert an upper
- * block triangular matrix.
- * As argument, it takes an array of array of LinearOperators @p block_matrix
- * representing an upper block triangular matrix and an array of LinearOperators
- * @p inverse_diagonal representing inverses of diagonal blocks of @p block_matrix.
+ * This function implements back substitution to invert an upper block
+ * triangular matrix. As arguments, it takes a BlockLinearOperator
+ * @p block_operator representing an upper block triangular matrix, as well
+ * as a BlockLinearOperator @p diagonal_inverse representing inverses of
+ * diagonal blocks of @p block_operator.
  *
  * Let us assume we have a linear system with the following block structure:
  *
@@ -454,85 +700,80 @@ block_forward_substitution(const std::array<std::array<LinearOperator<typename R
  * and therefore:
  *    x0 = A00^-1 ( y0 - A0n xn - ... - A01 x1 )
  *
- * @note Notice that we are not using the whole matrix:
- *       just the upper triangular block matrix obtained
- *       from @p block_matrix it is used.
+ * @note We are not using all blocks of the BlockLinearOperator arguments:
+ * Just the upper triangular block matrix of @p block_operator is used as
+ * well as the diagonal of @p diagonal_inverse.
+
  *
  * @ingroup LAOperators
  */
-template < unsigned int n,
-           typename Range = BlockVector<double>>
-LinearOperator<Range, Range>
-block_back_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &,
-                        const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &);
-
-template < unsigned int n,
-           typename Range>
-LinearOperator<Range, Range>
-block_back_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &block_matrix,
-                        const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &inverse_diagonal)
+template <typename Range, typename Domain>
+LinearOperator<Domain, Range>
+block_back_substitution(const BlockLinearOperator<Range, Domain> &block_operator,
+                        const BlockLinearOperator<Domain, Range> &diagonal_inverse)
 {
   LinearOperator<Range, Range> return_op;
 
-  return_op.reinit_range_vector = [inverse_diagonal](Range &v, bool fast)
-  {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(n);
-
-    // And reinitialize every individual block with reinit_range_vectors:
-    for (unsigned int i = 0; i < n; ++i)
-      inverse_diagonal[i].reinit_range_vector(v.block(i), fast);
-
-    v.collect_sizes();
-  };
+  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
+  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
 
-  return_op.reinit_domain_vector = [inverse_diagonal](Range &v, bool fast)
+  return_op.vmult = [block_operator, diagonal_inverse](Range &v, const Range &u)
   {
-    // Reinitialize the block vector to m blocks:
-    v.reinit(n);
-
-    // And reinitialize every individual block with reinit_domain_vectors:
-    for (unsigned int i = 0; i < n; ++i)
-      inverse_diagonal[i].reinit_domain_vector(v.block(i), fast);
+    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));
 
-    v.collect_sizes();
-  };
+    if (m == 0)
+      return;
 
-  return_op.vmult = [block_matrix, inverse_diagonal](Range &v, const Range &u)
-  {
-    static GrowingVectorMemory<typename  Range::BlockType> vector_memory;
-    typename  Range::BlockType *tmp = vector_memory.alloc();
+    diagonal_inverse.block(m-1, m-1).vmult(v.block(m-1),u.block(m-1));
 
-    inverse_diagonal[n-1].vmult(v.block(n-1),u.block(n-1));
-    for (int i=n-2; i>=0; --i)
+    for (int i = m - 2; i >= 0; --i)
       {
-        inverse_diagonal[i].reinit_range_vector(*tmp, /*bool fast=*/ true);
-        *tmp = u.block(i);
-        *tmp *= -1.;
-        for (int j=i+1; j<n; ++j)
-          block_matrix[i][j].vmult_add(*tmp,v.block(j));
-        *tmp *= -1.;
-        inverse_diagonal[i].vmult(v.block(i),*tmp);
+        auto &dst = v.block(i);
+        dst = u.block(i);
+        dst *= -1.;
+        for (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
       }
-
-    vector_memory.free(tmp);
   };
 
-  return_op.vmult_add = [&block_matrix, inverse_diagonal](Range &v, const Range &u)
+  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));
     static GrowingVectorMemory<typename  Range::BlockType> vector_memory;
     typename  Range::BlockType *tmp = vector_memory.alloc();
 
-    inverse_diagonal[n-1].vmult_add(v.block(n-1),u.block(n-1));
-    for (int i=n-2; i>=0; --i)
+    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)
       {
-        inverse_diagonal[i].reinit_range_vector(*tmp, /*bool fast=*/ true);
+        diagonal_inverse.block(i, i).reinit_range_vector(*tmp, /*bool fast=*/ true);
         *tmp = u.block(i);
         *tmp *= -1.;
-        for (int j=i+1; j<n; ++j)
-          block_matrix[i][j].vmult_add(*tmp,v.block(j));
+        for (int j = i + 1; j < m; ++j)
+          block_operator.block(i, j).vmult_add(*tmp,v.block(j));
         *tmp *= -1.;
-        inverse_diagonal[i].vmult_add(v.block(i),*tmp);
+        diagonal_inverse.block(i, i).vmult_add(v.block(i),*tmp);
       }
 
     vector_memory.free(tmp);
@@ -541,56 +782,6 @@ block_back_substitution(const std::array<std::array<LinearOperator<typename Rang
   return return_op;
 }
 
-/**
- * @relates LinearOperator
- *
- * This function uses above functions block_back_substitution and
- * block_forward_substitution to invert triangular matrices. It takes as input
- * a triangular block matrix @p block_matrix, an array of LinearOperators @p
- * inverse_diagonal representing inverses of block_matrix, and an optional
- * bool @p lower used to specify if block_matrix should be considered as lower
- * triangular matrix (true) or as upper triangular matrix (false). @p lower is
- * equal to true by default.
- *
- */
-
-// 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
-
-template <unsigned int n,
-          typename BlockMatrix,
-          typename Range = BlockVector<double>>
-LinearOperator<Range, Range>
-block_triangular_inverse(const BlockMatrix &,
-                         const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &,
-                         bool lower = true);
-
-template <unsigned int n,
-          typename BlockMatrix,
-          typename Range>
-LinearOperator<Range, Range>
-block_triangular_inverse(const BlockMatrix &block_matrix,
-                         const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &inverse_diagonal,
-                         bool lower)
-{
-  Assert(block_matrix.n_block_rows() == n,
-         ExcDimensionMismatch(block_matrix.n_block_rows(), n));
-  Assert(block_matrix.n_block_rows() == block_matrix.n_block_cols(),
-         ExcDimensionMismatch(block_matrix.n_block_rows(),
-                              block_matrix.n_block_cols()));
-
-  std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> M;
-  for (unsigned int i = 0; i<n; ++i)
-    for (unsigned int j = 0; j<n; ++j)
-      M[i][j] = linear_operator<typename Range::BlockType, typename Range::BlockType>(block_matrix.block(i,j));
-
-  if (lower)
-    return block_forward_substitution<n, Range>(M, inverse_diagonal);
-  else
-    return block_back_substitution<n, Range>(M, inverse_diagonal);
-}
-
 //@}
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/lac/block_linear_operator_02.cc b/tests/lac/block_linear_operator_02.cc
index 89b2827819..56b9a614e9 100644
--- a/tests/lac/block_linear_operator_02.cc
+++ b/tests/lac/block_linear_operator_02.cc
@@ -78,7 +78,7 @@ int main()
     a.block(2, 2).set(i, i, 3.);
 
 
-  auto op_a = linear_operator<BlockVector<double>>(a);
+  auto op_a = block_operator(a);
 
   auto op_b0 = linear_operator(a.block(0, 0));
   auto op_b1 = linear_operator(a.block(1, 1));
@@ -88,68 +88,71 @@ int main()
   std::array<decltype(op_b0), 3> temp{op_b0, op_b1, op_b2};
   auto op_b = block_diagonal_operator<3, BlockVector<double>>(temp);
 
-  // vmult:
 
-  BlockVector<double> u;
-  op_a.reinit_domain_vector(u, false);
-  for (unsigned int i = 0; i < u.size(); ++i) {
-    u[i] = (double)(i+1);
-  }
+  {
+
+    BlockVector<double> u;
+    op_a.reinit_domain_vector(u, false);
+    for (unsigned int i = 0; i < u.size(); ++i) {
+      u[i] = (double)(i+1);
+    }
 
-  PRINTME("u", u);
+    PRINTME("u", u);
 
-  BlockVector<double> v;
-  op_a.reinit_domain_vector(v, false);
-  BlockVector<double> w;
-  op_a.reinit_domain_vector(w, false);
-  BlockVector<double> x;
-  op_a.reinit_domain_vector(x, false);
+    BlockVector<double> v;
+    op_a.reinit_domain_vector(v, false);
+    BlockVector<double> w;
+    op_a.reinit_domain_vector(w, false);
+    BlockVector<double> x;
+    op_a.reinit_domain_vector(x, false);
 
-  op_a.vmult(v, u);
-  PRINTME("Au", v);
+    op_a.vmult(v, u);
+    PRINTME("Au", v);
 
-  op_b.vmult(w, u);
-  PRINTME("Bu", w);
+    op_b.vmult(w, u);
+    PRINTME("Bu", w);
 
-  x = v;
-  x -= w;
-  PRINTME("Au-Bu", x);
+    x = v;
+    x -= w;
+    PRINTME("Au-Bu", x);
 
-  // Test that both objects give the same results:
+    // Test that both objects give the same results:
 
-  auto op_x = op_a - op_b;
+    auto op_x = op_a - op_b;
 
-  op_x.vmult(x, u);
-  PRINTME("(A-B).vmult", x);
+    op_x.vmult(x, u);
+    PRINTME("(A-B).vmult", x);
 
-  x = 0.;
-  op_x.vmult_add(x, u);
-  PRINTME("(A-B).vmult_add", x);
+    x = 0.;
+    op_x.vmult_add(x, u);
+    PRINTME("(A-B).vmult_add", x);
 
-  op_x.Tvmult(x, u);
-  PRINTME("(A-B).Tvmult", x);
+    op_x.Tvmult(x, u);
+    PRINTME("(A-B).Tvmult", x);
 
-  x = 0.;
-  op_x.Tvmult_add(x, u);
-  PRINTME("(A-B).Tvmult_add", x);
+    x = 0.;
+    op_x.Tvmult_add(x, u);
+    PRINTME("(A-B).Tvmult_add", x);
 
 
-  // Test vector reinitalization:
+    // Test vector reinitalization:
 
-  op_x = op_b * op_b * op_b;
-  op_x.vmult(x, u);
-  PRINTME("(B*B*B) vmult", x);
+    op_x = op_b * op_b * op_b;
+    op_x.vmult(x, u);
+    PRINTME("(B*B*B) vmult", x);
 
-  x = 0.;
-  op_x.vmult_add(x, u);
-  PRINTME("(B*B*B) vmult_add", x);
+    x = 0.;
+    op_x.vmult_add(x, u);
+    PRINTME("(B*B*B) vmult_add", x);
 
-  op_x.Tvmult(x, u);
-  PRINTME("(B*B*B) Tvmult", x);
+    op_x.Tvmult(x, u);
+    PRINTME("(B*B*B) Tvmult", x);
 
-  x = 0.;
-  op_x.Tvmult_add(x, u);
-  PRINTME("(B*B*B) Tvmult_add", x);
+    x = 0.;
+    op_x.Tvmult_add(x, u);
+    PRINTME("(B*B*B) Tvmult_add", x);
+
+  }
 
   // And finally the other block_diagonal_operator variant:
 
@@ -158,30 +161,32 @@ int main()
 
   auto op_d = block_diagonal_operator<5, BlockVector<double>>(op_b0);
 
-  op_c.reinit_domain_vector(u, false);
-  for (unsigned int i = 0; i < u.size(); ++i) {
-    u[i] = (double)(i+1);
-  }
-  PRINTME("u", u);
-  op_c.reinit_domain_vector(x, false);
-
-  op_x = op_c - op_d;
-
-  op_x.vmult(x, u);
-  PRINTME("(C-D) vmult", x);
+  {
+    BlockVector<double> u;
+    op_c.reinit_domain_vector(u, false);
+    for (unsigned int i = 0; i < u.size(); ++i) {
+      u[i] = (double)(i+1);
+    }
+    PRINTME("u", u);
 
-  x = 0.;
-  op_x.vmult_add(x, u);
-  PRINTME("(C-D) vmult_add", x);
+    BlockVector<double> x;
+    op_c.reinit_range_vector(x, false);
 
-  op_x.Tvmult(x, u);
-  PRINTME("(C-D) Tvmult", x);
+    auto op_x = op_c - op_d;
 
-  x = 0.;
-  op_x.Tvmult_add(x, u);
-  PRINTME("(C-D) Tvmult_add", x);
+    op_x.vmult(x, u);
+    PRINTME("(C-D) vmult", x);
 
-}
+    x = 0.;
+    op_x.vmult_add(x, u);
+    PRINTME("(C-D) vmult_add", x);
 
+    op_x.Tvmult(x, u);
+    PRINTME("(C-D) Tvmult", x);
 
+    x = 0.;
+    op_x.Tvmult_add(x, u);
+    PRINTME("(C-D) Tvmult_add", x);
+  }
 
+}
diff --git a/tests/lac/block_linear_operator_03.cc b/tests/lac/block_linear_operator_03.cc
index 7c8378cb2a..4ec6e5ecf8 100644
--- a/tests/lac/block_linear_operator_03.cc
+++ b/tests/lac/block_linear_operator_03.cc
@@ -13,237 +13,55 @@
 //
 // ---------------------------------------------------------------------
 
-// Test block_back_substitution and block_forward_substitution:
+// Test block_operator's constructors and copy assignment variants
 
 #include "../tests.h"
 
-#include <deal.II/lac/block_linear_operator.h>
-#include <deal.II/lac/block_sparse_matrix.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/vector.h>
-
-#define PRINTME(name, var) \
-  deallog << "Block vector: " name << ":" << std::endl; \
-  for (unsigned int i = 0; i < var.n_blocks(); ++i) \
-    deallog << "[block " << i << " ]  " << var.block(i);
-
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_linear_operator.h>
 
 using namespace dealii;
 
 int main()
 {
   initlog();
-  deallog << std::setprecision(12);
 
-  // BlockSparseMatrix:
   {
-    BlockDynamicSparsityPattern dsp(3, 3);
-    for (unsigned int i = 0; i < 3; ++i)
-      for (unsigned int j = 0; j < 3; ++j)
-        dsp.block(i, j).reinit (1, 1);
-    dsp.collect_sizes ();
-
-    BlockSparsityPattern sparsity_pattern;
-    sparsity_pattern.copy_from(dsp);
-    sparsity_pattern.compress();
-
-    BlockSparseMatrix<double> a (sparsity_pattern);
-    for (unsigned int i = 0; i < 3; ++i)
-    {
-      a.block(i,i).set(0, 0, i+i +1);
-      for (unsigned int j = 0; j < i; ++j)
-        a.block(i,j).set(0, 0, 10);
-    }
-
-    BlockSparseMatrix<double> d(sparsity_pattern);
-    for (unsigned int i = 0; i < 3; ++i)
-        d.block(i,i).set(0,0, 1.0 / (i+i +1) );
-
-    auto op_a         = linear_operator< BlockVector<double> >(a);
-
-    auto a00 = linear_operator< Vector<double>, Vector<double> >(a.block(0,0));
-    auto a01 = linear_operator< Vector<double>, Vector<double> >(a.block(0,1));
-    auto a02 = linear_operator< Vector<double>, Vector<double> >(a.block(0,2));
-    auto a10 = linear_operator< Vector<double>, Vector<double> >(a.block(1,0));
-    auto a11 = linear_operator< Vector<double>, Vector<double> >(a.block(1,1));
-    auto a12 = linear_operator< Vector<double>, Vector<double> >(a.block(1,2));
-    auto a20 = linear_operator< Vector<double>, Vector<double> >(a.block(2,0));
-    auto a21 = linear_operator< Vector<double>, Vector<double> >(a.block(2,1));
-    auto a22 = linear_operator< Vector<double>, Vector<double> >(a.block(2,2));
-
-    auto d00 = linear_operator< Vector<double>, Vector<double> >(d.block(0,0));
-    auto d11 = linear_operator< Vector<double>, Vector<double> >(d.block(1,1));
-    auto d22 = linear_operator< Vector<double>, Vector<double> >(d.block(2,2));
-
-    auto inverse_op_a = block_forward_substitution< 3, BlockVector<double> >(
-      {{
-        {{a00, a01, a02}},
-        {{a10, a11, a12}},
-        {{a20, a21, a22}}
-      }},
-      { {d00, d11, d22}});
-
-    auto identity = inverse_op_a * op_a;
-
-    BlockVector<double> u;
-    BlockVector<double> v;
-
-    deallog << " -- Matrix -- " << std::endl;
-    op_a.reinit_domain_vector(u, false);
-    op_a.reinit_range_vector(v, false);
-    for(unsigned int j = 0; j<3; ++j)
-    {
-      for(unsigned int i = 0; i<3; ++i)
-      {
-        u.block(i)[0] = 0;;
-        v.block(i)[0] = 0;
-      }
-      u.block(j)[0] = 1;
-
-      op_a.vmult(v, u);
-
-      PRINTME("v", v);
-    }
-
-    deallog << " -- Inverse -- " << std::endl;
-    inverse_op_a.reinit_domain_vector(u, false);
-    inverse_op_a.reinit_range_vector(v, true);
-    for(unsigned int j = 0; j<3; ++j)
-    {
-      for(unsigned int i = 0; i<3; ++i)
-      {
-        u.block(i)[0] = 0;
-        v.block(i)[0] = 0;
-      }
-      u.block(j)[0] = 1;;
-
-      inverse_op_a.vmult(v, u);
-
-      PRINTME("v", v);
-    }
-
-    deallog << " -- Identity -- " << std::endl;
-    identity.reinit_domain_vector(u, false);
-    identity.reinit_range_vector(v, false);
-    for(unsigned int j = 0; j<3; ++j)
-    {
-      for(unsigned int i = 0; i<3; ++i)
-      {
-        u.block(i)[0] = 0;;
-        v.block(i)[0] = 0;
-      }
-      u.block(j)[0] = 1;;
-
-      identity.vmult(v, u);
-
-      PRINTME("v", v);
-    }
+    BlockSparseMatrix<double> a;
+
+    BlockLinearOperator<> op_b(a);
+    auto op_c = block_operator<BlockVector<double>>(a);
+    op_c = a;
+
+    auto op_d = block_operator(a);
+    op_d = a;
   }
 
+  auto op_a = LinearOperator<>();
 
   {
-    BlockDynamicSparsityPattern dsp(3, 3);
-    for (unsigned int i = 0; i < 3; ++i)
-      for (unsigned int j = 0; j < 3; ++j)
-        dsp.block(i, j).reinit (1, 1);
-    dsp.collect_sizes ();
-
-    BlockSparsityPattern sparsity_pattern;
-    sparsity_pattern.copy_from(dsp);
-    sparsity_pattern.compress();
-
-    BlockSparseMatrix<double> a (sparsity_pattern);
-    for (unsigned int i = 0; i < 3; ++i)
-    {
-      a.block(i,i).set(0, 0, i+i +1);
-      for (unsigned int j = i+1; j < 3; ++j)
-        a.block(i,j).set(0, 0, 10);
-    }
-
-    BlockSparseMatrix<double> d(sparsity_pattern);
-    for (unsigned int i = 0; i < 3; ++i)
-        d.block(i,i).set(0,0, 1.0 / (i+i +1) );
-
-        auto op_a         = linear_operator< BlockVector<double> >(a);
-
-        auto a00 = linear_operator< Vector<double>, Vector<double> >(a.block(0,0));
-        auto a01 = linear_operator< Vector<double>, Vector<double> >(a.block(0,1));
-        auto a02 = linear_operator< Vector<double>, Vector<double> >(a.block(0,2));
-        auto a10 = linear_operator< Vector<double>, Vector<double> >(a.block(1,0));
-        auto a11 = linear_operator< Vector<double>, Vector<double> >(a.block(1,1));
-        auto a12 = linear_operator< Vector<double>, Vector<double> >(a.block(1,2));
-        auto a20 = linear_operator< Vector<double>, Vector<double> >(a.block(2,0));
-        auto a21 = linear_operator< Vector<double>, Vector<double> >(a.block(2,1));
-        auto a22 = linear_operator< Vector<double>, Vector<double> >(a.block(2,2));
-
-        auto d00 = linear_operator< Vector<double>, Vector<double> >(d.block(0,0));
-        auto d11 = linear_operator< Vector<double>, Vector<double> >(d.block(1,1));
-        auto d22 = linear_operator< Vector<double>, Vector<double> >(d.block(2,2));
-
-        auto inverse_op_a = block_back_substitution< 3, BlockVector<double> >(
-          {{
-            {{a00, a01, a02}},
-            {{a10, a11, a12}},
-            {{a20, a21, a22}}
-          }},
-          { {d00, d11, d22}});
-
-        auto identity = inverse_op_a * op_a;
-
-    BlockVector<double> u;
-    BlockVector<double> v;
-
-    deallog << " -- Matrix -- " << std::endl;
-    op_a.reinit_domain_vector(u, false);
-    op_a.reinit_range_vector(v, false);
-    for(unsigned int j = 0; j<3; ++j)
-    {
-      for(unsigned int i = 0; i<3; ++i)
-      {
-        u.block(i)[0] = 0;;
-        v.block(i)[0] = 0;
-      }
-      u.block(j)[0] = 1;
-
-      op_a.vmult(v, u);
-
-      PRINTME("v", v);
-    }
-
-    deallog << " -- Inverse -- " << std::endl;
-    inverse_op_a.reinit_domain_vector(u, false);
-    inverse_op_a.reinit_range_vector(v, true);
-    for(unsigned int j = 0; j<3; ++j)
-    {
-      for(unsigned int i = 0; i<3; ++i)
-      {
-        u.block(i)[0] = 0;
-        v.block(i)[0] = 0;
-      }
-      u.block(j)[0] = 1;;
-
-      inverse_op_a.vmult(v, u);
-
-      PRINTME("v", v);
-    }
-
-    deallog << " -- Identity -- " << std::endl;
-    identity.reinit_domain_vector(u, false);
-    identity.reinit_range_vector(v, false);
-    for(unsigned int j = 0; j<3; ++j)
-    {
-      for(unsigned int i = 0; i<3; ++i)
-      {
-        u.block(i)[0] = 0;;
-        v.block(i)[0] = 0;
-      }
-      u.block(j)[0] = 1;;
-
-      identity.vmult(v, u);
-
-      PRINTME("v", v);
-    }
+    std::array<std::array<decltype(op_a), 2>, 2> temp {op_a, op_a, op_a, op_a};
+
+    BlockLinearOperator<> op_b(temp);
+    auto op_c = block_operator<2, 2, BlockVector<double>>(temp);
+    op_c = temp;
+
+    auto op_d = block_operator(temp);
+    op_d = temp;
   }
+
+  {
+    std::array<decltype(op_a), 2> temp {op_a, op_a};
+
+    auto op_c = block_diagonal_operator<2, BlockVector<double>>(temp);
+    op_c = temp;
+
+    auto op_d = block_diagonal_operator(temp);
+    op_d = temp;
+  }
+
+  deallog << "OK!" << std::endl;
 }
+
diff --git a/tests/lac/block_linear_operator_03.with_cxx11=on.output b/tests/lac/block_linear_operator_03.with_cxx11=on.output
index 1e80e98cd0..5cfb783b8f 100644
--- a/tests/lac/block_linear_operator_03.with_cxx11=on.output
+++ b/tests/lac/block_linear_operator_03.with_cxx11=on.output
@@ -1,79 +1,2 @@
 
-DEAL:: -- Matrix -- 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  1.00000000000 
-DEAL::[block 1 ]  10.0000000000 
-DEAL::[block 2 ]  10.0000000000 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  0.00000000000 
-DEAL::[block 1 ]  3.00000000000 
-DEAL::[block 2 ]  10.0000000000 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  0.00000000000 
-DEAL::[block 1 ]  0.00000000000 
-DEAL::[block 2 ]  5.00000000000 
-DEAL:: -- Inverse -- 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  1.00000000000 
-DEAL::[block 1 ]  -3.33333333333 
-DEAL::[block 2 ]  4.66666666667 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  0.00000000000 
-DEAL::[block 1 ]  0.333333333333 
-DEAL::[block 2 ]  -0.666666666667 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  0.00000000000 
-DEAL::[block 1 ]  0.00000000000 
-DEAL::[block 2 ]  0.200000000000 
-DEAL:: -- Identity -- 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  1.00000000000 
-DEAL::[block 1 ]  0.00000000000 
-DEAL::[block 2 ]  0.00000000000 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  0.00000000000 
-DEAL::[block 1 ]  1.00000000000 
-DEAL::[block 2 ]  0.00000000000 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  0.00000000000 
-DEAL::[block 1 ]  0.00000000000 
-DEAL::[block 2 ]  1.00000000000 
-DEAL:: -- Matrix -- 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  1.00000000000 
-DEAL::[block 1 ]  0.00000000000 
-DEAL::[block 2 ]  0.00000000000 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  10.0000000000 
-DEAL::[block 1 ]  3.00000000000 
-DEAL::[block 2 ]  0.00000000000 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  10.0000000000 
-DEAL::[block 1 ]  10.0000000000 
-DEAL::[block 2 ]  5.00000000000 
-DEAL:: -- Inverse -- 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  1.00000000000 
-DEAL::[block 1 ]  0.00000000000 
-DEAL::[block 2 ]  0.00000000000 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  -3.33333333333 
-DEAL::[block 1 ]  0.333333333333 
-DEAL::[block 2 ]  0.00000000000 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  4.66666666667 
-DEAL::[block 1 ]  -0.666666666667 
-DEAL::[block 2 ]  0.200000000000 
-DEAL:: -- Identity -- 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  1.00000000000 
-DEAL::[block 1 ]  0.00000000000 
-DEAL::[block 2 ]  0.00000000000 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  0.00000000000 
-DEAL::[block 1 ]  1.00000000000 
-DEAL::[block 2 ]  0.00000000000 
-DEAL::Block vector: v:
-DEAL::[block 0 ]  0.00000000000 
-DEAL::[block 1 ]  0.00000000000 
-DEAL::[block 2 ]  1.00000000000 
+DEAL::OK!
diff --git a/tests/lac/block_linear_operator_04.cc b/tests/lac/block_linear_operator_04.cc
index e49548d508..07cbca9e8a 100644
--- a/tests/lac/block_linear_operator_04.cc
+++ b/tests/lac/block_linear_operator_04.cc
@@ -61,12 +61,14 @@ int main()
     for (unsigned int i = 0; i < 3; ++i)
         d.block(i,i).set(0,0, 1.0 / (i+i +1) );
 
-    auto op_a         = linear_operator< BlockVector<double> >(a);
-    auto d00 = linear_operator< Vector<double>, Vector<double> >(d.block(0,0));
-    auto d11 = linear_operator< Vector<double>, Vector<double> >(d.block(1,1));
-    auto d22 = linear_operator< Vector<double>, Vector<double> >(d.block(2,2));
+    auto op_a  = linear_operator<BlockVector<double>>(a);
+
+    auto op_b1 = block_operator(a);
+    auto op_b2 = block_diagonal_operator(d);
+
+    auto inverse_op_a =
+        block_forward_substitution<BlockVector<double>>(op_b1, op_b2);
 
-    auto inverse_op_a = block_triangular_inverse<3, BlockSparseMatrix<double > >(a, {{d00, d11, d22}});
     auto identity = inverse_op_a * op_a;
 
     BlockVector<double> u;
@@ -146,14 +148,16 @@ int main()
 
     BlockSparseMatrix<double> d(sparsity_pattern);
     for (unsigned int i = 0; i < 3; ++i)
-        d.block(i,i).set(0,0, 1.0 / (i+i +1) );
+      d.block(i, i).set(0, 0, 1.0 / (i + i + 1));
+
+    auto op_a = linear_operator<BlockVector<double>>(a);
+
+    auto op_b1 = block_operator(a);
+    auto op_b2 = block_diagonal_operator(d);
 
-    auto op_a         = linear_operator< BlockVector<double> >(a);
-    auto d00 = linear_operator< Vector<double>, Vector<double> >(d.block(0,0));
-    auto d11 = linear_operator< Vector<double>, Vector<double> >(d.block(1,1));
-    auto d22 = linear_operator< Vector<double>, Vector<double> >(d.block(2,2));
+    auto inverse_op_a =
+        block_back_substitution<BlockVector<double>>(op_b1, op_b2);
 
-    auto inverse_op_a = block_triangular_inverse< 3, BlockSparseMatrix<double > >(a, {{d00, d11, d22}}, false);
     auto identity = inverse_op_a * op_a;
 
     BlockVector<double> u;