From: Matthias Maier <tamiko@43-1.org>
Date: Thu, 3 May 2018 20:55:20 +0000 (-0500)
Subject: BlockLinearOperator: Add assert ensuring that u != v in vmult
X-Git-Tag: v9.0.0-rc1~25^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=908eee11cede69531e9513e7d900d110bd62f43c;p=dealii.git

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 ebc8ad45cf..3ed6ac2ff0 100644
--- 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));