BlockLinearOperator: Add assert ensuring that u != v in vmult

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.

diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h
index ebc8ad45cf2702ef1d941be3d9a8bb56130e00c6..3ed6ac2ff0e48be00c584b744b7f48b5a8fe23a7 100644 (file)
--- a/include/deal.II/lac/block_linear_operator.h
+++ b/include/deal.II/lac/block_linear_operator.h
@@ -356,6 +356,16 @@ namespace internal
         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));
@@ -371,6 +381,16 @@ namespace internal
         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));
@@ -383,6 +403,16 @@ namespace internal
         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));
@@ -398,6 +428,16 @@ namespace internal
         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));