From f94596e7f62faf15afcfa10972a45de7badca4ce Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Mon, 15 Aug 2016 17:30:09 +0200 Subject: [PATCH] Fix TrilinosWrappers::locally_owned_elements --- include/deal.II/lac/trilinos_vector.h | 7 +- include/deal.II/lac/trilinos_vector_base.h | 29 ++---- source/lac/trilinos_vector.cc | 93 +++++++++++++------ source/lac/trilinos_vector_base.cc | 6 +- .../mpi/trilinos_locally_owned_elements_02.cc | 88 ++++++++++++++++++ ...ents_02.with_trilinos=true.mpirun=2.output | 17 ++++ 6 files changed, 186 insertions(+), 54 deletions(-) create mode 100644 tests/mpi/trilinos_locally_owned_elements_02.cc create mode 100644 tests/mpi/trilinos_locally_owned_elements_02.with_trilinos=true.mpirun=2.output diff --git a/include/deal.II/lac/trilinos_vector.h b/include/deal.II/lac/trilinos_vector.h index ef4cbb8783..8e57c44931 100644 --- a/include/deal.II/lac/trilinos_vector.h +++ b/include/deal.II/lac/trilinos_vector.h @@ -661,6 +661,7 @@ namespace TrilinosWrappers { *this = Vector(parallel_partitioner.make_trilinos_map (communicator, true), v); + owned_elements = parallel_partitioner; } @@ -671,9 +672,7 @@ namespace TrilinosWrappers const dealii::Vector &v) { if (vector.get() == 0 || vector->Map().SameAs(parallel_partitioner) == false) - vector.reset (new Epetra_FEVector(parallel_partitioner)); - - has_ghosts = vector->Map().UniqueGIDs()==false; + reinit(parallel_partitioner); const int size = parallel_partitioner.NumMyElements(); @@ -681,6 +680,7 @@ namespace TrilinosWrappers // that a direct access is not possible. for (int i=0; i nonlocal_vector; + /** + * An IndexSet storing the indices this vector owns exclusively. + */ + IndexSet owned_elements; + /** * Make the reference class a friend. */ @@ -1082,27 +1087,7 @@ namespace TrilinosWrappers IndexSet VectorBase::locally_owned_elements() const { - IndexSet is (size()); - - // easy case: local range is contiguous - if (vector->Map().LinearMap()) - { - const std::pair x = local_range(); - is.add_range (x.first, x.second); - } - else if (vector->Map().NumMyElements() > 0) - { - const size_type n_indices = vector->Map().NumMyElements(); -#ifndef DEAL_II_WITH_64BIT_INDICES - unsigned int *vector_indices = (unsigned int *)vector->Map().MyGlobalElements(); -#else - size_type *vector_indices = (size_type *)vector->Map().MyGlobalElements64(); -#endif - is.add_indices(vector_indices, vector_indices+n_indices); - is.compress(); - } - - return is; + return owned_elements; } @@ -1218,6 +1203,8 @@ namespace TrilinosWrappers if (v.nonlocal_vector.get() != 0) nonlocal_vector.reset(new Epetra_MultiVector(v.nonlocal_vector->Map(), 1)); + + owned_elements = v.locally_owned_elements(); } diff --git a/source/lac/trilinos_vector.cc b/source/lac/trilinos_vector.cc index 9bdddcb538..7c3f6868e3 100644 --- a/source/lac/trilinos_vector.cc +++ b/source/lac/trilinos_vector.cc @@ -122,6 +122,7 @@ namespace TrilinosWrappers last_action = Zero; vector.reset (new Epetra_FEVector(*v.vector)); has_ghosts = v.has_ghosts; + owned_elements = v.locally_owned_elements(); } @@ -132,6 +133,7 @@ namespace TrilinosWrappers // initialize a minimal, valid object and swap last_action = Zero; vector.reset(new Epetra_FEVector(Epetra_Map(0,0,0,Utilities::Trilinos::comm_self()))); + owned_elements.clear(); swap(v); } @@ -186,9 +188,7 @@ namespace TrilinosWrappers : VectorBase() { - IndexSet parallel_partitioning = local; - parallel_partitioning.add_indices(ghost); - reinit(parallel_partitioning, communicator); + reinit(local, ghost, communicator, false); } @@ -207,6 +207,30 @@ namespace TrilinosWrappers vector.reset (new Epetra_FEVector(input_map)); has_ghosts = vector->Map().UniqueGIDs()==false; + + { + owned_elements.clear(); + owned_elements.set_size(size()); + + // easy case: local range is contiguous + if (vector->Map().LinearMap()) + { + const std::pair x = local_range(); + owned_elements.add_range (x.first, x.second); + } + else if (vector->Map().NumMyElements() > 0) + { + const size_type n_indices = vector->Map().NumMyElements(); +#ifndef DEAL_II_WITH_64BIT_INDICES + unsigned int *vector_indices = (unsigned int *)vector->Map().MyGlobalElements(); +#else + size_type *vector_indices = (size_type *)vector->Map().MyGlobalElements64(); +#endif + owned_elements.add_indices(vector_indices, vector_indices+n_indices); + owned_elements.compress(); + } + } + last_action = Zero; } @@ -215,13 +239,20 @@ namespace TrilinosWrappers void Vector::reinit (const IndexSet ¶llel_partitioner, const MPI_Comm &communicator, - const bool omit_zeroing_entries) + const bool /*omit_zeroing_entries*/) { nonlocal_vector.reset(); Epetra_Map map = parallel_partitioner.make_trilinos_map (communicator, true); - reinit (map, omit_zeroing_entries); + + vector.reset (new Epetra_FEVector(map)); + + has_ghosts = vector->Map().UniqueGIDs()==false; + + owned_elements = parallel_partitioner; + + last_action = Zero; } @@ -255,6 +286,7 @@ namespace TrilinosWrappers vector.reset (new Epetra_FEVector(v.vector->Map())); has_ghosts = v.has_ghosts; last_action = Zero; + owned_elements = v.locally_owned_elements(); } else if (omit_zeroing_entries == false) { @@ -303,6 +335,7 @@ namespace TrilinosWrappers const bool import_data) { nonlocal_vector.reset(); + owned_elements.clear(); // In case we do not allow to have different maps, this call means that // we have to reset the vector. So clear the vector, initialize our map @@ -316,6 +349,7 @@ namespace TrilinosWrappers for (size_type block=0; block global_ids (n_elements, -1); + size_type max_size = 0; for (size_type block=0; blockMap().UniqueGIDs()==false; + + last_action = Zero; } else { @@ -438,6 +480,7 @@ namespace TrilinosWrappers vector.reset (new Epetra_FEVector(*v.vector)); last_action = Zero; has_ghosts = v.has_ghosts; + owned_elements = v.locally_owned_elements(); } if (v.nonlocal_vector.get() != 0) @@ -519,20 +562,14 @@ namespace TrilinosWrappers Vector::Vector (const size_type n) { - last_action = Zero; - Epetra_LocalMap map ((TrilinosWrappers::types::int_type)n, 0, Utilities::Trilinos::comm_self()); - vector.reset (new Epetra_FEVector (map)); + reinit(n); } Vector::Vector (const Epetra_Map &input_map) { - last_action = Zero; - Epetra_LocalMap map (n_global_elements(input_map), - input_map.IndexBase(), - input_map.Comm()); - vector.reset (new Epetra_FEVector(map)); + reinit(input_map); } @@ -540,16 +577,7 @@ namespace TrilinosWrappers Vector::Vector (const IndexSet &partitioning, const MPI_Comm &communicator) { - last_action = Zero; - Epetra_LocalMap map (static_cast(partitioning.size()), - 0, -#ifdef DEAL_II_WITH_MPI - Epetra_MpiComm(communicator)); -#else - Epetra_SerialComm()); - (void)communicator; -#endif - vector.reset (new Epetra_FEVector(map)); + reinit (partitioning, communicator); } @@ -569,7 +597,6 @@ namespace TrilinosWrappers } else reinit (v, false, true); - } @@ -578,6 +605,9 @@ namespace TrilinosWrappers Vector::reinit (const size_type n, const bool /*omit_zeroing_entries*/) { + owned_elements.clear(); + owned_elements.set_size(n); + owned_elements.add_range(0, n); Epetra_LocalMap map ((TrilinosWrappers::types::int_type)n, 0, Utilities::Trilinos::comm_self()); vector.reset (new Epetra_FEVector (map)); @@ -594,7 +624,9 @@ namespace TrilinosWrappers Epetra_LocalMap map (n_global_elements(input_map), input_map.IndexBase(), input_map.Comm()); - vector.reset (new Epetra_FEVector (map)); + vector.reset (new Epetra_FEVector(map)); + owned_elements.set_size(n_global_elements(input_map)); + owned_elements.add_range(0, n_global_elements(input_map)); last_action = Zero; } @@ -617,6 +649,7 @@ namespace TrilinosWrappers vector.reset (new Epetra_FEVector(map)); last_action = Zero; + owned_elements = partitioning; } @@ -650,6 +683,7 @@ namespace TrilinosWrappers v.vector->Map().IndexBase(), v.vector->Comm()); vector.reset (new Epetra_FEVector(map)); + owned_elements = v.locally_owned_elements(); } else if (omit_zeroing_entries) { @@ -723,6 +757,7 @@ namespace TrilinosWrappers v.vector->Map().IndexBase(), v.vector->Comm()); vector.reset (new Epetra_FEVector(map)); + owned_elements = v.locally_owned_elements(); } const int ierr = vector->Update(1.0, *v.vector, 0.0); diff --git a/source/lac/trilinos_vector_base.cc b/source/lac/trilinos_vector_base.cc index fd05084a02..1d7ebe5899 100644 --- a/source/lac/trilinos_vector_base.cc +++ b/source/lac/trilinos_vector_base.cc @@ -101,7 +101,10 @@ namespace TrilinosWrappers compressed (true), has_ghosts (v.has_ghosts), vector(new Epetra_FEVector(*v.vector)) - {} + { + owned_elements.set_size(v.size()); + owned_elements.add_range(0,v.size()); + } @@ -139,6 +142,7 @@ namespace TrilinosWrappers last_action = Zero; vector.reset (new Epetra_FEVector(*v.vector)); has_ghosts = v.has_ghosts; + owned_elements = v.locally_owned_elements(); } else { diff --git a/tests/mpi/trilinos_locally_owned_elements_02.cc b/tests/mpi/trilinos_locally_owned_elements_02.cc new file mode 100644 index 0000000000..7cd1ace7dc --- /dev/null +++ b/tests/mpi/trilinos_locally_owned_elements_02.cc @@ -0,0 +1,88 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2009 - 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. +// +// --------------------------------------------------------------------- + + + +// test TrilinosVector::locally_owned_elements + + +#include "../tests.h" +#include +#include + +#include +#include + + +void test () +{ + const int n_proc = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + const int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + //All processes should own 10 entries + const int entries_per_process = 10; + + IndexSet locally_owned(entries_per_process*n_proc); + const int begin_index = my_id*entries_per_process; + const int end_index = (my_id+1)*entries_per_process; + locally_owned.add_range(begin_index, end_index); + + IndexSet locally_relevant(entries_per_process*n_proc); + const int local_begin = std::max(0, begin_index-entries_per_process/2); + const int local_end = entries_per_process*n_proc; + locally_relevant.add_range (local_begin, local_end); + + TrilinosWrappers::MPI::Vector ghosted, distributed; + distributed.reinit(locally_owned, MPI_COMM_WORLD); + ghosted.reinit (locally_owned, locally_relevant, MPI_COMM_WORLD); + + IndexSet locally_owned_elements_distributed = distributed.locally_owned_elements(); + IndexSet locally_owned_elements_ghosted = ghosted.locally_owned_elements(); + + const types::global_dof_index local_range_begin_ghosted = ghosted.local_range().first; + const types::global_dof_index local_range_end_ghosted = ghosted.local_range().second; + + const types::global_dof_index local_range_begin_distributed = distributed.local_range().first; + const types::global_dof_index local_range_end_distributed = distributed.local_range().second; + + deallog << "locally_owned_elements_distributed: "; + locally_owned_elements_distributed.print(deallog); + deallog << "locally_owned_elements_ghosted: "; + locally_owned_elements_ghosted.print(deallog); + deallog << "local_range_begin_ghosted: " + << local_range_begin_ghosted << std::endl; + deallog << "local_range_end_ghosted: " + << local_range_end_ghosted << std::endl; + deallog << "local_range_begin_distributed: " + << local_range_begin_distributed << std::endl; + deallog << "local_range_end_distributed: " + << local_range_end_distributed << std::endl; + + AssertThrow (locally_owned_elements_distributed == locally_owned_elements_ghosted, + ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + + MPILogInitAll log; + + test(); +} diff --git a/tests/mpi/trilinos_locally_owned_elements_02.with_trilinos=true.mpirun=2.output b/tests/mpi/trilinos_locally_owned_elements_02.with_trilinos=true.mpirun=2.output new file mode 100644 index 0000000000..ce5508c5b9 --- /dev/null +++ b/tests/mpi/trilinos_locally_owned_elements_02.with_trilinos=true.mpirun=2.output @@ -0,0 +1,17 @@ + +DEAL:0::locally_owned_elements_distributed: {[0,9]} +DEAL:0::locally_owned_elements_ghosted: {[0,9]} +DEAL:0::local_range_begin_ghosted: 0 +DEAL:0::local_range_end_ghosted: 20 +DEAL:0::local_range_begin_distributed: 0 +DEAL:0::local_range_end_distributed: 10 +DEAL:0::OK + +DEAL:1::locally_owned_elements_distributed: {[10,19]} +DEAL:1::locally_owned_elements_ghosted: {[10,19]} +DEAL:1::local_range_begin_ghosted: 5 +DEAL:1::local_range_end_ghosted: 20 +DEAL:1::local_range_begin_distributed: 10 +DEAL:1::local_range_end_distributed: 20 +DEAL:1::OK + -- 2.39.5