BlockLinearOperator: Allow identical destination and source

diff --git a/doc/news/changes/minor/20180610DanielArndt b/doc/news/changes/minor/20180610DanielArndt
new file mode 100644 (file)
index 0000000..9b1a77a
--- /dev/null
+++ b/doc/news/changes/minor/20180610DanielArndt
@@ -0,0 +1,5 @@
+Fixed: BlockLinearOperator::vmult(), BlockLinearOperator::vmult_add(),
+BlockLinearOperator::Tvmult() and BlockLinearOperator::Tvmult_add()
+can be used with identical source and destination vectors.
+<br>
+(Daniel Arndt, 2018/06/10)

diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h
index 574d162d2ddca03431ffd0a0cd1ff9653b484e59..6161d630494b0c7b689b9c7a1d2580212ba741b4 100644 (file)
--- a/include/deal.II/lac/block_linear_operator.h
+++ b/include/deal.II/lac/block_linear_operator.h
@@ -351,6 +351,43 @@ namespace internal
 {
   namespace BlockLinearOperatorImplementation
   {
+    // A helper function to apply a given vmult, or Tvmult to a vector with
+    // intermediate storage, similar to the corresponding helper
+    // function for LinearOperator. Here, two operators are used.
+    // The first one takes care of the first "column" and typically doesn't add.
+    // On the other hand, the second operator is normally an adding one.
+    template <typename Function1,
+              typename Function2,
+              typename Range,
+              typename Domain>
+    void
+    apply_with_intermediate_storage(const Function1 &first_op,
+                                    const Function2 &loop_op,
+                                    Range &          v,
+                                    const Domain &   u,
+                                    bool             add)
+    {
+      GrowingVectorMemory<Range> vector_memory;
+
+      typename VectorMemory<Range>::Pointer tmp(vector_memory);
+      tmp->reinit(v, /*bool omit_zeroing_entries =*/true);
+
+      const unsigned int n = u.n_blocks();
+      const unsigned int m = v.n_blocks();
+
+      for (unsigned int i = 0; i < m; ++i)
+        {
+          first_op(*tmp, u, i, 0);
+          for (unsigned int j = 1; j < n; ++j)
+            loop_op(*tmp, u, i, j);
+        }
+
+      if (add)
+        v += *tmp;
+      else
+        v = *tmp;
+    }
+
     // Populate the LinearOperator interfaces with the help of the
     // BlockLinearOperator functions
     template <typename Range, typename Domain, typename BlockPayload>
@@ -384,96 +421,138 @@ namespace internal
         v.collect_sizes();
       };
 
-      op.vmult = [=](Range &v, const Domain &u) {
+      op.vmult = [&op](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)
+        if (PointerComparison::equal(&v, &u))
           {
-            op.block(i, 0).vmult(v.block(i), u.block(0));
-            for (unsigned int j = 1; j < n; ++j)
+            const auto first_op = [&op](Range &            v,
+                                        const Domain &     u,
+                                        const unsigned int i,
+                                        const unsigned int j) {
+              op.block(i, j).vmult(v.block(i), u.block(j));
+            };
+
+            const auto loop_op = [&op](Range &            v,
+                                       const Domain &     u,
+                                       const unsigned int i,
+                                       const unsigned int j) {
               op.block(i, j).vmult_add(v.block(i), u.block(j));
+            };
+
+            apply_with_intermediate_storage(first_op, loop_op, v, u, false);
+          }
+        else
+          {
+            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) {
+      op.vmult_add = [&op](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"));
+        if (PointerComparison::equal(&v, &u))
+          {
+            const auto first_op = [&op](Range &            v,
+                                        const Domain &     u,
+                                        const unsigned int i,
+                                        const unsigned int j) {
+              op.block(i, j).vmult(v.block(i), u.block(j));
+            };
+
+            const auto loop_op = [&op](Range &            v,
+                                       const Domain &     u,
+                                       const unsigned int i,
+                                       const unsigned int j) {
+              op.block(i, j).vmult_add(v.block(i), u.block(j));
+            };
 
-        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));
+            apply_with_intermediate_storage(first_op, loop_op, v, u, true);
+          }
+        else
+          {
+            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) {
+      op.Tvmult = [&op](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)
+        if (PointerComparison::equal(&v, &u))
           {
-            op.block(0, i).Tvmult(v.block(i), u.block(0));
-            for (unsigned int j = 1; j < m; ++j)
+            const auto first_op = [&op](Range &            v,
+                                        const Domain &     u,
+                                        const unsigned int i,
+                                        const unsigned int j) {
+              op.block(j, i).Tvmult(v.block(i), u.block(j));
+            };
+
+            const auto loop_op = [&op](Range &            v,
+                                       const Domain &     u,
+                                       const unsigned int i,
+                                       const unsigned int j) {
               op.block(j, i).Tvmult_add(v.block(i), u.block(j));
+            };
+
+            apply_with_intermediate_storage(first_op, loop_op, v, u, false);
+          }
+        else
+          {
+            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) {
+      op.Tvmult_add = [&op](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"));
+        if (PointerComparison::equal(&v, &u))
+          {
+            const auto first_op = [&op](Range &            v,
+                                        const Domain &     u,
+                                        const unsigned int i,
+                                        const unsigned int j) {
+              op.block(j, i).Tvmult(v.block(i), u.block(j));
+            };
+
+            const auto loop_op = [&op](Range &            v,
+                                       const Domain &     u,
+                                       const unsigned int i,
+                                       const unsigned int j) {
+              op.block(j, i).Tvmult_add(v.block(i), u.block(j));
+            };
 
-        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));
+            apply_with_intermediate_storage(first_op, loop_op, v, u, true);
+          }
+        else
+          {
+            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));
+          }
       };
     }
 

diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h
index c6fca2275895a6548659dd65ee6b1b9aa9af37c9..fa1299c04e34f32af2e1eacc5983558706053789 100644 (file)
--- a/include/deal.II/lac/linear_operator.h
+++ b/include/deal.II/lac/linear_operator.h
@@ -1172,7 +1172,7 @@ namespace
               [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
               v,
               u,
-              /*bool add =*/false);
+              /*bool add =*/true);
           }
         else
           {

diff --git a/tests/lac/block_linear_operator_06.cc b/tests/lac/block_linear_operator_06.cc
new file mode 100644 (file)
index 0000000..8329a10
--- /dev/null
+++ b/tests/lac/block_linear_operator_06.cc
@@ -0,0 +1,191 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Tests for block_operator for opertaions with identical source destinations
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q_generic.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.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 "../tests.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);
+
+
+using namespace dealii;
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  static const int dim = 2;
+
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation);
+  triangulation.refine_global(2);
+
+  MappingQGeneric<dim> mapping_q1(1);
+  FESystem<dim>        fe(FE_Q<dim>(2), 1, FE_Q<dim>(1), 1);
+  DoFHandler<dim>      dof_handler(triangulation);
+
+  dof_handler.distribute_dofs(fe);
+
+  std::vector<types::global_dof_index> dofs_per_component(2);
+  DoFTools::count_dofs_per_component(dof_handler, dofs_per_component);
+  const unsigned int n_u = dofs_per_component[0], n_p = dofs_per_component[1];
+
+  BlockDynamicSparsityPattern dsp(2, 2);
+  dsp.block(0, 0).reinit(n_u, n_u);
+  dsp.block(1, 0).reinit(n_p, n_u);
+  dsp.block(0, 1).reinit(n_u, n_p);
+  dsp.block(1, 1).reinit(n_p, n_p);
+
+  for (unsigned int i = 0; i < n_u; ++i)
+    {
+      for (unsigned int j = 0; j < n_u; ++j)
+        dsp.block(0, 0).add(i, j);
+      for (unsigned int j = 0; j < n_p; ++j)
+        {
+          dsp.block(0, 1).add(i, j);
+          dsp.block(1, 0).add(j, i);
+        }
+    }
+  for (unsigned int i = 0; i < n_p; ++i)
+    for (unsigned int j = 0; j < n_p; ++j)
+      dsp.block(1, 1).add(i, j);
+  dsp.collect_sizes();
+
+  BlockSparsityPattern sparsity_pattern;
+  sparsity_pattern.copy_from(dsp);
+  sparsity_pattern.compress();
+
+  BlockSparseMatrix<double> a(sparsity_pattern);
+
+  // Come up with some simple, but non-trivial structure:
+
+  for (unsigned int i = 0; i < a.block(0, 0).n(); ++i)
+    a.block(0, 0).set(i, i, 10.);
+
+  for (unsigned int i = 0; i < a.block(1, 1).n(); ++i)
+    a.block(1, 1).set(i, i, 5.);
+
+  for (unsigned int i = 0; i < a.block(1, 0).m(); ++i)
+    a.block(0, 1).set(i, i, 3.);
+
+  for (unsigned int i = 0; i < a.block(0, 1).n(); ++i)
+    a.block(1, 0).set(i, i, 1.);
+
+
+  auto op_a = linear_operator<BlockVector<double>>(a);
+
+  auto op_b00 = linear_operator(a.block(0, 0));
+  auto op_b01 = linear_operator(a.block(0, 1));
+  auto op_b10 = linear_operator(a.block(1, 0));
+  auto op_b11 = linear_operator(a.block(1, 1));
+
+  std::array<std::array<decltype(op_b00), 2>, 2> temp{
+    {{op_b00, op_b01}, {op_b10, op_b11}}};
+  auto op_b = block_operator<2, 2, 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_copy;
+  op_a.reinit_domain_vector(u_copy, false);
+  u_copy = u;
+
+  BlockVector<double> result;
+  op_a.reinit_domain_vector(result, false);
+
+  PRINTME("u", u);
+  PRINTME("u_copy", u);
+  {
+    op_a.vmult(u, u);
+    PRINTME("A.vmult", u);
+    result = u;
+    u      = u_copy;
+
+    op_b.vmult(u, u);
+    PRINTME("B.vmult", u);
+    result -= u;
+    u = u_copy;
+
+    AssertThrow(result.linfty_norm() < 1.e-10, ExcInternalError());
+  }
+
+  {
+    op_a.Tvmult(u, u);
+    PRINTME("A.Tvmult", u);
+    result = u;
+    u      = u_copy;
+
+    op_b.Tvmult(u, u);
+    PRINTME("B.Tvmult", u);
+    result -= u;
+    u = u_copy;
+
+    AssertThrow(result.linfty_norm() < 1.e-10, ExcInternalError());
+  }
+
+  {
+    op_a.vmult_add(u, u);
+    u -= u_copy;
+    PRINTME("A.vmult_add", u);
+    result = u;
+    u      = u_copy;
+
+    op_b.vmult_add(u, u);
+    u -= u_copy;
+    PRINTME("B.vmult_add", u);
+    result -= u;
+    u = u_copy;
+
+    AssertThrow(result.linfty_norm() < 1.e-10, ExcInternalError());
+  }
+
+  {
+    op_a.Tvmult_add(u, u);
+    PRINTME("A.Tvmult_add", u);
+    result = u;
+    u      = u_copy;
+
+    op_b.Tvmult_add(u, u);
+    PRINTME("B.Tvmult_add", u);
+    result -= u;
+
+    AssertThrow(result.linfty_norm() < 1.e-10, ExcInternalError());
+  }
+}

diff --git a/tests/lac/block_linear_operator_06.output b/tests/lac/block_linear_operator_06.output
new file mode 100644 (file)
index 0000000..c3d40d3
--- /dev/null
+++ b/tests/lac/block_linear_operator_06.output
@@ -0,0 +1,31 @@
+
+DEAL::Block vector: u:
+DEAL::[block 0 ]  1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000 51.00000000 52.00000000 53.00000000 54.00000000 55.00000000 56.00000000 57.00000000 58.00000000 59.00000000 60.00000000 61.00000000 62.00000000 63.00000000 64.00000000 65.00000000 66.00000000 67.00000000 68.00000000 69.00000000 70.00000000 71.00000000 72.00000000 73.00000000 74.00000000 75.00000000 76.00000000 77.00000000 78.00000000 79.00000000 80.00000000 81.00000000 
+DEAL::[block 1 ]  82.00000000 83.00000000 84.00000000 85.00000000 86.00000000 87.00000000 88.00000000 89.00000000 90.00000000 91.00000000 92.00000000 93.00000000 94.00000000 95.00000000 96.00000000 97.00000000 98.00000000 99.00000000 100.0000000 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 
+DEAL::Block vector: u_copy:
+DEAL::[block 0 ]  1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000 51.00000000 52.00000000 53.00000000 54.00000000 55.00000000 56.00000000 57.00000000 58.00000000 59.00000000 60.00000000 61.00000000 62.00000000 63.00000000 64.00000000 65.00000000 66.00000000 67.00000000 68.00000000 69.00000000 70.00000000 71.00000000 72.00000000 73.00000000 74.00000000 75.00000000 76.00000000 77.00000000 78.00000000 79.00000000 80.00000000 81.00000000 
+DEAL::[block 1 ]  82.00000000 83.00000000 84.00000000 85.00000000 86.00000000 87.00000000 88.00000000 89.00000000 90.00000000 91.00000000 92.00000000 93.00000000 94.00000000 95.00000000 96.00000000 97.00000000 98.00000000 99.00000000 100.0000000 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 
+DEAL::Block vector: A.vmult:
+DEAL::[block 0 ]  256.0000000 269.0000000 282.0000000 295.0000000 308.0000000 321.0000000 334.0000000 347.0000000 360.0000000 373.0000000 386.0000000 399.0000000 412.0000000 425.0000000 438.0000000 451.0000000 464.0000000 477.0000000 490.0000000 503.0000000 516.0000000 529.0000000 542.0000000 555.0000000 568.0000000 260.0000000 270.0000000 280.0000000 290.0000000 300.0000000 310.0000000 320.0000000 330.0000000 340.0000000 350.0000000 360.0000000 370.0000000 380.0000000 390.0000000 400.0000000 410.0000000 420.0000000 430.0000000 440.0000000 450.0000000 460.0000000 470.0000000 480.0000000 490.0000000 500.0000000 510.0000000 520.0000000 530.0000000 540.0000000 550.0000000 560.0000000 570.0000000 580.0000000 590.0000000 600.0000000 610.0000000 620.0000000 630.0000000 640.0000000 650.0000000 660.0000000 670.0000000 680.0000000 690.0000000 700.0000000 710.0000000 720.0000000 730.0000000 740.0000000 750.0000000 760.0000000 770.0000000 780.0000000 790.0000000 800.0000000 810.0000000 
+DEAL::[block 1 ]  411.0000000 417.0000000 423.0000000 429.0000000 435.0000000 441.0000000 447.0000000 453.0000000 459.0000000 465.0000000 471.0000000 477.0000000 483.0000000 489.0000000 495.0000000 501.0000000 507.0000000 513.0000000 519.0000000 525.0000000 531.0000000 537.0000000 543.0000000 549.0000000 555.0000000 
+DEAL::Block vector: B.vmult:
+DEAL::[block 0 ]  256.0000000 269.0000000 282.0000000 295.0000000 308.0000000 321.0000000 334.0000000 347.0000000 360.0000000 373.0000000 386.0000000 399.0000000 412.0000000 425.0000000 438.0000000 451.0000000 464.0000000 477.0000000 490.0000000 503.0000000 516.0000000 529.0000000 542.0000000 555.0000000 568.0000000 260.0000000 270.0000000 280.0000000 290.0000000 300.0000000 310.0000000 320.0000000 330.0000000 340.0000000 350.0000000 360.0000000 370.0000000 380.0000000 390.0000000 400.0000000 410.0000000 420.0000000 430.0000000 440.0000000 450.0000000 460.0000000 470.0000000 480.0000000 490.0000000 500.0000000 510.0000000 520.0000000 530.0000000 540.0000000 550.0000000 560.0000000 570.0000000 580.0000000 590.0000000 600.0000000 610.0000000 620.0000000 630.0000000 640.0000000 650.0000000 660.0000000 670.0000000 680.0000000 690.0000000 700.0000000 710.0000000 720.0000000 730.0000000 740.0000000 750.0000000 760.0000000 770.0000000 780.0000000 790.0000000 800.0000000 810.0000000 
+DEAL::[block 1 ]  411.0000000 417.0000000 423.0000000 429.0000000 435.0000000 441.0000000 447.0000000 453.0000000 459.0000000 465.0000000 471.0000000 477.0000000 483.0000000 489.0000000 495.0000000 501.0000000 507.0000000 513.0000000 519.0000000 525.0000000 531.0000000 537.0000000 543.0000000 549.0000000 555.0000000 
+DEAL::Block vector: A.Tvmult:
+DEAL::[block 0 ]  92.00000000 103.0000000 114.0000000 125.0000000 136.0000000 147.0000000 158.0000000 169.0000000 180.0000000 191.0000000 202.0000000 213.0000000 224.0000000 235.0000000 246.0000000 257.0000000 268.0000000 279.0000000 290.0000000 301.0000000 312.0000000 323.0000000 334.0000000 345.0000000 356.0000000 260.0000000 270.0000000 280.0000000 290.0000000 300.0000000 310.0000000 320.0000000 330.0000000 340.0000000 350.0000000 360.0000000 370.0000000 380.0000000 390.0000000 400.0000000 410.0000000 420.0000000 430.0000000 440.0000000 450.0000000 460.0000000 470.0000000 480.0000000 490.0000000 500.0000000 510.0000000 520.0000000 530.0000000 540.0000000 550.0000000 560.0000000 570.0000000 580.0000000 590.0000000 600.0000000 610.0000000 620.0000000 630.0000000 640.0000000 650.0000000 660.0000000 670.0000000 680.0000000 690.0000000 700.0000000 710.0000000 720.0000000 730.0000000 740.0000000 750.0000000 760.0000000 770.0000000 780.0000000 790.0000000 800.0000000 810.0000000 
+DEAL::[block 1 ]  413.0000000 421.0000000 429.0000000 437.0000000 445.0000000 453.0000000 461.0000000 469.0000000 477.0000000 485.0000000 493.0000000 501.0000000 509.0000000 517.0000000 525.0000000 533.0000000 541.0000000 549.0000000 557.0000000 565.0000000 573.0000000 581.0000000 589.0000000 597.0000000 605.0000000 
+DEAL::Block vector: B.Tvmult:
+DEAL::[block 0 ]  92.00000000 103.0000000 114.0000000 125.0000000 136.0000000 147.0000000 158.0000000 169.0000000 180.0000000 191.0000000 202.0000000 213.0000000 224.0000000 235.0000000 246.0000000 257.0000000 268.0000000 279.0000000 290.0000000 301.0000000 312.0000000 323.0000000 334.0000000 345.0000000 356.0000000 260.0000000 270.0000000 280.0000000 290.0000000 300.0000000 310.0000000 320.0000000 330.0000000 340.0000000 350.0000000 360.0000000 370.0000000 380.0000000 390.0000000 400.0000000 410.0000000 420.0000000 430.0000000 440.0000000 450.0000000 460.0000000 470.0000000 480.0000000 490.0000000 500.0000000 510.0000000 520.0000000 530.0000000 540.0000000 550.0000000 560.0000000 570.0000000 580.0000000 590.0000000 600.0000000 610.0000000 620.0000000 630.0000000 640.0000000 650.0000000 660.0000000 670.0000000 680.0000000 690.0000000 700.0000000 710.0000000 720.0000000 730.0000000 740.0000000 750.0000000 760.0000000 770.0000000 780.0000000 790.0000000 800.0000000 810.0000000 
+DEAL::[block 1 ]  413.0000000 421.0000000 429.0000000 437.0000000 445.0000000 453.0000000 461.0000000 469.0000000 477.0000000 485.0000000 493.0000000 501.0000000 509.0000000 517.0000000 525.0000000 533.0000000 541.0000000 549.0000000 557.0000000 565.0000000 573.0000000 581.0000000 589.0000000 597.0000000 605.0000000 
+DEAL::Block vector: A.vmult_add:
+DEAL::[block 0 ]  256.0000000 269.0000000 282.0000000 295.0000000 308.0000000 321.0000000 334.0000000 347.0000000 360.0000000 373.0000000 386.0000000 399.0000000 412.0000000 425.0000000 438.0000000 451.0000000 464.0000000 477.0000000 490.0000000 503.0000000 516.0000000 529.0000000 542.0000000 555.0000000 568.0000000 260.0000000 270.0000000 280.0000000 290.0000000 300.0000000 310.0000000 320.0000000 330.0000000 340.0000000 350.0000000 360.0000000 370.0000000 380.0000000 390.0000000 400.0000000 410.0000000 420.0000000 430.0000000 440.0000000 450.0000000 460.0000000 470.0000000 480.0000000 490.0000000 500.0000000 510.0000000 520.0000000 530.0000000 540.0000000 550.0000000 560.0000000 570.0000000 580.0000000 590.0000000 600.0000000 610.0000000 620.0000000 630.0000000 640.0000000 650.0000000 660.0000000 670.0000000 680.0000000 690.0000000 700.0000000 710.0000000 720.0000000 730.0000000 740.0000000 750.0000000 760.0000000 770.0000000 780.0000000 790.0000000 800.0000000 810.0000000 
+DEAL::[block 1 ]  411.0000000 417.0000000 423.0000000 429.0000000 435.0000000 441.0000000 447.0000000 453.0000000 459.0000000 465.0000000 471.0000000 477.0000000 483.0000000 489.0000000 495.0000000 501.0000000 507.0000000 513.0000000 519.0000000 525.0000000 531.0000000 537.0000000 543.0000000 549.0000000 555.0000000 
+DEAL::Block vector: B.vmult_add:
+DEAL::[block 0 ]  256.0000000 269.0000000 282.0000000 295.0000000 308.0000000 321.0000000 334.0000000 347.0000000 360.0000000 373.0000000 386.0000000 399.0000000 412.0000000 425.0000000 438.0000000 451.0000000 464.0000000 477.0000000 490.0000000 503.0000000 516.0000000 529.0000000 542.0000000 555.0000000 568.0000000 260.0000000 270.0000000 280.0000000 290.0000000 300.0000000 310.0000000 320.0000000 330.0000000 340.0000000 350.0000000 360.0000000 370.0000000 380.0000000 390.0000000 400.0000000 410.0000000 420.0000000 430.0000000 440.0000000 450.0000000 460.0000000 470.0000000 480.0000000 490.0000000 500.0000000 510.0000000 520.0000000 530.0000000 540.0000000 550.0000000 560.0000000 570.0000000 580.0000000 590.0000000 600.0000000 610.0000000 620.0000000 630.0000000 640.0000000 650.0000000 660.0000000 670.0000000 680.0000000 690.0000000 700.0000000 710.0000000 720.0000000 730.0000000 740.0000000 750.0000000 760.0000000 770.0000000 780.0000000 790.0000000 800.0000000 810.0000000 
+DEAL::[block 1 ]  411.0000000 417.0000000 423.0000000 429.0000000 435.0000000 441.0000000 447.0000000 453.0000000 459.0000000 465.0000000 471.0000000 477.0000000 483.0000000 489.0000000 495.0000000 501.0000000 507.0000000 513.0000000 519.0000000 525.0000000 531.0000000 537.0000000 543.0000000 549.0000000 555.0000000 
+DEAL::Block vector: A.Tvmult_add:
+DEAL::[block 0 ]  93.00000000 105.0000000 117.0000000 129.0000000 141.0000000 153.0000000 165.0000000 177.0000000 189.0000000 201.0000000 213.0000000 225.0000000 237.0000000 249.0000000 261.0000000 273.0000000 285.0000000 297.0000000 309.0000000 321.0000000 333.0000000 345.0000000 357.0000000 369.0000000 381.0000000 286.0000000 297.0000000 308.0000000 319.0000000 330.0000000 341.0000000 352.0000000 363.0000000 374.0000000 385.0000000 396.0000000 407.0000000 418.0000000 429.0000000 440.0000000 451.0000000 462.0000000 473.0000000 484.0000000 495.0000000 506.0000000 517.0000000 528.0000000 539.0000000 550.0000000 561.0000000 572.0000000 583.0000000 594.0000000 605.0000000 616.0000000 627.0000000 638.0000000 649.0000000 660.0000000 671.0000000 682.0000000 693.0000000 704.0000000 715.0000000 726.0000000 737.0000000 748.0000000 759.0000000 770.0000000 781.0000000 792.0000000 803.0000000 814.0000000 825.0000000 836.0000000 847.0000000 858.0000000 869.0000000 880.0000000 891.0000000 
+DEAL::[block 1 ]  495.0000000 504.0000000 513.0000000 522.0000000 531.0000000 540.0000000 549.0000000 558.0000000 567.0000000 576.0000000 585.0000000 594.0000000 603.0000000 612.0000000 621.0000000 630.0000000 639.0000000 648.0000000 657.0000000 666.0000000 675.0000000 684.0000000 693.0000000 702.0000000 711.0000000 
+DEAL::Block vector: B.Tvmult_add:
+DEAL::[block 0 ]  93.00000000 105.0000000 117.0000000 129.0000000 141.0000000 153.0000000 165.0000000 177.0000000 189.0000000 201.0000000 213.0000000 225.0000000 237.0000000 249.0000000 261.0000000 273.0000000 285.0000000 297.0000000 309.0000000 321.0000000 333.0000000 345.0000000 357.0000000 369.0000000 381.0000000 286.0000000 297.0000000 308.0000000 319.0000000 330.0000000 341.0000000 352.0000000 363.0000000 374.0000000 385.0000000 396.0000000 407.0000000 418.0000000 429.0000000 440.0000000 451.0000000 462.0000000 473.0000000 484.0000000 495.0000000 506.0000000 517.0000000 528.0000000 539.0000000 550.0000000 561.0000000 572.0000000 583.0000000 594.0000000 605.0000000 616.0000000 627.0000000 638.0000000 649.0000000 660.0000000 671.0000000 682.0000000 693.0000000 704.0000000 715.0000000 726.0000000 737.0000000 748.0000000 759.0000000 770.0000000 781.0000000 792.0000000 803.0000000 814.0000000 825.0000000 836.0000000 847.0000000 858.0000000 869.0000000 880.0000000 891.0000000 
+DEAL::[block 1 ]  495.0000000 504.0000000 513.0000000 522.0000000 531.0000000 540.0000000 549.0000000 558.0000000 567.0000000 576.0000000 585.0000000 594.0000000 603.0000000 612.0000000 621.0000000 630.0000000 639.0000000 648.0000000 657.0000000 666.0000000 675.0000000 684.0000000 693.0000000 702.0000000 711.0000000