From: Svenja Schoeder Date: Mon, 18 Jun 2018 14:24:22 +0000 (+0200) Subject: Introduce min and max to VectorOperation for parallel vectors X-Git-Tag: v9.1.0-rc1~945^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=23ffb7b74f7b32741fbccd347da461ce8c0c4842;p=dealii.git Introduce min and max to VectorOperation for parallel vectors Add tests for min max VectorOperation add changelog disallow complex numbers for VectorOperation::min and VectorOperation::max testing of VectorOperation::max if several processors have same ghost element but with different values fixed typo on cases --- diff --git a/doc/news/changes/minor/20180618SvenjaSchoeder b/doc/news/changes/minor/20180618SvenjaSchoeder new file mode 100644 index 0000000000..7c5f98b574 --- /dev/null +++ b/doc/news/changes/minor/20180618SvenjaSchoeder @@ -0,0 +1,6 @@ +New: The LinearAlgebra::distributed::Vector::compress() function has +gained two new combine operations VectorOperation::min and +VectorOperation::max to define the entry at the owning processor as +the minimum/maximum of its own value and the one coming from a ghost. +
+(Svenja Schoeder, 2018/06/18) diff --git a/include/deal.II/base/partitioner.templates.h b/include/deal.II/base/partitioner.templates.h index 2b5c8bfb16..87b5c601ff 100644 --- a/include/deal.II/base/partitioner.templates.h +++ b/include/deal.II/base/partitioner.templates.h @@ -329,6 +329,43 @@ namespace Utilities return a; } + // In the import_from_ghosted_array_finish we need to calculate maximal + // and minimal value on number types, which is not straight forward for + // complex numbers. Therfore, comparison of complex numbers is + // prohibited and throw an assert. + 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 max 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::min max not " + "implemented for complex numbers")); + return a; + } } // namespace internal @@ -403,11 +440,29 @@ namespace Utilities // local values. For insert, nothing is done here (but in debug mode // we assert that the specified value is either zero or matches with // the ones already present - if (vector_operation != dealii::VectorOperation::insert) + if (vector_operation == dealii::VectorOperation::add) for (; my_imports != import_indices_data.end(); ++my_imports) for (unsigned int j = my_imports->first; j < my_imports->second; j++) locally_owned_array[j] += *read_position++; + else if (vector_operation == dealii::VectorOperation::min) + for (; my_imports != import_indices_data.end(); ++my_imports) + for (unsigned int j = my_imports->first; j < my_imports->second; + j++) + { + locally_owned_array[j] = + internal::get_min(*read_position, locally_owned_array[j]); + read_position++; + } + else if (vector_operation == dealii::VectorOperation::max) + for (; my_imports != import_indices_data.end(); ++my_imports) + for (unsigned int j = my_imports->first; j < my_imports->second; + j++) + { + locally_owned_array[j] = + internal::get_max(*read_position, locally_owned_array[j]); + read_position++; + } else for (; my_imports != import_indices_data.end(); ++my_imports) for (unsigned int j = my_imports->first; j < my_imports->second; diff --git a/include/deal.II/lac/la_parallel_vector.h b/include/deal.II/lac/la_parallel_vector.h index 819d88cb5e..0a7fbc9dfd 100644 --- a/include/deal.II/lac/la_parallel_vector.h +++ b/include/deal.II/lac/la_parallel_vector.h @@ -402,7 +402,7 @@ namespace LinearAlgebra * @ref GlossCompress "Compressing distributed vectors and matrices" * in the glossary. * - * There are two variants for this function. If called with argument @p + * There are four variants for this function. If called with argument @p * VectorOperation::add adds all the data accumulated in ghost elements * to the respective elements on the owning processor and clears the * ghost array afterwards. If called with argument @p @@ -416,6 +416,13 @@ namespace LinearAlgebra * actually consistent between processors, i.e., whenever a non-zero * ghost element is found, it is compared to the value on the owning * processor and an exception is thrown if these elements do not agree. + * If called with VectorOperation::min or VectorOperation::max the + * minimum or maximum on all elements across the processors is set. If + * a processor does not set values for the ghosts, the value 0.0 is + * assumed, which can be the minimal or maximal value across all + * processors. The operations min and max are reasonable if all + * processors set relevant values or set values that do not alter the + * results, e.g. small values for the max operation. */ virtual void compress(::dealii::VectorOperation::values operation) override; diff --git a/include/deal.II/lac/read_write_vector.templates.h b/include/deal.II/lac/read_write_vector.templates.h index d577e4c7c7..62608caf10 100644 --- a/include/deal.II/lac/read_write_vector.templates.h +++ b/include/deal.II/lac/read_write_vector.templates.h @@ -250,7 +250,46 @@ namespace LinearAlgebra return *this; } + namespace internal + { + // In the import_from_ghosted_array_finish we need to calculate maximal + // and minimal value on number types, which is not straight forward for + // complex numbers. Therfore, comparison of complex numbers is + // prohibited and throw an assert. + 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 max 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::min max not " + "implemented for complex numbers")); + return a; + } + } // namespace internal template void @@ -288,6 +327,16 @@ namespace LinearAlgebra if (operation == VectorOperation::add) for (size_type i = 0; i < stored.n_elements(); ++i) local_element(i) += tmp_vector(stored.nth_index_in_set(i)); + else if (operation == VectorOperation::min) + for (size_type i = 0; i < stored.n_elements(); ++i) + local_element(i) = + internal::get_min(tmp_vector(stored.nth_index_in_set(i)), + local_element(i)); + else if (operation == VectorOperation::max) + for (size_type i = 0; i < stored.n_elements(); ++i) + local_element(i) = + internal::get_max(tmp_vector(stored.nth_index_in_set(i)), + local_element(i)); else for (size_type i = 0; i < stored.n_elements(); ++i) local_element(i) = tmp_vector(stored.nth_index_in_set(i)); @@ -438,6 +487,42 @@ namespace LinearAlgebra for (int i = 0; i < size; ++i) values[i] += new_values[i]; } + else if (operation == VectorOperation::min) + { + const int err = target_vector.Import(multivector, import, Add); + AssertThrow(err == 0, + ExcMessage("Epetra Import() failed with error code: " + + Utilities::to_string(err))); + + const double *new_values = target_vector.Values(); + const int size = target_vector.MyLength(); + Assert(size == 0 || values != nullptr, + ExcInternalError("Import failed.")); + + // To ensure that this code also compiles with complex + // numbers, we only compare the real part of the + // variable. Note that min/max do not make sense with complex + // numbers. + for (int i = 0; i < size; ++i) + if (std::real(new_values[i]) - std::real(values[i]) < 0.0) + values[i] = new_values[i]; + } + else if (operation == VectorOperation::max) + { + const int err = target_vector.Import(multivector, import, Add); + AssertThrow(err == 0, + ExcMessage("Epetra Import() failed with error code: " + + Utilities::to_string(err))); + + const double *new_values = target_vector.Values(); + const int size = target_vector.MyLength(); + Assert(size == 0 || values != nullptr, + ExcInternalError("Import failed.")); + + for (int i = 0; i < size; ++i) + if (std::real(new_values[i]) - std::real(values[i]) > 0.0) + values[i] = new_values[i]; + } else AssertThrow(false, ExcNotImplemented()); } @@ -500,7 +585,7 @@ namespace LinearAlgebra cudaMemcpyDeviceToHost); AssertCuda(error_code); } - else + else if (operation == VectorOperation::add) { // Copy the vector from the device to a temporary vector on the host std::vector tmp(n_elements); @@ -514,6 +599,40 @@ namespace LinearAlgebra for (unsigned int i = 0; i < n_elements; ++i) values[i] += tmp[i]; } + else if (operation == VectorOperation::min) + { + // Copy the vector from the device to a temporary vector on the host + std::vector tmp(n_elements); + cudaError_t error_code = cudaMemcpy(tmp.data(), + cuda_vec.get_values(), + n_elements * sizeof(Number), + cudaMemcpyDeviceToHost); + AssertCuda(error_code); + + // To ensure that this code also compiles with complex + // numbers, we only compare the real part of the + // variable. Note that min/max do not make sense with complex + // numbers. + for (unsigned int i = 0; i < n_elements; ++i) + if (std::real(tmp[i]) - std::real(values[i]) < 0.0) + values[i] = tmp[i]; + } + else if (operation == VectorOperation::max) + { + // Copy the vector from the device to a temporary vector on the host + std::vector tmp(n_elements); + cudaError_t error_code = cudaMemcpy(tmp.data(), + cuda_vec.get_values(), + n_elements * sizeof(Number), + cudaMemcpyDeviceToHost); + AssertCuda(error_code); + + for (unsigned int i = 0; i < n_elements; ++i) + if (std::real(tmp[i]) - std::real(values[i]) > 0.0) + values[i] = tmp[i]; + } + else + AssertThrow(false, ExcNotImplemented()); } #endif diff --git a/include/deal.II/lac/vector_operation.h b/include/deal.II/lac/vector_operation.h index 6e31772cfb..4b3fc163de 100644 --- a/include/deal.II/lac/vector_operation.h +++ b/include/deal.II/lac/vector_operation.h @@ -50,7 +50,15 @@ struct VectorOperation /** * The current operation is an addition. */ - add + add, + /** + * The current operation is a minimization. + */ + min, + /** + * The current operation is a maximization. + */ + max }; }; diff --git a/source/lac/cuda_vector.cu b/source/lac/cuda_vector.cu index 7f6b1c40e2..725dc483c3 100644 --- a/source/lac/cuda_vector.cu +++ b/source/lac/cuda_vector.cu @@ -615,7 +615,7 @@ namespace LinearAlgebra cudaMemcpyHostToDevice); AssertCuda(error_code); } - else + else if (operation == VectorOperation::add) { // Create a temporary vector on the device Number * tmp; @@ -644,6 +644,8 @@ namespace LinearAlgebra error_code = cudaFree(tmp); AssertCuda(error_code); } + else + AssertThrow(false, ExcNotImplemented()); } diff --git a/source/lac/trilinos_epetra_vector.cc b/source/lac/trilinos_epetra_vector.cc index 654a3668ed..36d38cf164 100644 --- a/source/lac/trilinos_epetra_vector.cc +++ b/source/lac/trilinos_epetra_vector.cc @@ -187,8 +187,14 @@ namespace LinearAlgebra if (operation == VectorOperation::insert) vector->Export(source_vector, import, Insert); - else + else if (operation == VectorOperation::add) vector->Export(source_vector, import, Add); + else if (operation == VectorOperation::max) + vector->Export(source_vector, import, Epetra_Max); + else if (operation == VectorOperation::min) + vector->Export(source_vector, import, Epetra_Min); + else + AssertThrow(false, ExcNotImplemented()); } diff --git a/source/lac/trilinos_sparse_matrix.cc b/source/lac/trilinos_sparse_matrix.cc index fe657772ef..1eee3a76cd 100644 --- a/source/lac/trilinos_sparse_matrix.cc +++ b/source/lac/trilinos_sparse_matrix.cc @@ -1155,6 +1155,10 @@ namespace TrilinosWrappers mode = Add; else if (operation == ::dealii::VectorOperation::insert) mode = Insert; + else + Assert( + false, + "compress() can only be called with VectorOperation add, insert, or unknown"); } else { diff --git a/source/lac/trilinos_vector.cc b/source/lac/trilinos_vector.cc index f2acfcec4b..f143f26209 100644 --- a/source/lac/trilinos_vector.cc +++ b/source/lac/trilinos_vector.cc @@ -560,6 +560,10 @@ namespace TrilinosWrappers mode = Add; else if (given_last_action == ::dealii::VectorOperation::insert) mode = Insert; + else + Assert( + false, + "compress() can only be called with VectorOperation add, insert, or unknown"); } else { diff --git a/tests/mpi/parallel_vector_22.cc b/tests/mpi/parallel_vector_22.cc new file mode 100644 index 0000000000..ecae1101db --- /dev/null +++ b/tests/mpi/parallel_vector_22.cc @@ -0,0 +1,114 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +// check LA::Vector::compress(VectorOperation::min/max) from ghosts + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + + +void +test() +{ + unsigned int myid = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + unsigned int numproc = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + + if (myid == 0) + deallog << "numproc=" << numproc << std::endl; + + + // each processor owns 2 indices and all + // are ghosting element 1 (the second) + IndexSet local_owned(numproc * 2); + local_owned.add_range(myid * 2, myid * 2 + 2); + IndexSet local_relevant(numproc * 2); + local_relevant = local_owned; + local_relevant.add_range(1, 2); + + // create vector + LinearAlgebra::distributed::Vector v(local_owned, + local_relevant, + MPI_COMM_WORLD); + + // set local values + v(myid * 2) = myid * 2.0; + v(myid * 2 + 1) = myid * 2.0 + 1.0; + v.compress(VectorOperation::add); + v *= 2.0; + + // check setup of vectors + deallog << myid << ":" + << "first owned entry: " << v(myid * 2) << std::endl; + deallog << myid << ":" + << "second owned entry: " << v(myid * 2 + 1) << std::endl; + + // set ghost dof on not owning processor and maximize + if (myid) + v(1) = 7. * myid; + v.compress(VectorOperation::max); + + // import ghosts onto all procs + v.update_ghost_values(); + + // check + deallog << myid << ":" + << "ghost entry: " << v(1) << std::endl; + + // ghosts are set to zero + v.zero_out_ghosts(); + + // minimize + v.compress(VectorOperation::min); + v.update_ghost_values(); + + // check + deallog << myid << ":" + << "ghost entry: " << v(1) << std::endl; + + // update of ghost value from owner and minimize + v.zero_out_ghosts(); + if (!myid) + v(1) = -1.; + v.compress(VectorOperation::min); + v.update_ghost_values(); + + // check + deallog << myid << ":" + << "ghost entry: " << v(1) << std::endl; + + 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()); + + MPILogInitAll log; + + test(); +} diff --git a/tests/mpi/parallel_vector_22.mpirun=10.output b/tests/mpi/parallel_vector_22.mpirun=10.output new file mode 100644 index 0000000000..8a13c7f82a --- /dev/null +++ b/tests/mpi/parallel_vector_22.mpirun=10.output @@ -0,0 +1,71 @@ + +DEAL:0::numproc=10 +DEAL:0::0:first owned entry: 0.00000 +DEAL:0::0:second owned entry: 2.00000 +DEAL:0::0:ghost entry: 63.0000 +DEAL:0::0:ghost entry: 0.00000 +DEAL:0::0:ghost entry: -1.00000 +DEAL:0::OK + +DEAL:1::1:first owned entry: 4.00000 +DEAL:1::1:second owned entry: 6.00000 +DEAL:1::1:ghost entry: 63.0000 +DEAL:1::1:ghost entry: 0.00000 +DEAL:1::1:ghost entry: -1.00000 + + +DEAL:2::2:first owned entry: 8.00000 +DEAL:2::2:second owned entry: 10.0000 +DEAL:2::2:ghost entry: 63.0000 +DEAL:2::2:ghost entry: 0.00000 +DEAL:2::2:ghost entry: -1.00000 + + +DEAL:3::3:first owned entry: 12.0000 +DEAL:3::3:second owned entry: 14.0000 +DEAL:3::3:ghost entry: 63.0000 +DEAL:3::3:ghost entry: 0.00000 +DEAL:3::3:ghost entry: -1.00000 + + +DEAL:4::4:first owned entry: 16.0000 +DEAL:4::4:second owned entry: 18.0000 +DEAL:4::4:ghost entry: 63.0000 +DEAL:4::4:ghost entry: 0.00000 +DEAL:4::4:ghost entry: -1.00000 + + +DEAL:5::5:first owned entry: 20.0000 +DEAL:5::5:second owned entry: 22.0000 +DEAL:5::5:ghost entry: 63.0000 +DEAL:5::5:ghost entry: 0.00000 +DEAL:5::5:ghost entry: -1.00000 + + +DEAL:6::6:first owned entry: 24.0000 +DEAL:6::6:second owned entry: 26.0000 +DEAL:6::6:ghost entry: 63.0000 +DEAL:6::6:ghost entry: 0.00000 +DEAL:6::6:ghost entry: -1.00000 + + +DEAL:7::7:first owned entry: 28.0000 +DEAL:7::7:second owned entry: 30.0000 +DEAL:7::7:ghost entry: 63.0000 +DEAL:7::7:ghost entry: 0.00000 +DEAL:7::7:ghost entry: -1.00000 + + +DEAL:8::8:first owned entry: 32.0000 +DEAL:8::8:second owned entry: 34.0000 +DEAL:8::8:ghost entry: 63.0000 +DEAL:8::8:ghost entry: 0.00000 +DEAL:8::8:ghost entry: -1.00000 + + +DEAL:9::9:first owned entry: 36.0000 +DEAL:9::9:second owned entry: 38.0000 +DEAL:9::9:ghost entry: 63.0000 +DEAL:9::9:ghost entry: 0.00000 +DEAL:9::9:ghost entry: -1.00000 + diff --git a/tests/mpi/parallel_vector_22.mpirun=4.output b/tests/mpi/parallel_vector_22.mpirun=4.output new file mode 100644 index 0000000000..4c36987598 --- /dev/null +++ b/tests/mpi/parallel_vector_22.mpirun=4.output @@ -0,0 +1,29 @@ + +DEAL:0::numproc=4 +DEAL:0::0:first owned entry: 0.00000 +DEAL:0::0:second owned entry: 2.00000 +DEAL:0::0:ghost entry: 21.0000 +DEAL:0::0:ghost entry: 0.00000 +DEAL:0::0:ghost entry: -1.00000 +DEAL:0::OK + +DEAL:1::1:first owned entry: 4.00000 +DEAL:1::1:second owned entry: 6.00000 +DEAL:1::1:ghost entry: 21.0000 +DEAL:1::1:ghost entry: 0.00000 +DEAL:1::1:ghost entry: -1.00000 + + +DEAL:2::2:first owned entry: 8.00000 +DEAL:2::2:second owned entry: 10.0000 +DEAL:2::2:ghost entry: 21.0000 +DEAL:2::2:ghost entry: 0.00000 +DEAL:2::2:ghost entry: -1.00000 + + +DEAL:3::3:first owned entry: 12.0000 +DEAL:3::3:second owned entry: 14.0000 +DEAL:3::3:ghost entry: 21.0000 +DEAL:3::3:ghost entry: 0.00000 +DEAL:3::3:ghost entry: -1.00000 +