From 908eee11cede69531e9513e7d900d110bd62f43c Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Thu, 3 May 2018 15:55:20 -0500 Subject: [PATCH] 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. --- include/deal.II/lac/block_linear_operator.h | 40 +++++++++++++++++++++ 1 file changed, 40 insertions(+) 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)); -- 2.39.5