From cabeb1ccd1a4b186885fdfba0027fae3e7225bbd Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 21 Mar 2020 23:14:25 +0100 Subject: [PATCH] Add test --- include/deal.II/lac/la_sm_vector.h | 19 +- include/deal.II/lac/la_sm_vector.templates.h | 133 +++++++++--- source/lac/la_sm_vector.inst.in | 25 +-- tests/sm/vector_sm_01.cc | 2 +- tests/sm/vector_sm_01.mpirun=3.output | 1 + tests/sm/vector_sm_02.cc | 208 +++++++++++++++++++ tests/sm/vector_sm_02.mpirun=1.output | 10 + tests/sm/vector_sm_02.mpirun=3.output | 10 + 8 files changed, 362 insertions(+), 46 deletions(-) create mode 100644 tests/sm/vector_sm_01.mpirun=3.output create mode 100644 tests/sm/vector_sm_02.cc create mode 100644 tests/sm/vector_sm_02.mpirun=1.output create mode 100644 tests/sm/vector_sm_02.mpirun=3.output diff --git a/include/deal.II/lac/la_sm_vector.h b/include/deal.II/lac/la_sm_vector.h index 79105114ee..d911f3a532 100644 --- a/include/deal.II/lac/la_sm_vector.h +++ b/include/deal.II/lac/la_sm_vector.h @@ -715,7 +715,16 @@ namespace LinearAlgebra Vector::local_element( const size_type local_index) const { - Assert(false, ExcNotImplemented()); + Assert((std::is_same::value), + ExcMessage( + "This function is only implemented for the Host memory space")); + AssertIndexRange(local_index, + partitioner->local_size() + + partitioner->n_ghost_indices()); + // do not allow reading a vector which is not in ghost mode + Assert(local_index < local_size() || vector_is_ghosted == true, + ExcMessage("You tried to read a ghost element of this vector, " + "but it has not imported its ghost values.")); return data.values[local_index]; } @@ -726,7 +735,13 @@ namespace LinearAlgebra inline Number & Vector::local_element(const size_type local_index) { - Assert(false, ExcNotImplemented()); + Assert((std::is_same::value), + ExcMessage( + "This function is only implemented for the Host memory space")); + + AssertIndexRange(local_index, + partitioner->local_size() + + partitioner->n_ghost_indices()); return data.values[local_index]; } diff --git a/include/deal.II/lac/la_sm_vector.templates.h b/include/deal.II/lac/la_sm_vector.templates.h index 320d8248f5..54a8697c28 100644 --- a/include/deal.II/lac/la_sm_vector.templates.h +++ b/include/deal.II/lac/la_sm_vector.templates.h @@ -134,9 +134,34 @@ namespace LinearAlgebra const Vector &v, const bool omit_zeroing_entries) { - Assert(false, ExcNotImplemented()); - (void)v; - (void)omit_zeroing_entries; + clear_mpi_requests(); + Assert(v.partitioner.get() != nullptr, ExcNotInitialized()); + + // check whether the partitioners are + // different (check only if the are allocated + // differently, not if the actual data is + // different) + if (partitioner.get() != v.partitioner.get()) + { + partitioner = v.partitioner; + partitioner_sm = v.partitioner_sm; + const size_type new_allocated_size = + partitioner->local_size() + partitioner->n_ghost_indices(); + resize_val(new_allocated_size, + partitioner_sm->get_sm_mpi_communicator()); + } + + if (omit_zeroing_entries == false) + this->operator=(Number()); + else + zero_out_ghosts(); + + // do not reallocate import_data directly, but only upon request. It + // is only used as temporary storage for compress() and + // update_ghost_values, and we might have vectors where we never + // call these methods and hence do not need to have the storage. + import_data.values.reset(); + import_data.values_dev.reset(); } @@ -215,8 +240,11 @@ namespace LinearAlgebra , allocated_size(0) , vector_is_ghosted(false) { - Assert(false, ExcNotImplemented()); - (void)v; + reinit(v, true); + + const size_type this_size = local_size(); + if (this_size > 0) + std::memcpy(this->begin(), v.begin(), this_size * sizeof(Number)); } @@ -489,8 +517,22 @@ namespace LinearAlgebra Vector:: operator-=(const VectorSpaceVector &vv) { - Assert(false, ExcNotImplemented()); - (void)vv; + // Downcast. Throws an exception if invalid. + using VectorType = Vector; + Assert(dynamic_cast(&vv) != nullptr, + ExcVectorTypeNotCompatible()); + const VectorType &v = dynamic_cast(vv); + + AssertDimension(local_size(), v.local_size()); + + auto values = this->begin(); + const auto values_other = v.begin(); + + for (unsigned int i = 0; i < partitioner->local_size(); i++) + values[i] -= values_other[i]; + + if (vector_is_ghosted) + update_ghost_values(); return *this; } @@ -513,9 +555,24 @@ namespace LinearAlgebra const Number a, const VectorSpaceVector &vv) { - Assert(false, ExcNotImplemented()); - (void)a; - (void)vv; + // Downcast. Throws an exception if invalid. + using VectorType = Vector; + Assert(dynamic_cast(&vv) != nullptr, + ExcVectorTypeNotCompatible()); + const VectorType &v = dynamic_cast(vv); + + AssertIsFinite(a); + AssertDimension(local_size(), v.local_size()); + + // nothing to do if a is zero + if (a == Number(0.)) + return; + + auto values = this->begin(); + const auto values_other = v.begin(); + + for (unsigned int i = 0; i < partitioner->local_size(); i++) + values[i] += a * values_other[i]; } @@ -525,9 +582,10 @@ namespace LinearAlgebra Vector::add(const Number a, const VectorSpaceVector &vv) { - Assert(false, ExcNotImplemented()); - (void)a; - (void)vv; + add_local(a, vv); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -580,10 +638,24 @@ namespace LinearAlgebra const Number a, const VectorSpaceVector &vv) { - Assert(false, ExcNotImplemented()); - (void)x; - (void)a; - (void)vv; + // Downcast. Throws an exception if invalid. + using VectorType = Vector; + Assert(dynamic_cast(&vv) != nullptr, + ExcVectorTypeNotCompatible()); + const VectorType &v = dynamic_cast(vv); + + AssertIsFinite(a); + AssertDimension(local_size(), v.local_size()); + + // nothing to do if a is zero + if (a == Number(0.)) + return; + + auto values = this->begin(); + const auto values_other = v.begin(); + + for (unsigned int i = 0; i < partitioner->local_size(); i++) + values[i] = x * values[i] + a * values_other[i]; } @@ -594,10 +666,10 @@ namespace LinearAlgebra const Number a, const VectorSpaceVector &vv) { - Assert(false, ExcNotImplemented()); - (void)x; - (void)a; - (void)vv; + sadd_local(x, a, vv); + + if (vector_is_ghosted) + update_ghost_values(); } @@ -813,8 +885,14 @@ namespace LinearAlgebra typename Vector::real_type Vector::linfty_norm_local() const { - Assert(false, ExcNotImplemented()); - return 0; + real_type max = 0.; + + auto values = this->begin(); + + for (unsigned int i = 0; i < partitioner->local_size(); i++) + max = std::max(std::abs(values[i]), max); + + return max; } @@ -823,8 +901,12 @@ namespace LinearAlgebra inline typename Vector::real_type Vector::linfty_norm() const { - Assert(false, ExcNotImplemented()); - return 0; + const real_type local_result = linfty_norm_local(); + if (partitioner->n_mpi_processes() > 1) + return Utilities::MPI::max(local_result, + partitioner->get_mpi_communicator()); + else + return local_result; } @@ -867,7 +949,6 @@ namespace LinearAlgebra Vector::partitioners_are_compatible( const Utilities::MPI::Partitioner &part) const { - Assert(false, ExcNotImplemented()); return partitioner->is_compatible(part); } diff --git a/source/lac/la_sm_vector.inst.in b/source/lac/la_sm_vector.inst.in index d54cca376d..0bd0109b65 100644 --- a/source/lac/la_sm_vector.inst.in +++ b/source/lac/la_sm_vector.inst.in @@ -27,26 +27,17 @@ for (SCALAR : REAL_SCALARS) ::dealii::MemorySpace::Host>( const Vector &, VectorOperation::values); - \} - \} - } - -for (S1 : REAL_SCALARS; S2 : REAL_SCALARS) - { - namespace LinearAlgebra - \{ - namespace SharedMPI - \{ template void - Vector::reinit( - const Vector &, + Vector::reinit( + const Vector &, const bool); - template S1 - Vector::inner_product_local( - const Vector &) const; + template SCALAR + Vector::inner_product_local< + SCALAR>(const Vector &) const; template void - Vector::copy_locally_owned_data_from< - S2>(const Vector &); + Vector:: + copy_locally_owned_data_from( + const Vector &); \} \} } diff --git a/tests/sm/vector_sm_01.cc b/tests/sm/vector_sm_01.cc index ed2ebca093..3a06dde010 100644 --- a/tests/sm/vector_sm_01.cc +++ b/tests/sm/vector_sm_01.cc @@ -60,7 +60,7 @@ test(const int n_refinements, const int degree, const int group_size) dealii::MatrixFree matrix_free; matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data); - MPI_Comm comm = MPI_COMM_WORLD; + const MPI_Comm comm = MPI_COMM_WORLD; const unsigned int rank = Utilities::MPI::this_mpi_process(comm); diff --git a/tests/sm/vector_sm_01.mpirun=3.output b/tests/sm/vector_sm_01.mpirun=3.output new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/tests/sm/vector_sm_01.mpirun=3.output @@ -0,0 +1 @@ + diff --git a/tests/sm/vector_sm_02.cc b/tests/sm/vector_sm_02.cc new file mode 100644 index 0000000000..4ee0d69cf3 --- /dev/null +++ b/tests/sm/vector_sm_02.cc @@ -0,0 +1,208 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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 the correctness of the matrix-free cell loop with additional +// operation_before_loop and operation_after_loop lambdas + +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + +#include "../matrix_free/matrix_vector_mf.h" + + + +template +class Matrix +{ +public: + Matrix(const MatrixFree &data_in) + : data(data_in) + {} + + void + vmult(VectorType &dst, const VectorType &src) const + { + const std::function &, + VectorType &, + const VectorType &, + const std::pair &)> + wrap = helmholtz_operator; + dst = 0; + data.cell_loop(wrap, dst, src); + for (auto i : data.get_constrained_dofs()) + dst.local_element(i) = src.local_element(i); + } + + void + vmult_with_update_basic(VectorType &dst, + VectorType &src, + VectorType &other_vector) const + { + src.add(1.5, other_vector); + other_vector.add(0.7, dst); + dst = 0; + vmult(dst, src); + dst.sadd(0.63, -1.3, other_vector); + } + + void + vmult_with_update_merged(VectorType &dst, + VectorType &src, + VectorType &other_vector) const + { + const std::function &, + VectorType &, + const VectorType &, + const std::pair &)> + wrap = helmholtz_operator; + data.cell_loop( + wrap, + dst, + src, + // operation before cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + src.local_element(i) += 1.5 * other_vector.local_element(i); + other_vector.local_element(i) += 0.7 * dst.local_element(i); + dst.local_element(i) = 0; + } + }, + // operation after cell operation + [&](const unsigned int start_range, const unsigned int end_range) { + for (unsigned int i = start_range; i < end_range; ++i) + { + dst.local_element(i) = + 0.63 * dst.local_element(i) - 1.3 * other_vector.local_element(i); + } + }); + } + +private: + const MatrixFree &data; +}; + + + +template +void +test() +{ + const MPI_Comm comm = MPI_COMM_WORLD; + + parallel::distributed::Triangulation tria(comm); + GridGenerator::hyper_cube(tria); + tria.refine_global(7 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + MatrixFree mf_data; + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + mf_data.reinit(dof, constraints, quad, data); + } + + const unsigned int group_size = 1; + const unsigned int rank = Utilities::MPI::this_mpi_process(comm); + + MPI_Comm comm_sm; + MPI_Comm_split(comm, rank / group_size, rank, &comm_sm); + + using VectorType = LinearAlgebra::SharedMPI::Vector; + + Matrix mf(mf_data); + VectorType vec1, vec2, vec3; + mf_data.initialize_dof_vector(vec1, comm_sm); + mf_data.initialize_dof_vector(vec2, comm_sm); + mf_data.initialize_dof_vector(vec3, comm_sm); + + for (unsigned int i = 0; i < vec1.local_size(); ++i) + { + double entry = random_value(); + vec1.local_element(i) = entry; + entry = random_value(); + vec2.local_element(i) = entry; + entry = random_value(); + vec3.local_element(i) = entry; + } + + VectorType ref1 = vec1; + VectorType ref2 = vec2; + VectorType ref3 = vec3; + + mf.vmult_with_update_basic(ref3, ref2, ref1); + mf.vmult_with_update_basic(vec3, vec2, vec1); + + vec3 -= ref3; + deallog << "Error in 1x merged operation: " << vec3.linfty_norm() + << std::endl; + + ref3 = 0; + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + mf.vmult_with_update_basic(ref1, ref2, ref3); + mf.vmult_with_update_basic(vec1, vec2, vec3); + + vec3 -= ref3; + deallog << "Error in 2x merged operation: " << vec3.linfty_norm() + << std::endl; +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); + deallog << std::setprecision(3); + + test<2, 3, double>(); + test<2, 3, float>(); + test<3, 3, double>(); +} \ No newline at end of file diff --git a/tests/sm/vector_sm_02.mpirun=1.output b/tests/sm/vector_sm_02.mpirun=1.output new file mode 100644 index 0000000000..57a015933a --- /dev/null +++ b/tests/sm/vector_sm_02.mpirun=1.output @@ -0,0 +1,10 @@ + +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<3>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 diff --git a/tests/sm/vector_sm_02.mpirun=3.output b/tests/sm/vector_sm_02.mpirun=3.output new file mode 100644 index 0000000000..57a015933a --- /dev/null +++ b/tests/sm/vector_sm_02.mpirun=3.output @@ -0,0 +1,10 @@ + +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<2>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 +DEAL::Testing FE_Q<3>(3) +DEAL::Error in 1x merged operation: 0.00 +DEAL::Error in 2x merged operation: 0.00 -- 2.39.5