--- /dev/null
+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.
+<br>
+(Svenja Schoeder, 2018/06/18)
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 <typename Number>
+ Number
+ get_min(const Number a, const Number b)
+ {
+ return std::min(a, b);
+ }
+
+ template <typename Number>
+ std::complex<Number>
+ get_min(const std::complex<Number> a, const std::complex<Number>)
+ {
+ AssertThrow(false,
+ ExcMessage("VectorOperation::min max not"
+ "implemented for complex numbers"));
+ return a;
+ }
+
+ template <typename Number>
+ Number
+ get_max(const Number a, const Number b)
+ {
+ return std::max(a, b);
+ }
+
+ template <typename Number>
+ std::complex<Number>
+ get_max(const std::complex<Number> a, const std::complex<Number>)
+ {
+ AssertThrow(false,
+ ExcMessage("VectorOperation::min max not "
+ "implemented for complex numbers"));
+ return a;
+ }
} // namespace internal
// 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;
* @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
* 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;
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 <typename Number>
+ Number
+ get_min(const Number a, const Number b)
+ {
+ return std::min(a, b);
+ }
+ template <typename Number>
+ std::complex<Number>
+ get_min(const std::complex<Number> a, const std::complex<Number>)
+ {
+ AssertThrow(false,
+ ExcMessage("VectorOperation::min max not"
+ "implemented for complex numbers"));
+ return a;
+ }
+
+ template <typename Number>
+ Number
+ get_max(const Number a, const Number b)
+ {
+ return std::max(a, b);
+ }
+
+ template <typename Number>
+ std::complex<Number>
+ get_max(const std::complex<Number> a, const std::complex<Number>)
+ {
+ AssertThrow(false,
+ ExcMessage("VectorOperation::min max not "
+ "implemented for complex numbers"));
+ return a;
+ }
+ } // namespace internal
template <typename Number>
void
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));
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());
}
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<Number> tmp(n_elements);
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<Number> 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<Number> 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
/**
* The current operation is an addition.
*/
- add
+ add,
+ /**
+ * The current operation is a minimization.
+ */
+ min,
+ /**
+ * The current operation is a maximization.
+ */
+ max
};
};
cudaMemcpyHostToDevice);
AssertCuda(error_code);
}
- else
+ else if (operation == VectorOperation::add)
{
// Create a temporary vector on the device
Number * tmp;
error_code = cudaFree(tmp);
AssertCuda(error_code);
}
+ else
+ AssertThrow(false, ExcNotImplemented());
}
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());
}
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
{
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
{
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/index_set.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <iostream>
+#include <vector>
+
+#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<double> 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();
+}
--- /dev/null
+
+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
+
--- /dev/null
+
+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
+
projects / dealii.git / commitdiff