From: Bruno Turcksin Date: Wed, 1 Aug 2018 18:09:51 +0000 (+0000) Subject: Fix a bug in distributed::Vector::import X-Git-Tag: v9.1.0-rc1~854^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2b7336b93edd0da72c5d6b8abdac5df1dbf2e82a;p=dealii.git Fix a bug in distributed::Vector::import The import function would fail if the ReadWriteVector would only contains ghosted values --- diff --git a/include/deal.II/lac/la_parallel_vector.templates.h b/include/deal.II/lac/la_parallel_vector.templates.h index d4d233d820..2898463552 100644 --- a/include/deal.II/lac/la_parallel_vector.templates.h +++ b/include/deal.II/lac/la_parallel_vector.templates.h @@ -36,6 +36,49 @@ namespace LinearAlgebra { namespace distributed { + namespace internal + { + // In the import_from_ghosted_array_finish we might need to calculate the + // maximal and minimal value for the given number type, which is not + // straight forward for complex numbers. Therefore, comparison of complex + // numbers is prohibited and throws an exception. + template + Number + get_min(const Number a, const Number b) + { + return std::min(a, b); + } + + template + std::complex + get_min(const std::complex a, const std::complex) + { + AssertThrow(false, + ExcMessage("VectorOperation::min not " + "implemented for complex numbers")); + return a; + } + + template + Number + get_max(const Number a, const Number b) + { + return std::max(a, b); + } + + template + std::complex + get_max(const std::complex a, const std::complex) + { + AssertThrow(false, + ExcMessage("VectorOperation::max not " + "implemented for complex numbers")); + return a; + } + } // namespace internal + + + template void Vector::clear_mpi_requests() @@ -669,8 +712,7 @@ namespace LinearAlgebra // then create the communication pattern IndexSet ghost_indices(V.get_stored_elements()); ghost_indices.subtract_set(locally_owned_elem); - IndexSet local_indices(V.get_stored_elements()); - local_indices.subtract_set(ghost_indices); + IndexSet local_indices(locally_owned_elem); comm_pattern = std::make_shared( local_indices, ghost_indices, get_mpi_communicator()); } @@ -683,28 +725,61 @@ namespace LinearAlgebra ExcMessage("The communication pattern is not of type " "Utilities::MPI::Partitioner.")); } + IndexSet ghost_indices(V.get_stored_elements()); + ghost_indices.subtract_set(locally_owned_elem); + IndexSet local_indices(locally_owned_elem); Vector tmp_vector(comm_pattern); + std::copy(begin(), end(), tmp_vector.begin()); // fill entries from ReadWriteVector into the distributed vector, // including ghost entries. this is not really efficient right now // because indices are translated twice, once by nth_index_in_set(i) and // once for operator() of tmp_vector - const IndexSet &v_stored = V.get_stored_elements(); - for (size_type i = 0; i < v_stored.n_elements(); ++i) - tmp_vector(v_stored.nth_index_in_set(i)) = V.local_element(i); + const IndexSet &v_stored = V.get_stored_elements(); + const size_type v_n_elements = v_stored.n_elements(); + switch (operation) + { + case VectorOperation::insert: + { + for (size_type i = 0; i < v_n_elements; ++i) + tmp_vector(v_stored.nth_index_in_set(i)) = V.local_element(i); - tmp_vector.compress(operation); + break; + } + case VectorOperation::add: + { + for (size_type i = 0; i < v_n_elements; ++i) + tmp_vector(v_stored.nth_index_in_set(i)) += V.local_element(i); - // Copy the local elements of tmp_vector to the right place in val - IndexSet tmp_index_set = tmp_vector.locally_owned_elements(); - for (size_type i = 0; i < tmp_index_set.n_elements(); ++i) - { - values[locally_owned_elem.index_within_set( - tmp_index_set.nth_index_in_set(i))] = tmp_vector.local_element(i); - } - } + break; + } + case VectorOperation::min: + { + for (size_type i = 0; i < v_n_elements; ++i) + tmp_vector(v_stored.nth_index_in_set(i)) = + internal::get_min(tmp_vector(v_stored.nth_index_in_set(i)), + V.local_element(i)); + + break; + } + case VectorOperation::max: + { + for (size_type i = 0; i < v_n_elements; ++i) + tmp_vector(v_stored.nth_index_in_set(i)) = + internal::get_max(tmp_vector(v_stored.nth_index_in_set(i)), + V.local_element(i)); + break; + } + default: + { + Assert(false, ExcMessage("This operation is not supported.")); + } + } + tmp_vector.compress(operation); + std::copy(tmp_vector.begin(), tmp_vector.end(), begin()); + } template void diff --git a/tests/mpi/parallel_vector_24.cc b/tests/mpi/parallel_vector_24.cc new file mode 100644 index 0000000000..9d510bcd28 --- /dev/null +++ b/tests/mpi/parallel_vector_24.cc @@ -0,0 +1,143 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +// test based on parallel_vector_03a but using import function. This test was +// not working because we would only set ghost elements and no local elements +// which resulted in an error during the import + +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +void +check(const unsigned int myid, + const LinearAlgebra::distributed::Vector &v) +{ + v.print(std::cout); + if (myid == 0) + { + AssertThrow(v(10) == 10.0, ExcInternalError()); + AssertThrow(v(11) == 0., ExcInternalError()); + AssertThrow(v(12) == 0., ExcInternalError()); + AssertThrow(v(14) == 14., ExcInternalError()); + + AssertThrow(v(5) == 55., ExcInternalError()); + } + else + { + AssertThrow(v(4) == 0., ExcInternalError()); + AssertThrow(v(5) == 55., ExcInternalError()); + AssertThrow(v(6) == 66., ExcInternalError()); + } +} + + +void +test() +{ + unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + + Assert(numproc == 2, ExcNotImplemented()); + + const unsigned int size = 20; + IndexSet local_owned(size); + IndexSet local_nonzero(size); + IndexSet local_relevant(size); + if (myid == 0) + { + local_owned.add_range(0, 10); + local_nonzero.add_range(5, 10); + local_relevant = local_owned; + local_relevant.add_range(10, 13); + local_relevant.add_range(14, 15); + } + else + { + local_owned.add_range(10, size); + local_nonzero.add_range(10, 11); + local_nonzero.add_range(13, 15); + local_relevant = local_owned; + local_relevant.add_range(4, 7); + } + + LinearAlgebra::distributed::Vector v(local_owned, + local_relevant, + MPI_COMM_WORLD); + v = 0.; + + // set local values + for (unsigned int i = 0; i < local_nonzero.n_elements(); i++) + v(local_nonzero.nth_index_in_set(i)) = local_nonzero.nth_index_in_set(i); + + // set value from processor which does not own it: + v(5) = 55.; + v.compress(VectorOperation::insert); + + // add to value from processor which has it as a ghost + IndexSet ghost_set(size); + LinearAlgebra::ReadWriteVector rw_vector; + if (myid == 1) + { + ghost_set.add_index(6); + rw_vector.reinit(ghost_set); + rw_vector(6) = 60; + } + else + { + rw_vector.reinit(ghost_set); + } + v.import(rw_vector, VectorOperation::add); // 60 + 6 + // compress(insert) used to leave ghosts un-touched which resulted in + // the wrong 55+55 for this compress(add) operation. + + v.update_ghost_values(); + + check(myid, v); + + if (myid == 0) + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + deallog.push(Utilities::int_to_string(myid)); + + if (myid == 0) + { + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog << std::setprecision(4); + + test(); + } + else + test(); +} diff --git a/tests/mpi/parallel_vector_24.mpirun=2.output b/tests/mpi/parallel_vector_24.mpirun=2.output new file mode 100644 index 0000000000..be8d055f86 --- /dev/null +++ b/tests/mpi/parallel_vector_24.mpirun=2.output @@ -0,0 +1,2 @@ + +DEAL:0::OK