Fix Vector::operator=(TrilinosWrappers::MPI::Vector)

Data was being copied to the local vector directly from the distributed
vector (as opposed to its localised equivalent). This lead to only a
partial copy of data, and all of the data was offset incorrectly in the
returned vector (giving each process a different view of the localised
vector).

diff --git a/doc/news/changes/minor/20180307Jean-PaulPelteret b/doc/news/changes/minor/20180307Jean-PaulPelteret
new file mode 100644 (file)
index 0000000..804b36b
--- /dev/null
+++ b/doc/news/changes/minor/20180307Jean-PaulPelteret
@@ -0,0 +1,7 @@
+Fixed: Copying a Trilinos::MPI::Vector to a local deal.II Vector using 
+operator=() resulted in only a partial copy of data to the local vector. In 
+addition, the copied  elements  were offset incorrectly on each process. 
+Similar held for copying a  Trilinos::MPI::BlockVector to a local deal.II 
+BlockVector. This has been fixed and now works as expected.
+<br>
+(Jean-Paul Pelteret, 2018/03/07)

diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h
index 393817135f0ed88e0491984cc1ec94351d41eba4..ebb917a35e152f5bf6d360d988926febd399e352 100644 (file)
--- a/include/deal.II/lac/vector.templates.h
+++ b/include/deal.II/lac/vector.templates.h
@@ -180,6 +180,9 @@ Vector<Number>::Vector (const TrilinosWrappers::MPI::Vector &v)
       localized_vector.reinit(complete_index_set(vec_size), v.get_mpi_communicator());
       localized_vector.reinit (v, false, true);
 
+      Assert(localized_vector.size() == vec_size,
+             ExcDimensionMismatch(localized_vector.size(), vec_size));
+
       // get a representation of the vector
       // and copy it
       TrilinosScalar **start_ptr;
@@ -825,22 +828,29 @@ template <typename Number>
 Vector<Number> &
 Vector<Number>::operator= (const TrilinosWrappers::MPI::Vector &v)
 {
-  // Generate a localized version
-  // of the Trilinos vectors and
-  // then call the other =
-  // operator.
-  TrilinosWrappers::MPI::Vector localized_vector;
-  localized_vector.reinit(complete_index_set(v.size()), v.get_mpi_communicator());
-  localized_vector.reinit(v, false, true);
-
   if (v.size() != vec_size)
     reinit (v.size(), true);
   if (vec_size != 0)
     {
+      // Copy the distributed vector to
+      // a local one at all processors
+      // that know about the original vector.
+      // TODO: There could
+      // be a better solution than
+      // this, but it has not yet been
+      // found.
+      TrilinosWrappers::MPI::Vector localized_vector;
+      localized_vector.reinit(complete_index_set(vec_size), v.get_mpi_communicator());
+      localized_vector.reinit (v, false, true);
+
+      Assert(localized_vector.size() == vec_size,
+             ExcDimensionMismatch(localized_vector.size(), vec_size));
+
       // get a representation of the vector
       // and copy it
       TrilinosScalar **start_ptr;
-      int ierr = v.trilinos_vector().ExtractView (&start_ptr);
+
+      int ierr = localized_vector.trilinos_vector().ExtractView (&start_ptr);
       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
       std::copy (start_ptr[0], start_ptr[0]+vec_size, begin());

diff --git a/tests/trilinos/parallel_block_vector_copy_01.cc b/tests/trilinos/parallel_block_vector_copy_01.cc
new file mode 100644 (file)
index 0000000..f61401c
--- /dev/null
+++ b/tests/trilinos/parallel_block_vector_copy_01.cc
@@ -0,0 +1,75 @@
+// ---------------------------------------------------------------------
+//
+// 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// make sure that block vector copies correctly to a serial vector
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <iostream>
+
+
+int main (int argc,char **argv)
+{
+  typedef BlockVector<double> BlockVectorLocal;
+  typedef TrilinosWrappers::MPI::BlockVector BlockVector;
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+  MPILogInitAll log;
+
+  const MPI_Comm &mpi_communicator = MPI_COMM_WORLD;
+  const unsigned int this_mpi_process = Utilities::MPI::this_mpi_process(mpi_communicator);
+  const unsigned int n_mpi_processes = Utilities::MPI::n_mpi_processes(mpi_communicator);
+
+  const unsigned int n_blocks = 2;
+  const unsigned int n_dofs_per_block = 10;
+
+  std::vector<IndexSet> locally_owned_partitioning (n_blocks);
+  for (unsigned int b=0; b<n_blocks; ++b)
+    {
+      locally_owned_partitioning[b].set_size(n_dofs_per_block);
+      locally_owned_partitioning[b].add_range(this_mpi_process*(n_dofs_per_block/2),(this_mpi_process+1)*(n_dofs_per_block/2));
+    }
+
+  BlockVector parallel_vector;
+  parallel_vector.reinit(locally_owned_partitioning);
+  deallog << "Locally owned indices" << std::endl;
+  parallel_vector.locally_owned_elements().print(deallog.get_file_stream());
+
+  // Set entries in parallel vector
+  for (auto idx : parallel_vector.locally_owned_elements())
+    {
+      parallel_vector[idx] = 10.0*idx;
+    }
+  deallog << "Parallel vector" << std::endl;
+  parallel_vector.print(deallog.get_file_stream());
+
+  // Copy distributed vector to local vector
+  deallog << "Localized vector (copy constructor)" << std::endl;
+  const BlockVectorLocal local_vector_1 (parallel_vector);
+  local_vector_1.print(deallog.get_file_stream());
+
+  // Copy distributed vector to local vector
+  deallog << "Localized vector (operator =)" << std::endl;
+  const BlockVectorLocal local_vector_2 = parallel_vector;
+  local_vector_2.print(deallog.get_file_stream());
+
+  Assert(local_vector_1 == local_vector_2, ExcMessage("Vectors don't match"));
+
+  deallog << "OK" << std::endl;
+}

diff --git a/tests/trilinos/parallel_block_vector_copy_01.mpirun=2.output b/tests/trilinos/parallel_block_vector_copy_01.mpirun=2.output
new file mode 100644 (file)
index 0000000..9a2236b
--- /dev/null
+++ b/tests/trilinos/parallel_block_vector_copy_01.mpirun=2.output
@@ -0,0 +1,47 @@
+
+DEAL:0::Locally owned indices
+{[0,4], [10,14]}
+DEAL:0::Parallel vector
+C0:size:10 local_size:5 :
+[0]: 0.000e+00
+[1]: 1.000e+01
+[2]: 2.000e+01
+[3]: 3.000e+01
+[4]: 4.000e+01
+C1:size:10 local_size:5 :
+[0]: 1.000e+02
+[1]: 1.100e+02
+[2]: 1.200e+02
+[3]: 1.300e+02
+[4]: 1.400e+02
+DEAL:0::Localized vector (copy constructor)
+C0:0.000e+00 1.000e+01 2.000e+01 3.000e+01 4.000e+01 5.000e+01 6.000e+01 7.000e+01 8.000e+01 9.000e+01 
+C1:1.000e+02 1.100e+02 1.200e+02 1.300e+02 1.400e+02 1.500e+02 1.600e+02 1.700e+02 1.800e+02 1.900e+02 
+DEAL:0::Localized vector (operator =)
+C0:0.000e+00 1.000e+01 2.000e+01 3.000e+01 4.000e+01 5.000e+01 6.000e+01 7.000e+01 8.000e+01 9.000e+01 
+C1:1.000e+02 1.100e+02 1.200e+02 1.300e+02 1.400e+02 1.500e+02 1.600e+02 1.700e+02 1.800e+02 1.900e+02 
+DEAL:0::OK
+
+DEAL:1::Locally owned indices
+{[5,9], [15,19]}
+DEAL:1::Parallel vector
+C0:size:10 local_size:5 :
+[5]: 5.000e+01
+[6]: 6.000e+01
+[7]: 7.000e+01
+[8]: 8.000e+01
+[9]: 9.000e+01
+C1:size:10 local_size:5 :
+[5]: 1.500e+02
+[6]: 1.600e+02
+[7]: 1.700e+02
+[8]: 1.800e+02
+[9]: 1.900e+02
+DEAL:1::Localized vector (copy constructor)
+C0:0.000e+00 1.000e+01 2.000e+01 3.000e+01 4.000e+01 5.000e+01 6.000e+01 7.000e+01 8.000e+01 9.000e+01 
+C1:1.000e+02 1.100e+02 1.200e+02 1.300e+02 1.400e+02 1.500e+02 1.600e+02 1.700e+02 1.800e+02 1.900e+02 
+DEAL:1::Localized vector (operator =)
+C0:0.000e+00 1.000e+01 2.000e+01 3.000e+01 4.000e+01 5.000e+01 6.000e+01 7.000e+01 8.000e+01 9.000e+01 
+C1:1.000e+02 1.100e+02 1.200e+02 1.300e+02 1.400e+02 1.500e+02 1.600e+02 1.700e+02 1.800e+02 1.900e+02 
+DEAL:1::OK
+