From: David Wells Date: Wed, 28 Dec 2016 17:39:22 +0000 (-0500) Subject: Add an 'add' method to the native parallel vectors. X-Git-Tag: v8.5.0-rc1~278^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2b6af4a08cee8dc90b31ebd8bfc586d80726df2b;p=dealii.git Add an 'add' method to the native parallel vectors. This is necessary for an upcoming patch that uses the method inside ConstraintMatrix::distribute_local_to_global. --- diff --git a/include/deal.II/lac/la_parallel_block_vector.h b/include/deal.II/lac/la_parallel_block_vector.h index f96e64052e..bbfeae8134 100644 --- a/include/deal.II/lac/la_parallel_block_vector.h +++ b/include/deal.II/lac/la_parallel_block_vector.h @@ -458,6 +458,13 @@ namespace LinearAlgebra virtual void add(const Number a, const VectorSpaceVector &V, const Number b, const VectorSpaceVector &W); + /** + * A collective add operation: This function adds a whole set of values + * stored in @p values to the vector components specified by @p indices. + */ + virtual void add (const std::vector &indices, + const std::vector &values); + /** * Scaling and simple addition of a multiple of a vector, i.e. *this = * s*(*this)+a*V. diff --git a/include/deal.II/lac/la_parallel_block_vector.templates.h b/include/deal.II/lac/la_parallel_block_vector.templates.h index 1f5ea732de..116c34d205 100644 --- a/include/deal.II/lac/la_parallel_block_vector.templates.h +++ b/include/deal.II/lac/la_parallel_block_vector.templates.h @@ -497,6 +497,18 @@ namespace LinearAlgebra + template + void + BlockVector::add (const std::vector &indices, + const std::vector &values) + { + for (size_type i=0; i bool BlockVector::all_zero () const diff --git a/include/deal.II/lac/la_parallel_vector.h b/include/deal.II/lac/la_parallel_vector.h index 43fa05a264..b73fd4ea7e 100644 --- a/include/deal.II/lac/la_parallel_vector.h +++ b/include/deal.II/lac/la_parallel_vector.h @@ -584,6 +584,13 @@ namespace LinearAlgebra virtual void add(const Number a, const VectorSpaceVector &V, const Number b, const VectorSpaceVector &W); + /** + * A collective add operation: This function adds a whole set of values + * stored in @p values to the vector components specified by @p indices. + */ + virtual void add (const std::vector &indices, + const std::vector &values); + /** * Scaling and simple addition of a multiple of a vector, i.e. *this = * s*(*this)+a*V. diff --git a/include/deal.II/lac/la_parallel_vector.templates.h b/include/deal.II/lac/la_parallel_vector.templates.h index 37d8433162..48a24d59da 100644 --- a/include/deal.II/lac/la_parallel_vector.templates.h +++ b/include/deal.II/lac/la_parallel_vector.templates.h @@ -1044,6 +1044,20 @@ namespace LinearAlgebra + template + void + Vector::add (const std::vector &indices, + const std::vector &values) + { + for (std::size_t i=0; ioperator()(indices[i]) += values[i]; + } + } + + + + template void Vector::sadd (const Number x, diff --git a/tests/mpi/parallel_block_vector_04.cc b/tests/mpi/parallel_block_vector_04.cc new file mode 100644 index 0000000000..d6eb2f5f20 --- /dev/null +++ b/tests/mpi/parallel_block_vector_04.cc @@ -0,0 +1,200 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 global reduction operations (norms, operator==, operator!=) with the +// parallel block vector class. This is the same as parallel_vector_01 except +// that this test uses the version of the add() method that takes a std::vector +// of indices and a std::vector of values to build the block vector instead of +// vector assignment. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include + + +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 from processor 1 to 8 + // owns 2 indices (the other processors do + // not own any dof), and all processors are + // ghosting element 1 (the second) + IndexSet local_owned(std::min(16U,numproc*2)); + if (myid < 8) + local_owned.add_range(myid*2,myid*2+2); + IndexSet local_relevant(numproc*2); + local_relevant = local_owned; + local_relevant.add_range(1,2); + + LinearAlgebra::distributed::Vector v(local_owned, local_relevant, + MPI_COMM_WORLD); + std::vector indices; + std::vector values; + + // set local values + if (myid < 8) + { + // unlike the first test, we scale values here since these are added to + // w directly + indices.push_back(myid*2); + values.push_back(myid*4.0); + indices.push_back(myid*2 + 1); + values.push_back(myid*4.0 + 2.0); + } + + LinearAlgebra::distributed::BlockVector w(3); + for (unsigned int i=0; i<3; ++i) + w.block(i) = v; + w.collect_sizes(); + + // set up v for comparison purposes below + v.add(indices, values); + v.compress(VectorOperation::add); + // in analogy to the first test, v is scaled differently than w + + if (myid < 8) + { + for (unsigned int i = 0; i<3; ++i) + { + // update the values in w + w.add(indices, values); + // update indices to point to the next block + indices[0] += v.size(); + indices[1] += v.size(); + } + } + w.compress(VectorOperation::add); + + // check l2 norm + { + const double l2_norm = w.l2_norm(); + if (myid == 0) + deallog << "l2 norm: " << l2_norm << std::endl; + AssertThrow (std::abs(v.l2_norm()*std::sqrt(3.)-w.l2_norm()) < 1e-13, + ExcInternalError()); + } + + // check l1 norm + { + const double l1_norm = w.l1_norm(); + if (myid == 0) + deallog << "l1 norm: " << l1_norm << std::endl; + AssertThrow (std::abs(v.l1_norm()*3.-w.l1_norm()) < 1e-14, + ExcInternalError()); + } + + // check linfty norm + { + const double linfty_norm = w.linfty_norm(); + if (myid == 0) + deallog << "linfty norm: " << linfty_norm << std::endl; + AssertThrow (v.linfty_norm()==w.linfty_norm(), + ExcInternalError()); + } + + // check lp norm + { + const double lp_norm = w.lp_norm(2.2); + if (myid == 0) + deallog << "l2.2 norm: " << lp_norm << std::endl; + + AssertThrow (std::fabs (w.l2_norm() - w.lp_norm(2.0)) < 1e-14, + ExcInternalError()); + } + + // check mean value (should be equal to l1 + // norm divided by vector size here since we + // have no negative entries) + { + const double mean = w.mean_value(); + if (myid == 0) + deallog << "Mean value: " << mean << std::endl; + + AssertThrow (std::fabs (mean * w.size() - w.l1_norm()) < 1e-15, + ExcInternalError()); + } + // check inner product + { + const double norm_sqr = w.l2_norm() * w.l2_norm(); + AssertThrow (std::fabs(w * w - norm_sqr) < 1e-12, + ExcInternalError()); + LinearAlgebra::distributed::BlockVector w2; + w2 = w; + AssertThrow (std::fabs(w2 * w - norm_sqr) < 1e-12, + ExcInternalError()); + + if (myid<8) + w2.block(0).local_element(0) = -1; + const double inner_prod = w * w2; + if (myid == 0) + deallog << "Inner product: " << inner_prod << std::endl; + } + + // check all_zero + { + bool allzero = w.all_zero(); + if (myid == 0) + deallog << " v==0 ? " << allzero << std::endl; + LinearAlgebra::distributed::BlockVector w2; + w2.reinit (w); + allzero = w2.all_zero(); + if (myid == 0) + deallog << " v2==0 ? " << allzero << std::endl; + + // now change one element to nonzero + if (myid == 0) + w2.block(1).local_element(1) = 1; + allzero = w2.all_zero(); + if (myid == 0) + deallog << " v2==0 ? " << allzero << 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()); + + 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); + deallog.threshold_double(1.e-10); + + test(); + } + else + test(); + +} diff --git a/tests/mpi/parallel_block_vector_04.mpirun=1.output b/tests/mpi/parallel_block_vector_04.mpirun=1.output new file mode 100644 index 0000000000..3a3e5186c7 --- /dev/null +++ b/tests/mpi/parallel_block_vector_04.mpirun=1.output @@ -0,0 +1,12 @@ + +DEAL:0::numproc=1 +DEAL:0::l2 norm: 3.464 +DEAL:0::l1 norm: 6.000 +DEAL:0::linfty norm: 2.000 +DEAL:0::l2.2 norm: 3.295 +DEAL:0::Mean value: 1.000 +DEAL:0::Inner product: 12.00 +DEAL:0:: v==0 ? 0 +DEAL:0:: v2==0 ? 1 +DEAL:0:: v2==0 ? 0 +DEAL:0::OK diff --git a/tests/mpi/parallel_block_vector_04.mpirun=10.output b/tests/mpi/parallel_block_vector_04.mpirun=10.output new file mode 100644 index 0000000000..ba129a4d03 --- /dev/null +++ b/tests/mpi/parallel_block_vector_04.mpirun=10.output @@ -0,0 +1,12 @@ + +DEAL:0::numproc=10 +DEAL:0::l2 norm: 122.0 +DEAL:0::l1 norm: 720.0 +DEAL:0::linfty norm: 30.00 +DEAL:0::l2.2 norm: 104.6 +DEAL:0::Mean value: 15.00 +DEAL:0::Inner product: 1.253e+04 +DEAL:0:: v==0 ? 0 +DEAL:0:: v2==0 ? 1 +DEAL:0:: v2==0 ? 0 +DEAL:0::OK diff --git a/tests/mpi/parallel_block_vector_04.mpirun=4.output b/tests/mpi/parallel_block_vector_04.mpirun=4.output new file mode 100644 index 0000000000..9ffa1364ef --- /dev/null +++ b/tests/mpi/parallel_block_vector_04.mpirun=4.output @@ -0,0 +1,12 @@ + +DEAL:0::numproc=4 +DEAL:0::l2 norm: 40.99 +DEAL:0::l1 norm: 168.0 +DEAL:0::linfty norm: 14.00 +DEAL:0::l2.2 norm: 36.31 +DEAL:0::Mean value: 7.000 +DEAL:0::Inner product: 1432. +DEAL:0:: v==0 ? 0 +DEAL:0:: v2==0 ? 1 +DEAL:0:: v2==0 ? 0 +DEAL:0::OK diff --git a/tests/mpi/parallel_vector_19.cc b/tests/mpi/parallel_vector_19.cc new file mode 100644 index 0000000000..fcd5f0e2e9 --- /dev/null +++ b/tests/mpi/parallel_vector_19.cc @@ -0,0 +1,96 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// use same partition as in parallel_vector_06 to check the version of the +// add() method that takes a std::vector of indices and a std::vector of +// values. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include + + +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 from processor 1 to 8 + // owns 2 indices (the other processors do + // not own any dof), and all processors are + // ghosting element 1 (the second) + IndexSet local_owned(std::min(16U,numproc*2)); + if (myid < 8) + local_owned.add_range(myid*2,myid*2+2); + IndexSet local_relevant(numproc*2); + local_relevant = local_owned; + local_relevant.add_range(1,2); + + LinearAlgebra::distributed::Vector v(local_owned, local_owned, MPI_COMM_WORLD); + + // set local values + if (myid < 8) + { + types::global_dof_index n_elements = 2; + std::vector indices(n_elements); + indices[0] = myid*2; + indices[1] = myid*2+1; + std::vector values(n_elements); + values[0] = myid*2.0; + values[1] = myid*2.0+1.0; + v.add (indices, values); + } + v.compress(VectorOperation::insert); + v*=2.0; + if (myid < 8) + { + AssertThrow(v(myid*2) == myid*4.0, ExcInternalError()); + AssertThrow(v(myid*2+1) == myid*4.0+2.0, ExcInternalError()); + } + + 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); + deallog.threshold_double(1.e-10); + + test(); + } + else + test(); + +} diff --git a/tests/mpi/parallel_vector_19.mpirun=1.output b/tests/mpi/parallel_vector_19.mpirun=1.output new file mode 100644 index 0000000000..bb1593da60 --- /dev/null +++ b/tests/mpi/parallel_vector_19.mpirun=1.output @@ -0,0 +1,3 @@ + +DEAL:0::numproc=1 +DEAL:0::OK diff --git a/tests/mpi/parallel_vector_19.mpirun=10.output b/tests/mpi/parallel_vector_19.mpirun=10.output new file mode 100644 index 0000000000..999baf1a98 --- /dev/null +++ b/tests/mpi/parallel_vector_19.mpirun=10.output @@ -0,0 +1,3 @@ + +DEAL:0::numproc=10 +DEAL:0::OK diff --git a/tests/mpi/parallel_vector_19.mpirun=4.output b/tests/mpi/parallel_vector_19.mpirun=4.output new file mode 100644 index 0000000000..6f5d2775d9 --- /dev/null +++ b/tests/mpi/parallel_vector_19.mpirun=4.output @@ -0,0 +1,3 @@ + +DEAL:0::numproc=4 +DEAL:0::OK