From 64063715f048fe7d6a6aabc29ff818ecefd73ead Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Mon, 15 Feb 2016 11:18:45 -0500 Subject: [PATCH] Make tests work again. --- include/deal.II/lac/la_vector.h | 7 ++++ include/deal.II/lac/la_vector.templates.h | 10 +++++ include/deal.II/lac/read_write_vector.h | 8 ++-- .../deal.II/lac/read_write_vector.templates.h | 11 ++--- include/deal.II/lac/trilinos_epetra_vector.h | 2 + include/deal.II/lac/vector_space_vector.h | 4 +- .../trilinos_epetra_communication_pattern.cc | 6 ++- source/lac/trilinos_epetra_vector.cc | 32 +++++++++++++-- tests/trilinos/epetra_vector_01.cc | 32 ++++++++------- tests/trilinos/epetra_vector_02.cc | 41 ++++++++++--------- tests/trilinos/readwritevector.cc | 29 ++++++++++++- 11 files changed, 131 insertions(+), 51 deletions(-) diff --git a/include/deal.II/lac/la_vector.h b/include/deal.II/lac/la_vector.h index 2ce6a28c82..fbdb0a75c7 100644 --- a/include/deal.II/lac/la_vector.h +++ b/include/deal.II/lac/la_vector.h @@ -132,6 +132,13 @@ namespace LinearAlgebra */ virtual Number operator* (const VectorSpaceVector &V) const override; + /** + * This function is not implemented and will throw an exception. + */ + virtual void import(const ReadWriteVector &V, + VectorOperation::values operation, + const CommunicationPatternBase *communication_pattern = NULL) override; + /** * Add @p a to all components. Note that @p a is a scalar not a vector. */ diff --git a/include/deal.II/lac/la_vector.templates.h b/include/deal.II/lac/la_vector.templates.h index bc014e9dae..a6ea7c6b52 100644 --- a/include/deal.II/lac/la_vector.templates.h +++ b/include/deal.II/lac/la_vector.templates.h @@ -114,6 +114,16 @@ namespace LinearAlgebra + template + void Vector::import(const ReadWriteVector &, + VectorOperation::values , + const CommunicationPatternBase *) + { + AssertThrow(false, ExcMessage("This function is not implemented.")); + } + + + template void Vector::add(const Number a, const VectorSpaceVector &V) { diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h index cb61c2c9a7..197ead3db0 100644 --- a/include/deal.II/lac/read_write_vector.h +++ b/include/deal.II/lac/read_write_vector.h @@ -233,7 +233,7 @@ namespace LinearAlgebra */ void import(const PETScWrappers::MPI::Vector &petsc_vec, VectorOperation::values operation, - const CommunicationPatternBase *communication_pattern = NULL); + const CommunicationPatternBase *communication_pattern = nullptr); #endif #ifdef DEAL_II_WITH_TRILINOS @@ -247,7 +247,7 @@ namespace LinearAlgebra */ void import(const TrilinosWrappers::MPI::Vector &trilinos_vec, VectorOperation::values operation, - const CommunicationPatternBase *communication_pattern = NULL); + const CommunicationPatternBase *communication_pattern = nullptr); /** * Imports all the elements present in the vector's IndexSet from the input @@ -259,7 +259,7 @@ namespace LinearAlgebra */ void import(const EpetraWrappers::Vector &epetra_vec, VectorOperation::values operation, - const CommunicationPatternBase *communication_pattern = NULL); + const CommunicationPatternBase *communication_pattern = nullptr); #endif /** @@ -446,7 +446,7 @@ namespace LinearAlgebra protected: #ifdef DEAL_II_WITH_TRILINOS /** - * Imports all the elements present in the vector's IndexSet from the input + * Import all the elements present in the vector's IndexSet from the input * vector @p multivector. This is an helper function and it should not be * used directly. */ diff --git a/include/deal.II/lac/read_write_vector.templates.h b/include/deal.II/lac/read_write_vector.templates.h index e3a7ea8a56..26c0c360de 100644 --- a/include/deal.II/lac/read_write_vector.templates.h +++ b/include/deal.II/lac/read_write_vector.templates.h @@ -191,10 +191,11 @@ namespace LinearAlgebra // The first time import is called, we create a communication pattern. // Check if the communication pattern already exists and if it can be // reused. - if (source_elements == source_stored_elements) + if ((source_elements.size() == source_stored_elements.size()) && + (source_elements == source_stored_elements)) { - epetra_comm_pattern.reset( - dynamic_cast (comm_pattern.get())); + epetra_comm_pattern = + std::dynamic_pointer_cast (comm_pattern); if (epetra_comm_pattern == nullptr) epetra_comm_pattern = std::make_shared( create_epetra_comm_pattern(source_elements, mpi_comm)); @@ -323,8 +324,8 @@ namespace LinearAlgebra source_stored_elements = source_index_set; EpetraWrappers::CommunicationPattern epetra_comm_pattern( source_stored_elements, stored_elements, mpi_comm); - comm_pattern = std::make_shared( - epetra_comm_pattern); + comm_pattern.reset(new EpetraWrappers::CommunicationPattern( + source_stored_elements, stored_elements, mpi_comm)); return epetra_comm_pattern; } diff --git a/include/deal.II/lac/trilinos_epetra_vector.h b/include/deal.II/lac/trilinos_epetra_vector.h index 67fd2efe71..b46e053009 100644 --- a/include/deal.II/lac/trilinos_epetra_vector.h +++ b/include/deal.II/lac/trilinos_epetra_vector.h @@ -202,6 +202,8 @@ namespace LinearAlgebra * calling separate methods means to load the calling vector @p this * twice. Since most vector operations are memory transfer limited, this * reduces the time by 25\% (or 50\% if @p W equals @p this). + * + * The vectors need to have the same layout. */ virtual double add_and_dot(const double a, const VectorSpaceVector &V, diff --git a/include/deal.II/lac/vector_space_vector.h b/include/deal.II/lac/vector_space_vector.h index 9f1f57ddc5..577da9b9d2 100644 --- a/include/deal.II/lac/vector_space_vector.h +++ b/include/deal.II/lac/vector_space_vector.h @@ -70,7 +70,7 @@ namespace LinearAlgebra virtual VectorSpaceVector &operator-= (const VectorSpaceVector &V) = 0; /** - * Imports all the elements present in the vector's IndexSet from the input + * Import all the elements present in the vector's IndexSet from the input * vector @p V. VectorOperation::values @p operation is used to decide if * the elements in @p V should be added to the current vector or replace the * current elements. The last parameter can be used if the same @@ -79,7 +79,7 @@ namespace LinearAlgebra */ virtual void import(const ReadWriteVector &V, VectorOperation::values operation, - const CommunicationPatternBase *communication_pattern = NULL); + const CommunicationPatternBase *communication_pattern = NULL) = 0; /** * Return the scalar product of two vectors. diff --git a/source/lac/trilinos_epetra_communication_pattern.cc b/source/lac/trilinos_epetra_communication_pattern.cc index ef93fc9db1..992d2d9c68 100644 --- a/source/lac/trilinos_epetra_communication_pattern.cc +++ b/source/lac/trilinos_epetra_communication_pattern.cc @@ -43,8 +43,10 @@ namespace LinearAlgebra { comm = std::make_shared(communicator); - Epetra_Map vector_space_vector_map = vector_space_vector_index_set.make_trilinos_map(*comm); - Epetra_Map read_write_vector_map = read_write_vector_index_set.make_trilinos_map(*comm); + Epetra_Map vector_space_vector_map = vector_space_vector_index_set.make_trilinos_map(*comm, + false); + Epetra_Map read_write_vector_map = read_write_vector_index_set.make_trilinos_map(*comm, + true); // Target map is read_write_vector_map // Source map is vector_space_vector_map. This map must have uniquely diff --git a/source/lac/trilinos_epetra_vector.cc b/source/lac/trilinos_epetra_vector.cc index b9ef42e3ae..b2b4cab8c6 100644 --- a/source/lac/trilinos_epetra_vector.cc +++ b/source/lac/trilinos_epetra_vector.cc @@ -33,7 +33,7 @@ namespace LinearAlgebra { Vector::Vector() : - vector(new Epetra_FEVector(Epetra_Map(0,1,0,Epetra_MpiComm(MPI_COMM_SELF)))), + vector(new Epetra_FEVector(Epetra_Map(0,0,0,Utilities::Trilinos::comm_self()))), epetra_comm_pattern(nullptr) {} @@ -110,7 +110,9 @@ namespace LinearAlgebra // The first time import is called, a communication pattern is created. // Check if the communication pattern already exists and if it can be // reused. - if (source_stored_elements == V.get_stored_elements()) + if ((source_stored_elements.size() != V.get_stored_elements().size()) || + ((source_stored_elements.size() == V.get_stored_elements().size()) && + (source_stored_elements != V.get_stored_elements()))) { create_epetra_comm_pattern(V.get_stored_elements(), dynamic_cast(vector->Comm()).Comm()); @@ -414,8 +416,30 @@ namespace LinearAlgebra ::dealii::IndexSet Vector::locally_owned_elements() const { - const Epetra_Map *map = dynamic_cast(&(vector->Map())); - return IndexSet(*map); + IndexSet is (size()); + + // easy case: local range is contiguous + if (vector->Map().LinearMap()) + { +#ifndef DEAL_II_WITH_64BIT_INDICES + is.add_range(vector->Map().MinMyGID(), vector->Map().MaxMyGID()+1); +#else + is.add_range(vector->Map().MinMyGID64(), vector->Map().MaxMyGID64()+1); +#endif + } + 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; } diff --git a/tests/trilinos/epetra_vector_01.cc b/tests/trilinos/epetra_vector_01.cc index 3059910224..1192ef5382 100644 --- a/tests/trilinos/epetra_vector_01.cc +++ b/tests/trilinos/epetra_vector_01.cc @@ -86,16 +86,18 @@ void test() } } - a = read_write_2; - b = read_write_1; - c = read_write_2; + a.import(read_write_2, VectorOperation::insert); + b.import(read_write_1, VectorOperation::insert); + c.import(read_write_2, VectorOperation::insert); - read_write_3 = a; + read_write_3.import(a, VectorOperation::insert); if (rank==0) { for (unsigned int i=0; i<5; ++i) - AssertThrow(read_write_2[i]==read_write_3[i], - ExcMessage("Vector a has been modified.")); + { + AssertThrow(read_write_2[i]==read_write_3[i], + ExcMessage("Vector a has been modified.")); + } } else { @@ -104,7 +106,7 @@ void test() ExcMessage("Vector a has been modified.")); } - read_write_3 = b; + read_write_3.import(b, VectorOperation::insert); if (rank==0) { for (unsigned int i=0; i<5; ++i) @@ -118,7 +120,7 @@ void test() ExcMessage("Vector b has been modified.")); } - read_write_3 = c; + read_write_3.import(c, VectorOperation::insert); if (rank==0) { for (unsigned int i=0; i<5; ++i) @@ -134,7 +136,7 @@ void test() a *= 2; - read_write_3 = a; + read_write_3.import(a, VectorOperation::insert); if (rank==0) { for (unsigned int i=0; i<5; ++i) @@ -149,7 +151,7 @@ void test() } c /= 2.; - read_write_3 = c; + read_write_3.import(c, VectorOperation::insert); if (rank==0) { for (unsigned int i=0; i<5; ++i) @@ -164,7 +166,7 @@ void test() } b += a; - read_write_3 = b; + read_write_3.import(b, VectorOperation::insert); if (rank==0) { for (unsigned int i=0; i<5; ++i) @@ -179,7 +181,7 @@ void test() } b -= c; - read_write_3 = b; + read_write_3.import(b, VectorOperation::insert); if (rank==0) { for (unsigned int i=0; i<5; ++i) @@ -193,8 +195,8 @@ void test() ExcMessage("Problem in operator -=.")); } - b = read_write_1; - c = read_write_1; + b.import(read_write_1, VectorOperation::insert); + c.import(read_write_1, VectorOperation::insert); const double val = b*c; AssertThrow(val==285., ExcMessage("Problem in operator *.")); } @@ -209,6 +211,8 @@ int main(int argc, char **argv) Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + test(); + deallog << "OK" < readwrite(readwrite_is); - readwrite = tril_vector; + readwrite.import(tril_vector,VectorOperation::insert); if (rank==0) { std::vector comp(4); @@ -78,6 +78,33 @@ void test() AssertThrow(readwrite.local_element(i)==comp[i], ExcMessage("Element not copied correctly")); } + + readwrite.import(tril_vector,VectorOperation::add); + + if (rank==0) + { + std::vector comp(4); + comp[0] = 0.; + comp[1] = 2.; + comp[2] = 12.; + comp[3] = 14.; + for (unsigned int i=0; i<4; ++i) + AssertThrow(readwrite.local_element(i)==comp[i], + ExcMessage("Element not copied correctly")); + } + MPI_Barrier(MPI_COMM_WORLD); + if (rank==1) + { + std::vector comp(4); + comp[0] = 4.; + comp[1] = 6.; + comp[2] = 8.; + comp[3] = 10.; + for (unsigned int i=0; i<4; ++i) + AssertThrow(readwrite.local_element(i)==comp[i], + ExcMessage("Element not copied correctly")); + } + deallog << "OK" <