From ba314705d507e7d318bf544c969c3799f03be5b2 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Mon, 27 Aug 2018 17:14:38 +0200 Subject: [PATCH] Implement TpetraWrappers::Vector --- cmake/config/template-arguments.in | 5 + cmake/configure/configure_2_trilinos.cmake | 3 +- doc/external-libs/trilinos.html | 2 + include/deal.II/base/index_set.h | 5 + include/deal.II/base/utilities.h | 3 + include/deal.II/dofs/dof_accessor.templates.h | 27 + .../fe/fe_tools_extrapolate.templates.h | 15 + .../fe/fe_tools_interpolate.templates.h | 16 + include/deal.II/lac/read_write_vector.h | 39 +- .../deal.II/lac/read_write_vector.templates.h | 121 +++- include/deal.II/lac/trilinos_sparse_matrix.h | 1 + .../trilinos_tpetra_communication_pattern.h | 106 +++ include/deal.II/lac/trilinos_tpetra_vector.h | 406 +++++++++++ include/deal.II/lac/vector_element_access.h | 70 +- include/deal.II/multigrid/mg_transfer.h | 32 + .../deal.II/numerics/vector_tools.templates.h | 1 + source/base/index_set.cc | 66 ++ source/base/time_stepping.cc | 1 + source/base/utilities.cc | 14 + source/dofs/dof_accessor_get.cc | 1 + source/dofs/dof_accessor_set.cc | 1 + source/fe/fe_values.cc | 1 + source/fe/mapping_fe_field.cc | 1 + source/lac/CMakeLists.txt | 2 + source/lac/solver.cc | 1 + source/lac/trilinos_sparse_matrix.cc | 50 ++ .../trilinos_tpetra_communication_pattern.cc | 104 +++ source/lac/trilinos_tpetra_vector.cc | 632 ++++++++++++++++++ source/lac/vector_memory.cc | 1 + source/multigrid/mg_base.cc | 1 + source/numerics/data_out_dof_data_codim.cc | 1 + source/numerics/derivative_approximation.cc | 1 + tests/trilinos/tpetra_vector_01.cc | 225 +++++++ .../trilinos/tpetra_vector_01.mpirun=2.output | 2 + tests/trilinos/tpetra_vector_02.cc | 216 ++++++ .../trilinos/tpetra_vector_02.mpirun=2.output | 2 + tests/trilinos/tpetra_vector_03.cc | 101 +++ .../trilinos/tpetra_vector_03.mpirun=2.output | 70 ++ .../trilinos/tpetra_vector_03.mpirun=4.output | 146 ++++ 39 files changed, 2488 insertions(+), 4 deletions(-) create mode 100644 include/deal.II/lac/trilinos_tpetra_communication_pattern.h create mode 100644 include/deal.II/lac/trilinos_tpetra_vector.h create mode 100644 source/lac/trilinos_tpetra_communication_pattern.cc create mode 100644 source/lac/trilinos_tpetra_vector.cc create mode 100644 tests/trilinos/tpetra_vector_01.cc create mode 100644 tests/trilinos/tpetra_vector_01.mpirun=2.output create mode 100644 tests/trilinos/tpetra_vector_02.cc create mode 100644 tests/trilinos/tpetra_vector_02.mpirun=2.output create mode 100644 tests/trilinos/tpetra_vector_03.cc create mode 100644 tests/trilinos/tpetra_vector_03.mpirun=2.output create mode 100644 tests/trilinos/tpetra_vector_03.mpirun=4.output diff --git a/cmake/config/template-arguments.in b/cmake/config/template-arguments.in index 551fdd08d3..17a9f2ba94 100644 --- a/cmake/config/template-arguments.in +++ b/cmake/config/template-arguments.in @@ -95,6 +95,7 @@ VECTOR_TYPES := { Vector; @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; @DEAL_II_EXPAND_EPETRA_VECTOR@; + @DEAL_II_EXPAND_TPETRA_VECTOR@; @DEAL_II_EXPAND_PETSC_MPI_VECTOR@; @DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR@; @@ -120,6 +121,7 @@ REAL_VECTOR_TYPES := { Vector; @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; @DEAL_II_EXPAND_EPETRA_VECTOR@; + @DEAL_II_EXPAND_TPETRA_VECTOR@; @DEAL_II_EXPAND_PETSC_MPI_VECTOR_REAL@; @DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR@; @@ -138,6 +140,7 @@ REAL_NONBLOCK_VECTORS := { Vector; @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; @DEAL_II_EXPAND_EPETRA_VECTOR@; + @DEAL_II_EXPAND_TPETRA_VECTOR@; @DEAL_II_EXPAND_PETSC_MPI_VECTOR_REAL@; } @@ -145,6 +148,7 @@ REAL_NONBLOCK_VECTORS := { Vector; EXTERNAL_PARALLEL_VECTORS := { @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; @DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR@; @DEAL_II_EXPAND_EPETRA_VECTOR@; + @DEAL_II_EXPAND_TPETRA_VECTOR@; @DEAL_II_EXPAND_PETSC_MPI_VECTOR@; @DEAL_II_EXPAND_PETSC_MPI_BLOCKVECTOR@ } @@ -163,6 +167,7 @@ VECTORS_WITH_MATRIX := { Vector; @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; @DEAL_II_EXPAND_EPETRA_VECTOR@; + @DEAL_II_EXPAND_TPETRA_VECTOR@; @DEAL_II_EXPAND_PETSC_MPI_VECTOR@; } diff --git a/cmake/configure/configure_2_trilinos.cmake b/cmake/configure/configure_2_trilinos.cmake index c53f105c09..abe66905d8 100644 --- a/cmake/configure/configure_2_trilinos.cmake +++ b/cmake/configure/configure_2_trilinos.cmake @@ -42,7 +42,7 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var) ) FOREACH(_module - Amesos Epetra Ifpack AztecOO Teuchos ML MueLu + Amesos Epetra Ifpack AztecOO Teuchos Tpetra ML MueLu ) ITEM_MATCHES(_module_found ${_module} ${Trilinos_PACKAGE_LIST}) IF(_module_found) @@ -227,6 +227,7 @@ MACRO(FEATURE_TRILINOS_CONFIGURE_EXTERNAL) SET(DEAL_II_EXPAND_TRILINOS_MPI_VECTOR "TrilinosWrappers::MPI::Vector") IF (TRILINOS_WITH_MPI) SET(DEAL_II_EXPAND_EPETRA_VECTOR "LinearAlgebra::EpetraWrappers::Vector") + SET(DEAL_II_EXPAND_TPETRA_VECTOR "LinearAlgebra::TpetraWrappers::Vector") ENDIF() IF(${DEAL_II_TRILINOS_WITH_SACADO}) # Note: Only CMake 3.0 and greater support line continuation with the "\" character diff --git a/doc/external-libs/trilinos.html b/doc/external-libs/trilinos.html index 9579c9a62a..d950f55b11 100644 --- a/doc/external-libs/trilinos.html +++ b/doc/external-libs/trilinos.html @@ -67,6 +67,7 @@
  • ROL (optional),
  • Sacado (optional),
  • Teuchos, +
  • Tpetra,
  • Zoltan (optional). @@ -92,6 +93,7 @@ -DTrilinos_ENABLE_MueLu=ON \ -DTrilinos_ENABLE_ML=ON \ -DTrilinos_ENABLE_ROL=ON \ + -DTrilinos_ENABLE_Tpetra=ON \ -DTrilinos_ENABLE_Zoltan=ON \ -DTrilinos_VERBOSE_CONFIGURE=OFF \ -DTPL_ENABLE_MPI=ON \ diff --git a/include/deal.II/base/index_set.h b/include/deal.II/base/index_set.h index fec38b883f..2e21fab644 100644 --- a/include/deal.II/base/index_set.h +++ b/include/deal.II/base/index_set.h @@ -31,6 +31,7 @@ #ifdef DEAL_II_WITH_TRILINOS # include +# include #endif #if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC) @@ -446,6 +447,10 @@ public: Epetra_Map make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD, const bool overlapping = false) const; + + Tpetra::Map<> + make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD, + const bool overlapping = false) const; #endif diff --git a/include/deal.II/base/utilities.h b/include/deal.II/base/utilities.h index f49c302132..71d05faf24 100644 --- a/include/deal.II/base/utilities.h +++ b/include/deal.II/base/utilities.h @@ -888,6 +888,9 @@ namespace Utilities const Epetra_Comm & comm_self(); + const Teuchos::RCP> & + tpetra_comm_self(); + /** * Given a communicator, duplicate it. If the given communicator is * serial, that means to just return a copy of itself. On the other hand, diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index a7c898c953..73b77c07a5 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -1534,6 +1534,33 @@ namespace internal + template + static void + extract_subvector_to(const LinearAlgebra::TpetraWrappers::Vector &values, + const types::global_dof_index *cache_begin, + const types::global_dof_index *cache_end, + ForwardIterator local_values_begin) + { + std::vector sorted_indices_pos = + sort_indices(cache_begin, cache_end); + const unsigned int cache_size = cache_end - cache_begin; + std::vector cache_indices(cache_size); + for (unsigned int i = 0; i < cache_size; ++i) + cache_indices[i] = *(cache_begin + sorted_indices_pos[i]); + + IndexSet index_set(cache_indices.back() + 1); + index_set.add_indices(cache_indices.begin(), cache_indices.end()); + index_set.compress(); + LinearAlgebra::ReadWriteVector read_write_vector(index_set); + read_write_vector.import(values, VectorOperation::insert); + + // Copy the elements from read_write_vector and reorder them. + for (unsigned int i = 0; i < cache_size; ++i, ++local_values_begin) + *local_values_begin = read_write_vector[sorted_indices_pos[i]]; + } + + + template static void extract_subvector_to(const LinearAlgebra::EpetraWrappers::Vector &values, diff --git a/include/deal.II/fe/fe_tools_extrapolate.templates.h b/include/deal.II/fe/fe_tools_extrapolate.templates.h index f2acc692d7..850a4c6216 100644 --- a/include/deal.II/fe/fe_tools_extrapolate.templates.h +++ b/include/deal.II/fe/fe_tools_extrapolate.templates.h @@ -1581,6 +1581,21 @@ namespace FETools # ifdef DEAL_II_WITH_MPI + template + void + reinit_distributed(const DoFHandler & dh, + LinearAlgebra::TpetraWrappers::Vector &vector) + { + const parallel::distributed::Triangulation *parallel_tria = + dynamic_cast< + const parallel::distributed::Triangulation *>( + &dh.get_triangulation()); + Assert(parallel_tria != nullptr, ExcNotImplemented()); + + const IndexSet &locally_owned_dofs = dh.locally_owned_dofs(); + vector.reinit(locally_owned_dofs, parallel_tria->get_communicator()); + } + template void reinit_distributed(const DoFHandler & dh, diff --git a/include/deal.II/fe/fe_tools_interpolate.templates.h b/include/deal.II/fe/fe_tools_interpolate.templates.h index 82460f14bb..e2e9b417fe 100644 --- a/include/deal.II/fe/fe_tools_interpolate.templates.h +++ b/include/deal.II/fe/fe_tools_interpolate.templates.h @@ -46,6 +46,7 @@ #include #include #include +#include #include #include @@ -511,6 +512,21 @@ namespace FETools { AssertThrow(false, ExcNotImplemented()); } + + template + void + back_interpolate( + const DoFHandler &, + const AffineConstraints< + typename LinearAlgebra::TpetraWrappers::Vector::value_type> &, + const LinearAlgebra::TpetraWrappers::Vector &, + const DoFHandler &, + const AffineConstraints< + typename LinearAlgebra::TpetraWrappers::Vector::value_type> &, + LinearAlgebra::TpetraWrappers::Vector &) + { + AssertThrow(false, ExcNotImplemented()); + } #endif diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h index 9e67f98178..70d89ade7a 100644 --- a/include/deal.II/lac/read_write_vector.h +++ b/include/deal.II/lac/read_write_vector.h @@ -35,6 +35,7 @@ #ifdef DEAL_II_WITH_TRILINOS # include # include +# include # include @@ -334,6 +335,21 @@ namespace LinearAlgebra std::shared_ptr()); # ifdef DEAL_II_WITH_MPI + /** + * Imports all the elements present in the vector's IndexSet from the input + * vector @p tpetra_vec. 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 + * communication pattern is used multiple times. This can be used to improve + * performance. + */ + void + import(const TpetraWrappers::Vector &tpetra_vec, + VectorOperation::values operation, + const std::shared_ptr + &communication_pattern = + std::shared_ptr()); + /** * Imports all the elements present in the vector's IndexSet from the input * vector @p epetra_vec. VectorOperation::values @p operation is used to @@ -593,6 +609,19 @@ namespace LinearAlgebra protected: #ifdef DEAL_II_WITH_TRILINOS + /** + * Import all the elements present in the vector's IndexSet from the input + * vector @p tpetra_vector. This is an helper function and it should not be + * used directly. + */ + void + import(const Tpetra::Vector<> &tpetra_vector, + const IndexSet & locally_owned_elements, + VectorOperation::values operation, + const MPI_Comm & mpi_comm, + const std::shared_ptr + &communication_pattern); + /** * 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 @@ -627,7 +656,15 @@ namespace LinearAlgebra #if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI) /** - * Return a EpetraWrappers::Communication pattern and store it for future + * Return a TpetraWrappers::CommunicationPattern and store it for future + * use. + */ + TpetraWrappers::CommunicationPattern + create_tpetra_comm_pattern(const IndexSet &source_index_set, + const MPI_Comm &mpi_comm); + + /** + * Return a EpetraWrappers::CommunicationPattern and store it for future * use. */ EpetraWrappers::CommunicationPattern diff --git a/include/deal.II/lac/read_write_vector.templates.h b/include/deal.II/lac/read_write_vector.templates.h index 9a80226481..b4437cfd3f 100644 --- a/include/deal.II/lac/read_write_vector.templates.h +++ b/include/deal.II/lac/read_write_vector.templates.h @@ -34,7 +34,6 @@ #ifdef DEAL_II_WITH_TRILINOS # include -# include # include # include @@ -507,6 +506,93 @@ namespace LinearAlgebra #if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI) + template + void + ReadWriteVector::import( + const Tpetra::Vector<> &vector, + const IndexSet & source_elements, + VectorOperation::values operation, + const MPI_Comm & mpi_comm, + const std::shared_ptr + &communication_pattern) + { + std::shared_ptr + tpetra_comm_pattern; + + // If no communication pattern is given, create one. Otherwise, use the one + // given. + if (communication_pattern == nullptr) + { + // 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.size() == source_stored_elements.size()) && + (source_elements == source_stored_elements)) + { + tpetra_comm_pattern = std::dynamic_pointer_cast< + const TpetraWrappers::CommunicationPattern>(comm_pattern); + if (tpetra_comm_pattern == nullptr) + tpetra_comm_pattern = + std::make_shared( + create_tpetra_comm_pattern(source_elements, mpi_comm)); + } + else + tpetra_comm_pattern = + std::make_shared( + create_tpetra_comm_pattern(source_elements, mpi_comm)); + } + else + { + tpetra_comm_pattern = + std::dynamic_pointer_cast( + communication_pattern); + AssertThrow(tpetra_comm_pattern != nullptr, + ExcMessage( + std::string("The communication pattern is not of type ") + + "LinearAlgebra::TpetraWrappers::CommunicationPattern.")); + } + + Tpetra::Export<> tpetra_export(tpetra_comm_pattern->get_tpetra_export()); + + Tpetra::Vector<> target_vector(tpetra_export.getSourceMap()); + target_vector.doImport(vector, tpetra_export, Tpetra::REPLACE); + + const auto *new_values = target_vector.getData().get(); + const auto size = target_vector.getLocalLength(); + + using size_type = std::decay::type; + + Assert(size == 0 || values != nullptr, ExcInternalError("Export failed.")); + AssertDimension(size, stored_elements.n_elements()); + + if (operation == VectorOperation::insert) + { + for (size_type i = 0; i < size; ++i) + values[i] = new_values[i]; + } + else if (operation == VectorOperation::add) + { + for (size_type i = 0; i < size; ++i) + values[i] += new_values[i]; + } + else if (operation == VectorOperation::min) + { + for (size_type 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) + { + for (size_type 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()); + } + + + template void ReadWriteVector::import( @@ -652,6 +738,23 @@ namespace LinearAlgebra + template + void + ReadWriteVector::import( + const LinearAlgebra::TpetraWrappers::Vector &trilinos_vec, + VectorOperation::values operation, + const std::shared_ptr + &communication_pattern) + { + import(trilinos_vec.trilinos_vector(), + trilinos_vec.locally_owned_elements(), + operation, + trilinos_vec.get_mpi_communicator(), + communication_pattern); + } + + + template void ReadWriteVector::import( @@ -793,6 +896,22 @@ namespace LinearAlgebra #if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI) + template + TpetraWrappers::CommunicationPattern + ReadWriteVector::create_tpetra_comm_pattern( + const IndexSet &source_index_set, + const MPI_Comm &mpi_comm) + { + source_stored_elements = source_index_set; + TpetraWrappers::CommunicationPattern epetra_comm_pattern( + source_stored_elements, stored_elements, mpi_comm); + comm_pattern = std::make_shared( + source_stored_elements, stored_elements, mpi_comm); + + return epetra_comm_pattern; + } + + template EpetraWrappers::CommunicationPattern ReadWriteVector::create_epetra_comm_pattern( diff --git a/include/deal.II/lac/trilinos_sparse_matrix.h b/include/deal.II/lac/trilinos_sparse_matrix.h index 25478951f6..006dc31712 100644 --- a/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_sparse_matrix.h @@ -27,6 +27,7 @@ # include # include # include +# include # include # include # include diff --git a/include/deal.II/lac/trilinos_tpetra_communication_pattern.h b/include/deal.II/lac/trilinos_tpetra_communication_pattern.h new file mode 100644 index 0000000000..8dc7cf7874 --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_communication_pattern.h @@ -0,0 +1,106 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_trilinos_tpetra_communication_pattern_h +#define dealii_trilinos_tpetra_communication_pattern_h + + +#include + +#ifdef DEAL_II_WITH_TRILINOS + +# ifdef DEAL_II_WITH_MPI + +# include + +# include + +# include + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + /** + * This class implements a wrapper to Tpetra::Import and Tpetra::Export. + */ + class CommunicationPattern : public CommunicationPatternBase + { + public: + /** + * Reinitialize the communication pattern. The first argument @p + * vector_space_vector_index_set is the index set associated to a + * VectorSpaceVector object. The second argument @p + * read_write_vector_index_set is the index set associated to a + * ReadWriteVector object. + */ + CommunicationPattern(const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator); + + /** + * Reinitialize the object. + */ + virtual void + reinit(const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator) override; + + /** + * Return the underlying MPI communicator. + */ + virtual const MPI_Comm & + get_mpi_communicator() const override; + + /** + * Return the underlying Tpetra::Import object. + */ + const Tpetra::Import<> & + get_tpetra_import() const; + + /** + * Return the underlying Tpetra::Export object. + */ + const Tpetra::Export<> & + get_tpetra_export() const; + + private: + /** + * Shared pointer to the MPI communicator used. + */ + std::shared_ptr comm; + + /** + * Shared pointer to the Tpetra::Import object used. + */ + std::unique_ptr> tpetra_import; + + /** + * Shared pointer to the Tpetra::Export object used. + */ + std::unique_ptr> tpetra_export; + }; + } // end of namespace TpetraWrappers +} // end of namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +# endif + +#endif + +#endif diff --git a/include/deal.II/lac/trilinos_tpetra_vector.h b/include/deal.II/lac/trilinos_tpetra_vector.h new file mode 100644 index 0000000000..df0b0c99c5 --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_vector.h @@ -0,0 +1,406 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_trilinos_tpetra_vector_h +#define dealii_trilinos_tpetra_vector_h + + +#include + +#if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI) + +# include +# include + +# include +# include +# include +# include + +# include +# include +# include +# include +# include +# include + +# include + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + // Forward declaration + template + class ReadWriteVector; + + namespace TpetraWrappers + { + /** + * This class implements a wrapper to the Trilinos distributed vector + * class Tpetra::Vector. This class is derived from the + * LinearAlgebra::VectorSpaceVector class. Note however that Tpetra only + * works with Number = double. + * + * @ingroup TrilinosWrappers + * @ingroup Vectors + * @author Daniel Arndt, 2018 + */ + class Vector : public VectorSpaceVector, public Subscriptor + { + public: + /** + * Constructor. Create a vector of dimension zero. + */ + Vector(); + + /** + * Copy constructor. Sets the dimension and the partitioning to that of + * the given vector and copies all elements. + */ + Vector(const Vector &V); + + /** + * This constructor takes an IndexSet that defines how to distribute the + * individual components among the MPI processors. Since it also + * includes information about the size of the vector, this is all we + * need to generate a %parallel vector. + */ + explicit Vector(const IndexSet ¶llel_partitioner, + const MPI_Comm &communicator); + + /** + * Reinit functionality. This function destroys the old vector content + * and generates a new one based on the input partitioning. The flag + * omit_zeroing_entries determines whether the vector should be + * filled with zero (false) or left untouched (true). + */ + void + reinit(const IndexSet ¶llel_partitioner, + const MPI_Comm &communicator, + const bool omit_zeroing_entries = false); + + /** + * Change the dimension to that of the vector V. The elements of V are not + * copied. + */ + virtual void + reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries = false) override; + + /** + * Copy function. This function takes a Vector and copies all the + * elements. The Vector will have the same parallel distribution as @p + * V. + */ + Vector & + operator=(const Vector &V); + + /** + * Sets all elements of the vector to the scalar @p s. This operation is + * only allowed if @p s is equal to zero. + */ + virtual Vector & + operator=(const double s) override; + + /** + * Imports 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 + * communication pattern is used multiple times. This can be used to + * improve performance. + */ + virtual void + import( + const ReadWriteVector & V, + VectorOperation::values operation, + std::shared_ptr communication_pattern = + std::shared_ptr()) override; + + /** + * Multiply the entire vector by a fixed factor. + */ + virtual Vector & + operator*=(const double factor) override; + + /** + * Divide the entire vector by a fixed factor. + */ + virtual Vector & + operator/=(const double factor) override; + + /** + * Add the vector @p V to the present one. + */ + virtual Vector & + operator+=(const VectorSpaceVector &V) override; + + /** + * Subtract the vector @p V from the present one. + */ + virtual Vector & + operator-=(const VectorSpaceVector &V) override; + + /** + * Return the scalar product of two vectors. The vectors need to have the + * same layout. + */ + virtual double + operator*(const VectorSpaceVector &V) const override; + + /** + * Add @p a to all components. Note that @p is a scalar not a vector. + */ + virtual void + add(const double a) override; + + /** + * Simple addition of a multiple of a vector, i.e. *this += + * a*V. The vectors need to have the same layout. + */ + virtual void + add(const double a, const VectorSpaceVector &V) override; + + /** + * Multiple addition of multiple of a vector, i.e. *this> += + * a*V+b*W. The vectors need to have the same layout. + */ + virtual void + add(const double a, + const VectorSpaceVector &V, + const double b, + const VectorSpaceVector &W) override; + + /** + * Scaling and simple addition of a multiple of a vector, i.e. *this + * = s*(*this)+a*V. + */ + virtual void + sadd(const double s, + const double a, + const VectorSpaceVector &V) override; + + /** + * Scale each element of this vector by the corresponding element in the + * argument. This function is mostly meant to simulate multiplication + * (and immediate re-assignment) by a diagonal scaling matrix. The + * vectors need to have the same layout. + */ + virtual void + scale(const VectorSpaceVector &scaling_factors) override; + + /** + * Assignment *this = a*V. + */ + virtual void + equ(const double a, const VectorSpaceVector &V) override; + + /** + * Return whether the vector contains only elements with value zero. + */ + virtual bool + all_zero() const override; + + /** + * Return the mean value of the element of this vector. + */ + virtual double + mean_value() const override; + + /** + * Return the l1 norm of the vector (i.e., the sum of the + * absolute values of all entries among all processors). + */ + virtual double + l1_norm() const override; + + /** + * Return the l2 norm of the vector (i.e., the square root of + * the sum of the square of all entries among all processors). + */ + virtual double + l2_norm() const override; + + /** + * Return the maximum norm of the vector (i.e., the maximum absolute value + * among all entries and among all processors). + */ + virtual double + linfty_norm() const override; + + /** + * Performs a combined operation of a vector addition and a subsequent + * inner product, returning the value of the inner product. In other + * words, the result of this function is the same as if the user called + * @code + * this->add(a, V); + * return_value = *this * W; + * @endcode + * + * The reason this function exists is that this operation involves less + * memory transfer than calling the two functions separately. This + * method only needs to load three vectors, @p this, @p V, @p W, whereas + * 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. + * + * For complex-valued vectors, the scalar product in the second step is + * implemented as + * $\left=\sum_i v_i \bar{w_i}$. + */ + virtual double + add_and_dot(const double a, + const VectorSpaceVector &V, + const VectorSpaceVector &W) override; + /** + * This function always returns false and is present only for backward + * compatibility. + */ + bool + has_ghost_elements() const; + + /** + * Return the global size of the vector, equal to the sum of the number of + * locally owned indices among all processors. + */ + virtual size_type + size() const override; + + /** + * Return the MPI communicator object in use with this object. + */ + MPI_Comm + get_mpi_communicator() const; + + /** + * Return an index set that describes which elements of this vector are + * owned by the current processor. As a consequence, the index sets + * returned on different processors if this is a distributed vector will + * form disjoint sets that add up to the complete index set. Obviously, if + * a vector is created on only one processor, then the result would + * satisfy + * @code + * vec.locally_owned_elements() == complete_index_set(vec.size()) + * @endcode + */ + virtual ::dealii::IndexSet + locally_owned_elements() const override; + + /** + * Return a const reference to the underlying Trilinos + * Tpetra::Vector class. + */ + const Tpetra::Vector<> & + trilinos_vector() const; + + /** + * Return a (modifyable) reference to the underlying Trilinos + * Tpetra::Vector class. + */ + Tpetra::Vector<> & + trilinos_vector(); + + /** + * Prints the vector to the output stream @p out. + */ + virtual void + print(std::ostream & out, + const unsigned int precision = 3, + const bool scientific = true, + const bool across = true) const override; + + /** + * Return the memory consumption of this class in bytes. + */ + virtual std::size_t + memory_consumption() const override; + + /** + * The vectors have different partitioning, i.e. they have use different + * IndexSet. + */ + DeclException0(ExcDifferentParallelPartitioning); + + /** + * Attempt to perform an operation between two incompatible vector types. + * + * @ingroup Exceptions + */ + DeclException0(ExcVectorTypeNotCompatible); + + /** + * Exception thrown by an error in Trilinos. + * + * @ingroup Exceptions + */ + DeclException1(ExcTrilinosError, + int, + << "An error with error number " << arg1 + << " occurred while calling a Trilinos function"); + + private: + /** + * Create the CommunicationPattern for the communication between the + * IndexSet @p source_index_set and the current vector based + * on the communicator @p mpi_comm. + */ + void + create_tpetra_comm_pattern(const IndexSet &source_index_set, + const MPI_Comm &mpi_comm); + + /** + * Pointer to the actual Tpetra vector object. + */ + std::unique_ptr> vector; + + /** + * IndexSet of the elements of the last imported vector. + */ + ::dealii::IndexSet source_stored_elements; + + /** + * CommunicationPattern for the communication between the + * source_stored_elements IndexSet and the current vector. + */ + std::shared_ptr + tpetra_comm_pattern; + }; + + + inline bool + Vector::has_ghost_elements() const + { + return false; + } + } // namespace TpetraWrappers +} // namespace LinearAlgebra + + +/** + * Declare dealii::LinearAlgebra::TpetraWrappers::Vector as distributed vector. + */ +template <> +struct is_serial_vector : std::false_type +{}; + +DEAL_II_NAMESPACE_CLOSE + +#endif + +#endif diff --git a/include/deal.II/lac/vector_element_access.h b/include/deal.II/lac/vector_element_access.h index 791f6894c3..b168260588 100644 --- a/include/deal.II/lac/vector_element_access.h +++ b/include/deal.II/lac/vector_element_access.h @@ -18,6 +18,7 @@ #include +#include DEAL_II_NAMESPACE_OPEN @@ -108,7 +109,6 @@ namespace internal vector[0][trilinos_i] = value; } - template <> inline double ElementAccess::get( @@ -122,6 +122,74 @@ namespace internal return vector[0][trilinos_i]; } + + + + template <> + inline void + ElementAccess::add( + const double value, + const types::global_dof_index i, + LinearAlgebra::TpetraWrappers::Vector &V) + { + // Extract local indices in the vector. + Tpetra::Vector<> vector = V.trilinos_vector(); + TrilinosWrappers::types::int_type trilinos_i = + vector.getMap()->getLocalElement( + static_cast(i)); + + vector.sync(); + auto vector_2d = vector.getLocalView(); + auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); + // We're going to modify the data on host. + vector.modify(); + vector_1d(trilinos_i) += value; + vector.sync::device_type::memory_space>(); + } + + + + template <> + inline void + ElementAccess::set( + const double value, + const types::global_dof_index i, + LinearAlgebra::TpetraWrappers::Vector &V) + { + // Extract local indices in the vector. + Tpetra::Vector<> vector = V.trilinos_vector(); + TrilinosWrappers::types::int_type trilinos_i = + vector.getMap()->getLocalElement( + static_cast(i)); + + vector.sync(); + auto vector_2d = vector.getLocalView(); + auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); + // We're going to modify the data on host. + vector.modify(); + vector_1d(trilinos_i) = value; + vector.sync::device_type::memory_space>(); + } + + + template <> + inline double + ElementAccess::get( + const LinearAlgebra::TpetraWrappers::Vector &V, + const types::global_dof_index i) + { + // Extract local indices in the vector. + Tpetra::Vector<> vector = V.trilinos_vector(); + TrilinosWrappers::types::int_type trilinos_i = + vector.getMap()->getLocalElement( + static_cast(i)); + + vector.sync(); + auto vector_2d = vector.getLocalView(); + auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); + // We're going to modify the data on host. + return vector_1d(trilinos_i); + } #endif } // namespace internal diff --git a/include/deal.II/multigrid/mg_transfer.h b/include/deal.II/multigrid/mg_transfer.h index a88bf02884..e4b4f3db2a 100644 --- a/include/deal.II/multigrid/mg_transfer.h +++ b/include/deal.II/multigrid/mg_transfer.h @@ -132,6 +132,38 @@ namespace internal }; # ifdef DEAL_II_WITH_MPI + template <> + struct MatrixSelector + { + using Sparsity = ::dealii::TrilinosWrappers::SparsityPattern; + using Matrix = ::dealii::TrilinosWrappers::SparseMatrix; + + static const bool requires_distributed_sparsity_pattern = false; + + template + static void + reinit(Matrix &matrix, + Sparsity &, + int level, + const SparsityPatternType &sp, + DoFHandlerType & dh) + { + const parallel::Triangulation + *dist_tria = dynamic_cast< + const parallel::Triangulation *>( + &(dh.get_triangulation())); + MPI_Comm communicator = + dist_tria != nullptr ? dist_tria->get_communicator() : MPI_COMM_SELF; + matrix.reinit(dh.locally_owned_mg_dofs(level + 1), + dh.locally_owned_mg_dofs(level), + sp, + communicator, + true); + } + }; + template <> struct MatrixSelector { diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 30e5fa5f1f..aebdafb337 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -70,6 +70,7 @@ #include #include #include +#include #include #include #include diff --git a/source/base/index_set.cc b/source/base/index_set.cc index 3e7f0914cf..6fe2920ffe 100644 --- a/source/base/index_set.cc +++ b/source/base/index_set.cc @@ -520,6 +520,72 @@ IndexSet::fill_index_vector(std::vector &indices) const #ifdef DEAL_II_WITH_TRILINOS +Tpetra::Map<> +IndexSet::make_tpetra_map(const MPI_Comm &communicator, + const bool overlapping) const +{ + compress(); + (void)communicator; + +# ifdef DEBUG + if (!overlapping) + { + const size_type n_global_elements = + Utilities::MPI::sum(n_elements(), communicator); + Assert(n_global_elements == size(), + ExcMessage("You are trying to create an Tpetra::Map object " + "that partitions elements of an index set " + "between processors. However, the union of the " + "index sets on different processors does not " + "contain all indices exactly once: the sum of " + "the number of entries the various processors " + "want to store locally is " + + Utilities::to_string(n_global_elements) + + " whereas the total size of the object to be " + "allocated is " + + Utilities::to_string(size()) + + ". In other words, there are " + "either indices that are not spoken for " + "by any processor, or there are indices that are " + "claimed by multiple processors.")); + } +# endif + + // Find out if the IndexSet is ascending and 1:1. This corresponds to a + // linear Tpetra::Map. Overlapping IndexSets are never 1:1. + const bool linear = + overlapping ? false : is_ascending_and_one_to_one(communicator); + if (linear) + return Tpetra::Map<>(size(), + n_elements(), + 0, +# ifdef DEAL_II_WITH_MPI + Teuchos::rcp(new Teuchos::MpiComm(communicator)) +# else + Teuchos::rcp(new Teuchos::Comm()) +# endif + ); + else + { + std::vector indices; + fill_index_vector(indices); + std::vector int_indices(indices.size()); + std::copy(indices.begin(), indices.end(), int_indices.begin()); + return Tpetra::Map<>(size(), + (n_elements() > 0 ? int_indices.data() : nullptr), + n_elements(), + 0, +# ifdef DEAL_II_WITH_MPI + Teuchos::rcp(new Teuchos::MpiComm(communicator)) +# else + Teuchos::rcp(new Teuchos::Comm()) +# endif + ); + } +} + + + Epetra_Map IndexSet::make_trilinos_map(const MPI_Comm &communicator, const bool overlapping) const diff --git a/source/base/time_stepping.cc b/source/base/time_stepping.cc index 9bb02320eb..3e4d6bed67 100644 --- a/source/base/time_stepping.cc +++ b/source/base/time_stepping.cc @@ -22,6 +22,7 @@ #include #include #include +#include #include #include diff --git a/source/base/utilities.cc b/source/base/utilities.cc index 27b798cc12..ff4d0d9573 100644 --- a/source/base/utilities.cc +++ b/source/base/utilities.cc @@ -979,6 +979,20 @@ namespace Utilities + const Teuchos::RCP> & + tpetra_comm_self() + { +# ifdef DEAL_II_WITH_MPI + static auto communicator = Teuchos::RCP>( + new Teuchos::MpiComm(MPI_COMM_SELF)); +# else + static auto communicator = + Teuchos::RCP>(new Teuchos::Comm()); +# endif + + return communicator; + } + const Epetra_Comm & comm_self() { diff --git a/source/dofs/dof_accessor_get.cc b/source/dofs/dof_accessor_get.cc index 4f1bdc2b03..53b47d2238 100644 --- a/source/dofs/dof_accessor_get.cc +++ b/source/dofs/dof_accessor_get.cc @@ -32,6 +32,7 @@ #include #include #include +#include #include #include diff --git a/source/dofs/dof_accessor_set.cc b/source/dofs/dof_accessor_set.cc index 556dcd48c8..a116ab0858 100644 --- a/source/dofs/dof_accessor_set.cc +++ b/source/dofs/dof_accessor_set.cc @@ -32,6 +32,7 @@ #include #include #include +#include #include #include diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index 920f57a61a..b2c44908a6 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -39,6 +39,7 @@ #include #include #include +#include #include #include #include diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc index 5289c9d1e7..1659ccfb3a 100644 --- a/source/fe/mapping_fe_field.cc +++ b/source/fe/mapping_fe_field.cc @@ -44,6 +44,7 @@ #include #include #include +#include #include #include diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index a6ce60360f..973837219e 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -132,6 +132,8 @@ IF(DEAL_II_WITH_TRILINOS) trilinos_solver.cc trilinos_sparse_matrix.cc trilinos_sparsity_pattern.cc + trilinos_tpetra_communication_pattern.cc + trilinos_tpetra_vector.cc trilinos_vector.cc ) diff --git a/source/lac/solver.cc b/source/lac/solver.cc index 6d432ae184..7bf9e30c3e 100644 --- a/source/lac/solver.cc +++ b/source/lac/solver.cc @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include diff --git a/source/lac/trilinos_sparse_matrix.cc b/source/lac/trilinos_sparse_matrix.cc index 291929c5c3..350006da3c 100644 --- a/source/lac/trilinos_sparse_matrix.cc +++ b/source/lac/trilinos_sparse_matrix.cc @@ -103,6 +103,36 @@ namespace TrilinosWrappers { return V.trilinos_vector()[0] + V.trilinos_vector().MyLength(); } + + template <> + double * + begin(LinearAlgebra::TpetraWrappers::Vector &V) + { + return V.trilinos_vector().getDataNonConst().get(); + } + + template <> + const double * + begin(const LinearAlgebra::TpetraWrappers::Vector &V) + { + return V.trilinos_vector().getData().get(); + } + + template <> + double * + end(LinearAlgebra::TpetraWrappers::Vector &V) + { + return V.trilinos_vector().getDataNonConst().get() + + V.trilinos_vector().getLocalLength(); + } + + template <> + const double * + end(const LinearAlgebra::TpetraWrappers::Vector &V) + { + return V.trilinos_vector().getData().get() + + V.trilinos_vector().getLocalLength(); + } # endif } // namespace internal @@ -3454,6 +3484,11 @@ namespace TrilinosWrappers dealii::LinearAlgebra::distributed::Vector &, const dealii::LinearAlgebra::distributed::Vector &) const; # ifdef DEAL_II_WITH_MPI + template void + SparseMatrix::vmult( + dealii::LinearAlgebra::TpetraWrappers::Vector &, + const dealii::LinearAlgebra::TpetraWrappers::Vector &) const; + template void SparseMatrix::vmult( dealii::LinearAlgebra::EpetraWrappers::Vector &, @@ -3469,6 +3504,11 @@ namespace TrilinosWrappers dealii::LinearAlgebra::distributed::Vector &, const dealii::LinearAlgebra::distributed::Vector &) const; # ifdef DEAL_II_WITH_MPI + template void + SparseMatrix::Tvmult( + dealii::LinearAlgebra::TpetraWrappers::Vector &, + const dealii::LinearAlgebra::TpetraWrappers::Vector &) const; + template void SparseMatrix::Tvmult( dealii::LinearAlgebra::EpetraWrappers::Vector &, @@ -3484,6 +3524,11 @@ namespace TrilinosWrappers dealii::LinearAlgebra::distributed::Vector &, const dealii::LinearAlgebra::distributed::Vector &) const; # ifdef DEAL_II_WITH_MPI + template void + SparseMatrix::vmult_add( + dealii::LinearAlgebra::TpetraWrappers::Vector &, + const dealii::LinearAlgebra::TpetraWrappers::Vector &) const; + template void SparseMatrix::vmult_add( dealii::LinearAlgebra::EpetraWrappers::Vector &, @@ -3499,6 +3544,11 @@ namespace TrilinosWrappers dealii::LinearAlgebra::distributed::Vector &, const dealii::LinearAlgebra::distributed::Vector &) const; # ifdef DEAL_II_WITH_MPI + template void + SparseMatrix::Tvmult_add( + dealii::LinearAlgebra::TpetraWrappers::Vector &, + const dealii::LinearAlgebra::TpetraWrappers::Vector &) const; + template void SparseMatrix::Tvmult_add( dealii::LinearAlgebra::EpetraWrappers::Vector &, diff --git a/source/lac/trilinos_tpetra_communication_pattern.cc b/source/lac/trilinos_tpetra_communication_pattern.cc new file mode 100644 index 0000000000..e78dfda3d3 --- /dev/null +++ b/source/lac/trilinos_tpetra_communication_pattern.cc @@ -0,0 +1,104 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include + +#include + +#ifdef DEAL_II_WITH_TRILINOS + +# ifdef DEAL_II_WITH_MPI + +# include + +# include + +# include + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + CommunicationPattern::CommunicationPattern( + const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator) + { + // virtual functions called in constructors and destructors never use the + // override in a derived class + // for clarity be explicit on which function is called + CommunicationPattern::reinit(vector_space_vector_index_set, + read_write_vector_index_set, + communicator); + } + + + + void + CommunicationPattern::reinit(const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator) + { + comm = std::make_shared(communicator); + + auto vector_space_vector_map = Teuchos::rcp(new Tpetra::Map<>( + vector_space_vector_index_set.make_tpetra_map(*comm, false))); + auto read_write_vector_map = Teuchos::rcp(new Tpetra::Map<>( + read_write_vector_index_set.make_tpetra_map(*comm, true))); + + // Target map is read_write_vector_map + // Source map is vector_space_vector_map. This map must have uniquely + // owned GID. + tpetra_import = + std_cxx14::make_unique>(read_write_vector_map, + vector_space_vector_map); + tpetra_export = + std_cxx14::make_unique>(read_write_vector_map, + vector_space_vector_map); + } + + + + const MPI_Comm & + CommunicationPattern::get_mpi_communicator() const + { + return *comm; + } + + + + const Tpetra::Import<> & + CommunicationPattern::get_tpetra_import() const + { + return *tpetra_import; + } + + + + const Tpetra::Export<> & + CommunicationPattern::get_tpetra_export() const + { + return *tpetra_export; + } + } // namespace TpetraWrappers +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +# endif + +#endif diff --git a/source/lac/trilinos_tpetra_vector.cc b/source/lac/trilinos_tpetra_vector.cc new file mode 100644 index 0000000000..0a52aa1972 --- /dev/null +++ b/source/lac/trilinos_tpetra_vector.cc @@ -0,0 +1,632 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#include + +#include + +#ifdef DEAL_II_WITH_TRILINOS + +# ifdef DEAL_II_WITH_MPI + +# include + +# include + +# include + +# include +# include +# include + +# include + + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + Vector::Vector() + : vector(new Tpetra::Vector<>(Teuchos::RCP>( + new Tpetra::Map<>(0, 0, Utilities::Trilinos::tpetra_comm_self())))) + {} + + + + Vector::Vector(const Vector &V) + : Subscriptor() + , vector(new Tpetra::Vector<>(V.trilinos_vector(), Teuchos::Copy)) + {} + + + + Vector::Vector(const IndexSet ¶llel_partitioner, + const MPI_Comm &communicator) + : vector(new Tpetra::Vector<>(Teuchos::rcp(new Tpetra::Map<>( + parallel_partitioner.make_tpetra_map(communicator, false))))) + {} + + + + void + Vector::reinit(const IndexSet ¶llel_partitioner, + const MPI_Comm &communicator, + const bool omit_zeroing_entries) + { + Tpetra::Map<> input_map = + parallel_partitioner.make_tpetra_map(communicator, false); + if (vector->getMap()->isSameAs(input_map) == false) + vector = std_cxx14::make_unique>( + Teuchos::rcp(new Tpetra::Map<>(input_map))); + else if (omit_zeroing_entries == false) + { + vector->putScalar(0.); + } + } + + + + void + Vector::reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries) + { + // Check that casting will work. + Assert(dynamic_cast(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast(V); + + reinit(down_V.locally_owned_elements(), + down_V.get_mpi_communicator(), + omit_zeroing_entries); + } + + + + Vector & + Vector::operator=(const Vector &V) + { + // Distinguish three cases: + // - First case: both vectors have the same layout. + // - Second case: both vectors have the same size but different layout. + // - Third case: the vectors have different size. + if (vector->getMap()->isSameAs(*(V.trilinos_vector().getMap()))) + *vector = V.trilinos_vector(); + else + { + if (size() == V.size()) + { + Tpetra::Import<> data_exchange(vector->getMap(), + V.trilinos_vector().getMap()); + + vector->doImport(V.trilinos_vector(), + data_exchange, + Tpetra::REPLACE); + } + else + vector = + std_cxx14::make_unique>(V.trilinos_vector()); + } + + return *this; + } + + + + Vector & + Vector::operator=(const double s) + { + Assert(s == 0., ExcMessage("Only 0 can be assigned to a vector.")); + + vector->putScalar(s); + + return *this; + } + + + + void + Vector::import( + const ReadWriteVector & V, + VectorOperation::values operation, + std::shared_ptr communication_pattern) + { + // If no communication pattern is given, create one. Otherwsie, use the + // one given. + if (communication_pattern == nullptr) + { + // 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.size() != + V.get_stored_elements().size()) || + (source_stored_elements != V.get_stored_elements())) + { + const Teuchos::MpiComm *mpi_comm = + dynamic_cast *>( + vector->getMap()->getComm().get()); + Assert(mpi_comm != nullptr, ExcInternalError()); + create_tpetra_comm_pattern(V.get_stored_elements(), + *(mpi_comm->getRawMpiComm())()); + } + } + else + { + tpetra_comm_pattern = std::dynamic_pointer_cast< + const TpetraWrappers::CommunicationPattern>(communication_pattern); + AssertThrow( + tpetra_comm_pattern != nullptr, + ExcMessage( + std::string("The communication pattern is not of type ") + + "LinearAlgebra::TpetraWrappers::CommunicationPattern.")); + } + + Tpetra::Export<> tpetra_export(tpetra_comm_pattern->get_tpetra_export()); + Tpetra::Vector<> source_vector(tpetra_export.getSourceMap()); + + source_vector.sync(); + auto x_2d = source_vector.getLocalView(); + auto x_1d = Kokkos::subview(x_2d, Kokkos::ALL(), 0); + source_vector.modify(); + const size_t localLength = source_vector.getLocalLength(); + auto values_it = V.begin(); + for (size_t k = 0; k < localLength; ++k) + x_1d(k) = *values_it++; + source_vector.sync::device_type::memory_space>(); + if (operation == VectorOperation::insert) + vector->doExport(source_vector, tpetra_export, Tpetra::REPLACE); + else if (operation == VectorOperation::add) + vector->doExport(source_vector, tpetra_export, Tpetra::ADD); + else + AssertThrow(false, ExcNotImplemented()); + } + + + + Vector & + Vector::operator*=(const double factor) + { + AssertIsFinite(factor); + vector->scale(factor); + + return *this; + } + + + + Vector & + Vector::operator/=(const double factor) + { + AssertIsFinite(factor); + Assert(factor != 0., ExcZero()); + *this *= 1. / factor; + + return *this; + } + + + + Vector & + Vector::operator+=(const VectorSpaceVector &V) + { + // Check that casting will work. + Assert(dynamic_cast(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast(V); + // If the maps are the same we can Update right away. + if (vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap()))) + { + vector->update(1., down_V.trilinos_vector(), 1.); + } + else + { + Assert(this->size() == down_V.size(), + ExcDimensionMismatch(this->size(), down_V.size())); + + // TODO: The code doesn't work as expected so we use a workaround. + /*Tpetra::Export<> data_exchange(vector->getMap(), + down_V.trilinos_vector().getMap()); + vector->doExport(down_V.trilinos_vector(), + data_exchange, + Tpetra::ADD);*/ + + Tpetra::Vector<> dummy(vector->getMap(), false); + Tpetra::Import<> data_exchange(dummy.getMap(), + down_V.trilinos_vector().getMap()); + + dummy.doExport(down_V.trilinos_vector(), + data_exchange, + Tpetra::REPLACE); + + vector->update(1.0, dummy, 1.0); + } + + return *this; + } + + + + Vector & + Vector::operator-=(const VectorSpaceVector &V) + { + this->add(-1., V); + + return *this; + } + + + + double Vector::operator*(const VectorSpaceVector &V) const + { + // Check that casting will work. + Assert(dynamic_cast(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast(V); + Assert(this->size() == down_V.size(), + ExcDimensionMismatch(this->size(), down_V.size())); + Assert(vector->getMap()->isSameAs(*down_V.trilinos_vector().getMap()), + ExcDifferentParallelPartitioning()); + + return vector->dot(down_V.trilinos_vector()); + } + + + + void + Vector::add(const double a) + { + AssertIsFinite(a); + + vector->sync(); + auto vector_2d = vector->getLocalView(); + auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); + vector->modify(); + const size_t localLength = vector->getLocalLength(); + for (size_t k = 0; k < localLength; ++k) + { + vector_1d(k) += a; + } + vector->sync::device_type::memory_space>(); + } + + + + void + Vector::add(const double a, const VectorSpaceVector &V) + { + // Check that casting will work. + Assert(dynamic_cast(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast(V); + AssertIsFinite(a); + Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())), + ExcDifferentParallelPartitioning()); + + vector->update(a, down_V.trilinos_vector(), 1.); + } + + + + void + Vector::add(const double a, + const VectorSpaceVector &V, + const double b, + const VectorSpaceVector &W) + { + // Check that casting will work. + Assert(dynamic_cast(&V) != nullptr, + ExcVectorTypeNotCompatible()); + // Check that casting will work. + Assert(dynamic_cast(&W) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast(V); + // Downcast W. If fails, throws an exception. + const Vector &down_W = dynamic_cast(W); + Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())), + ExcDifferentParallelPartitioning()); + Assert(vector->getMap()->isSameAs(*(down_W.trilinos_vector().getMap())), + ExcDifferentParallelPartitioning()); + AssertIsFinite(a); + AssertIsFinite(b); + + vector->update( + a, down_V.trilinos_vector(), b, down_W.trilinos_vector(), 1.); + } + + + + void + Vector::sadd(const double s, + const double a, + const VectorSpaceVector &V) + { + // Check that casting will work. + Assert(dynamic_cast(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + *this *= s; + // Downcast V. It fails, throws an exception. + const Vector &down_V = dynamic_cast(V); + Vector tmp(down_V); + tmp *= a; + *this += tmp; + } + + + + void + Vector::scale(const VectorSpaceVector &scaling_factors) + { + // Check that casting will work. + Assert(dynamic_cast(&scaling_factors) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast scaling_factors. If fails, throws an exception. + const Vector &down_scaling_factors = + dynamic_cast(scaling_factors); + Assert(vector->getMap()->isSameAs( + *(down_scaling_factors.trilinos_vector().getMap())), + ExcDifferentParallelPartitioning()); + + vector->elementWiseMultiply(1., + *down_scaling_factors.vector, + *vector, + 0.); + } + + + + void + Vector::equ(const double a, const VectorSpaceVector &V) + { + // Check that casting will work. + Assert(dynamic_cast(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast(V); + // If we don't have the same map, copy. + if (vector->getMap()->isSameAs(*down_V.trilinos_vector().getMap()) == + false) + this->sadd(0., a, V); + else + { + // Otherwise, just update + vector->update(a, down_V.trilinos_vector(), 0.); + } + } + + + + bool + Vector::all_zero() const + { + // get a representation of the vector and + // loop over all the elements + double * start_ptr = vector->getDataNonConst().get(); + const double *ptr = start_ptr, + *eptr = start_ptr + vector->getLocalLength(); + unsigned int flag = 0; + while (ptr != eptr) + { + if (*ptr != 0) + { + flag = 1; + break; + } + ++ptr; + } + + // Check that the vector is zero on _all_ processors. + const Teuchos::MpiComm *mpi_comm = + dynamic_cast *>( + vector->getMap()->getComm().get()); + Assert(mpi_comm != nullptr, ExcInternalError()); + unsigned int num_nonzero = + Utilities::MPI::sum(flag, *(mpi_comm->getRawMpiComm())()); + + return num_nonzero == 0; + } + + + + double + Vector::mean_value() const + { + return vector->meanValue(); + } + + + + double + Vector::l1_norm() const + { + return vector->norm1(); + } + + + + double + Vector::l2_norm() const + { + return vector->norm2(); + } + + + + double + Vector::linfty_norm() const + { + return vector->normInf(); + } + + + + double + Vector::add_and_dot(const double a, + const VectorSpaceVector &V, + const VectorSpaceVector &W) + { + this->add(a, V); + + return *this * W; + } + + + + Vector::size_type + Vector::size() const + { + return vector->getGlobalLength(); + } + + + + MPI_Comm + Vector::get_mpi_communicator() const + { + const auto tpetra_comm = dynamic_cast *>( + vector->getMap()->getComm().get()); + Assert(tpetra_comm != nullptr, ExcInternalError()); + return *(tpetra_comm->getRawMpiComm())(); + } + + + + ::dealii::IndexSet + Vector::locally_owned_elements() const + { + IndexSet is(size()); + + // easy case: local range is contiguous + if (vector->getMap()->isContiguous()) + { +# ifndef DEAL_II_WITH_64BIT_INDICES + is.add_range(vector->getMap()->getMinGlobalIndex(), + vector->getMap()->getMaxGlobalIndex() + 1); +# else + is.add_range(vector->getMap()->getMinGlobalIndex(), + vector->getMap()->getMaxGlobalIndex() + 1); +# endif + } + else if (vector->getLocalLength() > 0) + { + const size_type n_indices = vector->getLocalLength(); + auto vector_indices = vector->getMap()->getMyGlobalIndices(); + is.add_indices((unsigned int *)&vector_indices[0], + (unsigned int *)&vector_indices[0] + n_indices); + } + is.compress(); + + return is; + } + + + + const Tpetra::Vector<> & + Vector::trilinos_vector() const + { + return *vector; + } + + + + Tpetra::Vector<> & + Vector::trilinos_vector() + { + return *vector; + } + + + + void + Vector::print(std::ostream & out, + const unsigned int precision, + const bool scientific, + const bool across) const + { + AssertThrow(out, ExcIO()); + boost::io::ios_flags_saver restore_flags(out); + + // Get a representation of the vector and loop over all + // the elements + const auto val = vector->get1dView(); + + out.precision(precision); + if (scientific) + out.setf(std::ios::scientific, std::ios::floatfield); + else + out.setf(std::ios::fixed, std::ios::floatfield); + + vector->sync(); + auto vector_2d = vector->getLocalView(); + auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); + const size_t local_length = vector->getLocalLength(); + + if (across) + for (unsigned int i = 0; i < local_length; ++i) + out << i << "->" << vector->getMap()->getGlobalElement(i) << ": " + << vector_1d(i) << std::endl; + // out << val[i] << ' '; + else + for (unsigned int i = 0; i < local_length; ++i) + out << val[i] << std::endl; + out << std::endl; + + // restore the representation + // of the vector + AssertThrow(out, ExcIO()); + } + + + + std::size_t + Vector::memory_consumption() const + { + return sizeof(*this) + + vector->getLocalLength() * + (sizeof(double) + sizeof(TrilinosWrappers::types::int_type)); + } + + + + void + Vector::create_tpetra_comm_pattern(const IndexSet &source_index_set, + const MPI_Comm &mpi_comm) + { + source_stored_elements = source_index_set; + tpetra_comm_pattern = + std::make_shared( + locally_owned_elements(), source_index_set, mpi_comm); + } + } // namespace TpetraWrappers +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +# endif + +#endif diff --git a/source/lac/vector_memory.cc b/source/lac/vector_memory.cc index f15f32800e..3fa5affd9b 100644 --- a/source/lac/vector_memory.cc +++ b/source/lac/vector_memory.cc @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include diff --git a/source/multigrid/mg_base.cc b/source/multigrid/mg_base.cc index 8744c6cf40..a32b68a096 100644 --- a/source/multigrid/mg_base.cc +++ b/source/multigrid/mg_base.cc @@ -19,6 +19,7 @@ #include #include #include +#include #include #include diff --git a/source/numerics/data_out_dof_data_codim.cc b/source/numerics/data_out_dof_data_codim.cc index e085c3d382..676a33714e 100644 --- a/source/numerics/data_out_dof_data_codim.cc +++ b/source/numerics/data_out_dof_data_codim.cc @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include diff --git a/source/numerics/derivative_approximation.cc b/source/numerics/derivative_approximation.cc index 73f574bdc1..d062b4ec2b 100644 --- a/source/numerics/derivative_approximation.cc +++ b/source/numerics/derivative_approximation.cc @@ -39,6 +39,7 @@ #include #include #include +#include #include #include diff --git a/tests/trilinos/tpetra_vector_01.cc b/tests/trilinos/tpetra_vector_01.cc new file mode 100644 index 0000000000..c10ac87bb1 --- /dev/null +++ b/tests/trilinos/tpetra_vector_01.cc @@ -0,0 +1,225 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +// Check LinearAlgebra::TpetraWrappers::Vector assignment and import + + +void +test() +{ + IndexSet parallel_partitioner_1(10); + IndexSet parallel_partitioner_2(10); + unsigned int rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + if (rank == 0) + { + parallel_partitioner_1.add_range(0, 5); + parallel_partitioner_2.add_range(0, 3); + } + else + { + parallel_partitioner_1.add_range(5, 10); + parallel_partitioner_2.add_range(3, 10); + } + parallel_partitioner_1.compress(); + parallel_partitioner_2.compress(); + LinearAlgebra::TpetraWrappers::Vector a; + LinearAlgebra::TpetraWrappers::Vector b(parallel_partitioner_1, + MPI_COMM_WORLD); + LinearAlgebra::TpetraWrappers::Vector c(b); + + AssertThrow(a.size() == 0, ExcMessage("Vector has the wrong size.")); + AssertThrow(b.size() == 10, ExcMessage("Vector has the wrong size.")); + AssertThrow(c.size() == 10, ExcMessage("Vector has the wrong size.")); + + a.reinit(parallel_partitioner_2, MPI_COMM_WORLD); + AssertThrow(a.size() == 10, ExcMessage("Vector has the wrong size.")); + + AssertThrow(parallel_partitioner_1 == b.locally_owned_elements(), + ExcMessage("IndexSet has been modified.")); + AssertThrow(parallel_partitioner_2 == a.locally_owned_elements(), + ExcMessage("IndexSet has been modified.")); + + IndexSet read_write_index_set(10); + if (rank == 0) + read_write_index_set.add_range(0, 6); + else + read_write_index_set.add_range(4, 10); + read_write_index_set.compress(); + + LinearAlgebra::ReadWriteVector read_write_1(read_write_index_set); + LinearAlgebra::ReadWriteVector read_write_2(read_write_index_set); + LinearAlgebra::ReadWriteVector read_write_3(read_write_index_set); + if (rank == 0) + { + for (unsigned int i = 0; i < 6; ++i) + { + read_write_1[i] = i; + read_write_2[i] = 5. + i; + } + } + else + { + for (unsigned int i = 4; i < 10; ++i) + { + read_write_1[i] = i; + read_write_2[i] = 5. + i; + } + } + + a.import(read_write_2, VectorOperation::insert); + AssertThrow(a.size() == 10, ExcMessage("Vector has the wrong size.")); + + read_write_3.import(a, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 6; ++i) + { + AssertThrow(read_write_2[i] == read_write_3[i], + ExcMessage("Vector a has been modified.")); + } + } + else + { + for (unsigned int i = 4; i < 10; ++i) + AssertThrow(read_write_2[i] == read_write_3[i], + ExcMessage("Vector a has been modified.")); + } + + b.import(read_write_1, VectorOperation::insert); + AssertThrow(b.size() == 10, ExcMessage("Vector has the wrong size.")); + + read_write_3.import(b, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 6; ++i) + AssertThrow(read_write_1[i] == read_write_3[i], + ExcMessage("Vector b has been modified.")); + } + else + { + for (unsigned int i = 4; i < 10; ++i) + AssertThrow(read_write_1[i] == read_write_3[i], + ExcMessage("Vector b has been modified.")); + } + + c.import(read_write_2, VectorOperation::insert); + AssertThrow(c.size() == 10, ExcMessage("Vector has the wrong size.")); + + read_write_3.import(c, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 6; ++i) + AssertThrow(read_write_2[i] == read_write_3[i], + ExcMessage("Vector c has been modified.")); + } + else + { + for (unsigned int i = 4; i < 10; ++i) + AssertThrow(read_write_2[i] == read_write_3[i], + ExcMessage("Vector c has been modified.")); + } + + a *= 2; + read_write_3.import(a, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 6; ++i) + AssertThrow(2. * read_write_2[i] == read_write_3[i], + ExcMessage("Problem in operator *=.")); + } + else + { + for (unsigned int i = 4; i < 10; ++i) + AssertThrow(2. * read_write_2[i] == read_write_3[i], + ExcMessage("Problem in operator *=.")); + } + + c /= 2.; + read_write_3.import(c, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 6; ++i) + AssertThrow(0.5 * read_write_2[i] == read_write_3[i], + ExcMessage("Problem in operator /=.")); + } + else + { + for (unsigned int i = 4; i < 10; ++i) + AssertThrow(0.5 * read_write_2[i] == read_write_3[i], + ExcMessage("Problem in operator /=.")); + } + + b += a; + read_write_3.import(b, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 6; ++i) + AssertThrow(2. * read_write_2[i] + read_write_1[i] == read_write_3[i], + ExcMessage("Problem in operator +=.")); + } + else + { + for (unsigned int i = 4; i < 10; ++i) + AssertThrow(2. * read_write_2[i] + read_write_1[i] == read_write_3[i], + ExcMessage("Problem in operator +=.")); + } + + b -= c; + read_write_3.import(b, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 6; ++i) + AssertThrow(1.5 * read_write_2[i] + read_write_1[i] == read_write_3[i], + ExcMessage("Problem in operator -=.")); + } + else + { + for (unsigned int i = 4; i < 10; ++i) + AssertThrow(1.5 * read_write_2[i] + read_write_1[i] == read_write_3[i], + ExcMessage("Problem in operator -=.")); + } + + 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 *.")); +} + + +int +main(int argc, char **argv) +{ + initlog(); + deallog.depth_console(0); + + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + test(); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/trilinos/tpetra_vector_01.mpirun=2.output b/tests/trilinos/tpetra_vector_01.mpirun=2.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos/tpetra_vector_01.mpirun=2.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/trilinos/tpetra_vector_02.cc b/tests/trilinos/tpetra_vector_02.cc new file mode 100644 index 0000000000..8e48753530 --- /dev/null +++ b/tests/trilinos/tpetra_vector_02.cc @@ -0,0 +1,216 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 - 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +// Check LinearAlgebra::TpetraWrappers::Vector add and sadd. + +void +test() +{ + IndexSet parallel_partitioner_1(10); + IndexSet parallel_partitioner_2(10); + unsigned int rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + if (rank == 0) + { + parallel_partitioner_1.add_range(0, 5); + parallel_partitioner_2.add_range(0, 3); + } + else + { + parallel_partitioner_1.add_range(5, 10); + parallel_partitioner_2.add_range(3, 10); + } + parallel_partitioner_1.compress(); + parallel_partitioner_2.compress(); + LinearAlgebra::TpetraWrappers::Vector a(parallel_partitioner_1, + MPI_COMM_WORLD); + LinearAlgebra::TpetraWrappers::Vector b(parallel_partitioner_1, + MPI_COMM_WORLD); + LinearAlgebra::TpetraWrappers::Vector c(parallel_partitioner_2, + MPI_COMM_WORLD); + + IndexSet read_write_index_set(10); + if (rank == 0) + read_write_index_set.add_range(0, 5); + else + read_write_index_set.add_range(5, 10); + read_write_index_set.compress(); + + LinearAlgebra::ReadWriteVector read_write_1(read_write_index_set); + LinearAlgebra::ReadWriteVector read_write_2(read_write_index_set); + LinearAlgebra::ReadWriteVector read_write_3(read_write_index_set); + if (rank == 0) + { + for (unsigned int i = 0; i < 5; ++i) + { + read_write_1[i] = i; + read_write_2[i] = 5. + i; + } + } + else + { + for (unsigned int i = 5; i < 10; ++i) + { + read_write_1[i] = i; + read_write_2[i] = 5. + i; + } + } + + a.import(read_write_1, VectorOperation::insert); + b.import(read_write_2, VectorOperation::insert); + c.import(read_write_2, VectorOperation::insert); + + a.add(1.); + read_write_3.import(a, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 5; ++i) + AssertThrow(1. + read_write_1[i] == read_write_3[i], + ExcMessage("Problem in add(scalar).")); + } + else + { + for (unsigned int i = 5; i < 10; ++i) + AssertThrow(1. + read_write_1[i] == read_write_3[i], + ExcMessage("Problem in add(scalar).")); + } + + a.add(2., b); + read_write_3.import(a, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 5; ++i) + AssertThrow(1. + read_write_1[i] + 2. * read_write_2[i] == + read_write_3[i], + ExcMessage("Problem in add(scalar,Vector).")); + } + else + { + for (unsigned int i = 5; i < 10; ++i) + AssertThrow(1. + read_write_1[i] + 2. * read_write_2[i] == + read_write_3[i], + ExcMessage("Problem in add(scalar,Vector).")); + } + + + LinearAlgebra::TpetraWrappers::Vector d(a); + a.add(2., b, 3., d); + read_write_3.import(a, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 5; ++i) + AssertThrow(4. + 4. * read_write_1[i] + 10. * read_write_2[i] == + read_write_3[i], + ExcMessage("Problem in add(scalar,Vector,scalar,Vector).")); + } + else + { + for (unsigned int i = 5; i < 10; ++i) + AssertThrow(4. + 4. * read_write_1[i] + 10. * read_write_2[i] == + read_write_3[i], + ExcMessage("Problem in add(scalar,Vector,scalar,Vector).")); + } + + + a.import(read_write_1, VectorOperation::insert); + a.sadd(3., 2., c); + read_write_3.import(a, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 5; ++i) + AssertThrow(3. * read_write_1[i] + 2. * read_write_2[i] == + read_write_3[i], + ExcMessage("Problem in sadd(scalar,scalar,Vector).")); + } + else + { + for (unsigned int i = 5; i < 10; ++i) + AssertThrow(3. * read_write_1[i] + 2. * read_write_2[i] == + read_write_3[i], + ExcMessage("Problem in sadd(scalar,scalar,Vector).")); + } + + + a.import(read_write_1, VectorOperation::insert); + a.scale(b); + read_write_3.import(a, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 5; ++i) + AssertThrow(read_write_1[i] * read_write_2[i] == read_write_3[i], + ExcMessage("Problem in scale.")); + } + else + { + for (unsigned int i = 5; i < 10; ++i) + AssertThrow(read_write_1[i] * read_write_2[i] == read_write_3[i], + ExcMessage("Problem in scale.")); + } + + + a.equ(2., c); + read_write_3.import(a, VectorOperation::insert); + if (rank == 0) + { + for (unsigned int i = 0; i < 5; ++i) + AssertThrow(2. * read_write_2[i] == read_write_3[i], + ExcMessage("Problem in scale.")); + } + else + { + for (unsigned int i = 5; i < 10; ++i) + AssertThrow(2. * read_write_2[i] == read_write_3[i], + ExcMessage("Problem in equ.")); + } + + + AssertThrow(b.l1_norm() == 95., ExcMessage("Problem in l1_norm.")); + + const double eps = 1e-6; + AssertThrow(std::fabs(b.l2_norm() - 31.3847096) < eps, + ExcMessage("Problem in l2_norm")); + + AssertThrow(b.linfty_norm() == 14., ExcMessage("Problem in linfty_norm.")); + + a.import(read_write_1, VectorOperation::insert); + const double val = a.add_and_dot(2., a, b); + AssertThrow(val == 1530., ExcMessage("Problem in add_and_dot")); +} + + +int +main(int argc, char **argv) +{ + initlog(); + deallog.depth_console(0); + + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + test(); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/trilinos/tpetra_vector_02.mpirun=2.output b/tests/trilinos/tpetra_vector_02.mpirun=2.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/trilinos/tpetra_vector_02.mpirun=2.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/trilinos/tpetra_vector_03.cc b/tests/trilinos/tpetra_vector_03.cc new file mode 100644 index 0000000000..4927aaedb8 --- /dev/null +++ b/tests/trilinos/tpetra_vector_03.cc @@ -0,0 +1,101 @@ +// --------------------------------------------------------------------- +// +// 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Check LinearAlgebra::TpetraWrappers::Vector::import +// for VectorOperation::add and check LinearAlgebra::ReadWriteVector +// for VectorOperation::add/min/max. + +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +void +test() +{ + unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + + IndexSet locally_owned(n_procs * 2); + locally_owned.add_range(my_id * 2, my_id * 2 + 2); + locally_owned.compress(); + + LinearAlgebra::TpetraWrappers::Vector v(locally_owned, MPI_COMM_WORLD); + + IndexSet workaround_set(n_procs * 2); + unsigned int my_first_index = (my_id * 2 + 2) % (n_procs * 2); + unsigned int my_second_index = my_id * 2 + 1; + workaround_set.add_index(my_first_index); + workaround_set.add_index(my_second_index); + workaround_set.add_index(0); + workaround_set.compress(); + LinearAlgebra::ReadWriteVector rw_vector(workaround_set); + + rw_vector(my_first_index) = my_id + 10; + rw_vector(my_second_index) = my_id + 100; + rw_vector(0) = 1.; + // rw_vector(2) = 1.; + v.import(rw_vector, VectorOperation::add); + deallog << "Tpetra first import add:" << std::endl; + v.print(deallog.get_file_stream()); + rw_vector.print(deallog.get_file_stream()); + + rw_vector(my_first_index) = my_id + 20; + rw_vector(my_second_index) = my_id + 200; + rw_vector(0) = 2.; + // rw_vector(2) = 3.; + v.import(rw_vector, VectorOperation::add); + deallog << "Tpetra second import add:" << std::endl; + v.print(deallog.get_file_stream()); + rw_vector.print(deallog.get_file_stream()); + + rw_vector.import(v, VectorOperation::add); + deallog << "ReadWrite import add:" << std::endl; + rw_vector.print(deallog.get_file_stream()); + + rw_vector(my_first_index) = my_id + 100; + rw_vector(my_second_index) = 1; + rw_vector(0) = 4.; + // rw_vector(2) = 3.; + rw_vector.import(v, VectorOperation::min); + deallog << "ReadWrite import min:" << std::endl; + rw_vector.print(deallog.get_file_stream()); + + rw_vector(my_first_index) = my_id + 100; + rw_vector(my_second_index) = 1; + rw_vector(0) = 4.; + // rw_vector(2) = 3.; + rw_vector.import(v, VectorOperation::max); + deallog << "ReadWrite import max:" << std::endl; + rw_vector.print(deallog.get_file_stream()); +} + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + MPILogInitAll log; + test(); +} diff --git a/tests/trilinos/tpetra_vector_03.mpirun=2.output b/tests/trilinos/tpetra_vector_03.mpirun=2.output new file mode 100644 index 0000000000..5d5f9c0afb --- /dev/null +++ b/tests/trilinos/tpetra_vector_03.mpirun=2.output @@ -0,0 +1,70 @@ + +DEAL:0::Tpetra first import add: +0->0: 2.000e+00 +1->1: 1.000e+02 + +IndexSet: {[0,2]} + +[0]: 1.000e+00 +[1]: 1.000e+02 +[2]: 1.000e+01 +DEAL:0::Tpetra second import add: +0->0: 4.000e+00 +1->1: 2.000e+02 + +IndexSet: {[0,2]} + +[0]: 2.000e+00 +[1]: 2.000e+02 +[2]: 2.000e+01 +DEAL:0::ReadWrite import add: +IndexSet: {[0,2]} + +[0]: 6.000e+00 +[1]: 4.000e+02 +[2]: 5.000e+01 +DEAL:0::ReadWrite import min: +IndexSet: {[0,2]} + +[0]: 4.000e+00 +[1]: 1.000e+00 +[2]: 3.000e+01 +DEAL:0::ReadWrite import max: +IndexSet: {[0,2]} + +[0]: 4.000e+00 +[1]: 2.000e+02 +[2]: 1.000e+02 + +DEAL:1::Tpetra first import add: +0->2: 1.000e+01 +1->3: 1.010e+02 + +IndexSet: {0, 3} + +[0]: 1.000e+00 +[3]: 1.010e+02 +DEAL:1::Tpetra second import add: +0->2: 3.000e+01 +1->3: 2.010e+02 + +IndexSet: {0, 3} + +[0]: 2.000e+00 +[3]: 2.010e+02 +DEAL:1::ReadWrite import add: +IndexSet: {0, 3} + +[0]: 6.000e+00 +[3]: 4.020e+02 +DEAL:1::ReadWrite import min: +IndexSet: {0, 3} + +[0]: 4.000e+00 +[3]: 1.000e+00 +DEAL:1::ReadWrite import max: +IndexSet: {0, 3} + +[0]: 4.000e+00 +[3]: 2.010e+02 + diff --git a/tests/trilinos/tpetra_vector_03.mpirun=4.output b/tests/trilinos/tpetra_vector_03.mpirun=4.output new file mode 100644 index 0000000000..bb8d36da88 --- /dev/null +++ b/tests/trilinos/tpetra_vector_03.mpirun=4.output @@ -0,0 +1,146 @@ + +DEAL:0::Tpetra first import add: +0->0: 4.000e+00 +1->1: 1.000e+02 + +IndexSet: {[0,2]} + +[0]: 1.000e+00 +[1]: 1.000e+02 +[2]: 1.000e+01 +DEAL:0::Tpetra second import add: +0->0: 8.000e+00 +1->1: 2.000e+02 + +IndexSet: {[0,2]} + +[0]: 2.000e+00 +[1]: 2.000e+02 +[2]: 2.000e+01 +DEAL:0::ReadWrite import add: +IndexSet: {[0,2]} + +[0]: 1.000e+01 +[1]: 4.000e+02 +[2]: 5.000e+01 +DEAL:0::ReadWrite import min: +IndexSet: {[0,2]} + +[0]: 4.000e+00 +[1]: 1.000e+00 +[2]: 3.000e+01 +DEAL:0::ReadWrite import max: +IndexSet: {[0,2]} + +[0]: 8.000e+00 +[1]: 2.000e+02 +[2]: 1.000e+02 + +DEAL:1::Tpetra first import add: +0->2: 1.000e+01 +1->3: 1.010e+02 + +IndexSet: {0, [3,4]} + +[0]: 1.000e+00 +[3]: 1.010e+02 +[4]: 1.100e+01 +DEAL:1::Tpetra second import add: +0->2: 3.000e+01 +1->3: 2.010e+02 + +IndexSet: {0, [3,4]} + +[0]: 2.000e+00 +[3]: 2.010e+02 +[4]: 2.100e+01 +DEAL:1::ReadWrite import add: +IndexSet: {0, [3,4]} + +[0]: 1.000e+01 +[3]: 4.020e+02 +[4]: 5.300e+01 +DEAL:1::ReadWrite import min: +IndexSet: {0, [3,4]} + +[0]: 4.000e+00 +[3]: 1.000e+00 +[4]: 3.200e+01 +DEAL:1::ReadWrite import max: +IndexSet: {0, [3,4]} + +[0]: 8.000e+00 +[3]: 2.010e+02 +[4]: 1.010e+02 + + +DEAL:2::Tpetra first import add: +0->4: 1.100e+01 +1->5: 1.020e+02 + +IndexSet: {0, [5,6]} + +[0]: 1.000e+00 +[5]: 1.020e+02 +[6]: 1.200e+01 +DEAL:2::Tpetra second import add: +0->4: 3.200e+01 +1->5: 2.020e+02 + +IndexSet: {0, [5,6]} + +[0]: 2.000e+00 +[5]: 2.020e+02 +[6]: 2.200e+01 +DEAL:2::ReadWrite import add: +IndexSet: {0, [5,6]} + +[0]: 1.000e+01 +[5]: 4.040e+02 +[6]: 5.600e+01 +DEAL:2::ReadWrite import min: +IndexSet: {0, [5,6]} + +[0]: 4.000e+00 +[5]: 1.000e+00 +[6]: 3.400e+01 +DEAL:2::ReadWrite import max: +IndexSet: {0, [5,6]} + +[0]: 8.000e+00 +[5]: 2.020e+02 +[6]: 1.020e+02 + + +DEAL:3::Tpetra first import add: +0->6: 1.200e+01 +1->7: 1.030e+02 + +IndexSet: {0, 7} + +[0]: 1.000e+00 +[7]: 1.030e+02 +DEAL:3::Tpetra second import add: +0->6: 3.400e+01 +1->7: 2.030e+02 + +IndexSet: {0, 7} + +[0]: 2.000e+00 +[7]: 2.030e+02 +DEAL:3::ReadWrite import add: +IndexSet: {0, 7} + +[0]: 1.000e+01 +[7]: 4.060e+02 +DEAL:3::ReadWrite import min: +IndexSet: {0, 7} + +[0]: 4.000e+00 +[7]: 1.000e+00 +DEAL:3::ReadWrite import max: +IndexSet: {0, 7} + +[0]: 8.000e+00 +[7]: 2.030e+02 + -- 2.39.5