From 1bbf0198f12dc8e0238883ec75cf90b4592fb26c Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 14 Sep 2018 13:59:27 +0200 Subject: [PATCH] Working --- cmake/config/template-arguments.in | 15 +- cmake/configure/configure_2_trilinos.cmake | 3 +- include/deal.II/dofs/dof_accessor.templates.h | 13 +- .../fe/fe_tools_extrapolate.templates.h | 6 +- .../fe/fe_tools_interpolate.templates.h | 10 +- include/deal.II/lac/read_write_vector.h | 12 +- .../deal.II/lac/read_write_vector.templates.h | 16 +- include/deal.II/lac/trilinos_sparse_matrix.h | 8 + include/deal.II/lac/trilinos_tpetra_vector.h | 80 +++--- include/deal.II/lac/vector_element_access.h | 102 +++++-- include/deal.II/multigrid/mg_transfer.h | 4 +- source/fe/mapping_q1_eulerian.cc | 3 +- source/lac/trilinos_sparse_matrix.cc | 84 ++++-- source/lac/trilinos_tpetra_vector.cc | 267 ++++++++++-------- tests/trilinos/tpetra_vector_01.cc | 21 +- tests/trilinos/tpetra_vector_02.cc | 26 +- tests/trilinos/tpetra_vector_03.cc | 9 +- .../trilinos/tpetra_vector_03.mpirun=2.output | 59 ++++ .../trilinos/tpetra_vector_03.mpirun=4.output | 123 ++++++++ 19 files changed, 607 insertions(+), 254 deletions(-) diff --git a/cmake/config/template-arguments.in b/cmake/config/template-arguments.in index 17a9f2ba94..8a58e90596 100644 --- a/cmake/config/template-arguments.in +++ b/cmake/config/template-arguments.in @@ -95,7 +95,8 @@ VECTOR_TYPES := { Vector; @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; @DEAL_II_EXPAND_EPETRA_VECTOR@; - @DEAL_II_EXPAND_TPETRA_VECTOR@; + @DEAL_II_EXPAND_TPETRA_VECTOR_DOUBLE@; + @DEAL_II_EXPAND_TPETRA_VECTOR_FLOAT@; @DEAL_II_EXPAND_PETSC_MPI_VECTOR@; @DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR@; @@ -121,7 +122,8 @@ REAL_VECTOR_TYPES := { Vector; @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; @DEAL_II_EXPAND_EPETRA_VECTOR@; - @DEAL_II_EXPAND_TPETRA_VECTOR@; + @DEAL_II_EXPAND_TPETRA_VECTOR_DOUBLE@; + @DEAL_II_EXPAND_TPETRA_VECTOR_FLOAT@; @DEAL_II_EXPAND_PETSC_MPI_VECTOR_REAL@; @DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR@; @@ -140,7 +142,8 @@ REAL_NONBLOCK_VECTORS := { Vector; @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; @DEAL_II_EXPAND_EPETRA_VECTOR@; - @DEAL_II_EXPAND_TPETRA_VECTOR@; + @DEAL_II_EXPAND_TPETRA_VECTOR_DOUBLE@; + @DEAL_II_EXPAND_TPETRA_VECTOR_FLOAT@; @DEAL_II_EXPAND_PETSC_MPI_VECTOR_REAL@; } @@ -148,7 +151,8 @@ 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_TPETRA_VECTOR_DOUBLE@; + @DEAL_II_EXPAND_TPETRA_VECTOR_FLOAT@; @DEAL_II_EXPAND_PETSC_MPI_VECTOR@; @DEAL_II_EXPAND_PETSC_MPI_BLOCKVECTOR@ } @@ -167,7 +171,8 @@ VECTORS_WITH_MATRIX := { Vector; @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; @DEAL_II_EXPAND_EPETRA_VECTOR@; - @DEAL_II_EXPAND_TPETRA_VECTOR@; + @DEAL_II_EXPAND_TPETRA_VECTOR_DOUBLE@; + @DEAL_II_EXPAND_TPETRA_VECTOR_FLOAT@; @DEAL_II_EXPAND_PETSC_MPI_VECTOR@; } diff --git a/cmake/configure/configure_2_trilinos.cmake b/cmake/configure/configure_2_trilinos.cmake index abe66905d8..661c739b7c 100644 --- a/cmake/configure/configure_2_trilinos.cmake +++ b/cmake/configure/configure_2_trilinos.cmake @@ -227,7 +227,8 @@ 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") + SET(DEAL_II_EXPAND_TPETRA_VECTOR_DOUBLE "LinearAlgebra::TpetraWrappers::Vector") + #SET(DEAL_II_EXPAND_TPETRA_VECTOR_FLOAT "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/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index 73b77c07a5..9b89e6e31f 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -1534,12 +1534,13 @@ namespace internal - template + 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) + 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); @@ -1551,7 +1552,7 @@ namespace internal 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); + LinearAlgebra::ReadWriteVector read_write_vector(index_set); read_write_vector.import(values, VectorOperation::insert); // Copy the elements from read_write_vector and reorder them. diff --git a/include/deal.II/fe/fe_tools_extrapolate.templates.h b/include/deal.II/fe/fe_tools_extrapolate.templates.h index 850a4c6216..ea76059356 100644 --- a/include/deal.II/fe/fe_tools_extrapolate.templates.h +++ b/include/deal.II/fe/fe_tools_extrapolate.templates.h @@ -1581,10 +1581,10 @@ namespace FETools # ifdef DEAL_II_WITH_MPI - template + template void - reinit_distributed(const DoFHandler & dh, - LinearAlgebra::TpetraWrappers::Vector &vector) + reinit_distributed(const DoFHandler & dh, + LinearAlgebra::TpetraWrappers::Vector &vector) { const parallel::distributed::Triangulation *parallel_tria = dynamic_cast< diff --git a/include/deal.II/fe/fe_tools_interpolate.templates.h b/include/deal.II/fe/fe_tools_interpolate.templates.h index e2e9b417fe..9699c0a9b4 100644 --- a/include/deal.II/fe/fe_tools_interpolate.templates.h +++ b/include/deal.II/fe/fe_tools_interpolate.templates.h @@ -513,17 +513,17 @@ namespace FETools AssertThrow(false, ExcNotImplemented()); } - template + template void back_interpolate( const DoFHandler &, const AffineConstraints< - typename LinearAlgebra::TpetraWrappers::Vector::value_type> &, - const LinearAlgebra::TpetraWrappers::Vector &, + typename LinearAlgebra::TpetraWrappers::Vector::value_type> &, + const LinearAlgebra::TpetraWrappers::Vector &, const DoFHandler &, const AffineConstraints< - typename LinearAlgebra::TpetraWrappers::Vector::value_type> &, - LinearAlgebra::TpetraWrappers::Vector &) + typename LinearAlgebra::TpetraWrappers::Vector::value_type> &, + LinearAlgebra::TpetraWrappers::Vector &) { AssertThrow(false, ExcNotImplemented()); } diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h index 70d89ade7a..0ed6f9d829 100644 --- a/include/deal.II/lac/read_write_vector.h +++ b/include/deal.II/lac/read_write_vector.h @@ -344,8 +344,8 @@ namespace LinearAlgebra * performance. */ void - import(const TpetraWrappers::Vector &tpetra_vec, - VectorOperation::values operation, + import(const TpetraWrappers::Vector &tpetra_vec, + VectorOperation::values operation, const std::shared_ptr &communication_pattern = std::shared_ptr()); @@ -615,10 +615,10 @@ namespace LinearAlgebra * used directly. */ void - import(const Tpetra::Vector<> &tpetra_vector, - const IndexSet & locally_owned_elements, - VectorOperation::values operation, - const MPI_Comm & mpi_comm, + 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); diff --git a/include/deal.II/lac/read_write_vector.templates.h b/include/deal.II/lac/read_write_vector.templates.h index b4437cfd3f..87139aedd7 100644 --- a/include/deal.II/lac/read_write_vector.templates.h +++ b/include/deal.II/lac/read_write_vector.templates.h @@ -509,10 +509,10 @@ namespace LinearAlgebra template void ReadWriteVector::import( - const Tpetra::Vector<> &vector, - const IndexSet & source_elements, - VectorOperation::values operation, - const MPI_Comm & mpi_comm, + const Tpetra::Vector &vector, + const IndexSet & source_elements, + VectorOperation::values operation, + const MPI_Comm & mpi_comm, const std::shared_ptr &communication_pattern) { @@ -554,13 +554,13 @@ namespace LinearAlgebra Tpetra::Export<> tpetra_export(tpetra_comm_pattern->get_tpetra_export()); - Tpetra::Vector<> target_vector(tpetra_export.getSourceMap()); + 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; + using size_type = typename std::decay::type; Assert(size == 0 || values != nullptr, ExcInternalError("Export failed.")); AssertDimension(size, stored_elements.n_elements()); @@ -741,8 +741,8 @@ namespace LinearAlgebra template void ReadWriteVector::import( - const LinearAlgebra::TpetraWrappers::Vector &trilinos_vec, - VectorOperation::values operation, + const LinearAlgebra::TpetraWrappers::Vector &trilinos_vec, + VectorOperation::values operation, const std::shared_ptr &communication_pattern) { diff --git a/include/deal.II/lac/trilinos_sparse_matrix.h b/include/deal.II/lac/trilinos_sparse_matrix.h index 006dc31712..b158a224d4 100644 --- a/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_sparse_matrix.h @@ -1317,6 +1317,14 @@ namespace TrilinosWrappers * initialize the matrix with a sparsity pattern to fix the matrix * structure before inserting elements. */ + template + void + set(const size_type row, + const size_type n_cols, + const size_type *col_indices, + const Number * values, + const bool elide_zero_values = false); + void set(const size_type row, const size_type n_cols, diff --git a/include/deal.II/lac/trilinos_tpetra_vector.h b/include/deal.II/lac/trilinos_tpetra_vector.h index df0b0c99c5..798248703a 100644 --- a/include/deal.II/lac/trilinos_tpetra_vector.h +++ b/include/deal.II/lac/trilinos_tpetra_vector.h @@ -51,16 +51,20 @@ namespace LinearAlgebra /** * 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. + * LinearAlgebra::VectorSpaceVector class. * * @ingroup TrilinosWrappers * @ingroup Vectors * @author Daniel Arndt, 2018 */ - class Vector : public VectorSpaceVector, public Subscriptor + template + class Vector : public VectorSpaceVector, public Subscriptor { public: + using value_type = Number; + + using size_type = typename VectorSpaceVector::size_type; + /** * Constructor. Create a vector of dimension zero. */ @@ -97,7 +101,7 @@ namespace LinearAlgebra * copied. */ virtual void - reinit(const VectorSpaceVector &V, + reinit(const VectorSpaceVector &V, const bool omit_zeroing_entries = false) override; /** @@ -113,7 +117,7 @@ namespace LinearAlgebra * only allowed if @p s is equal to zero. */ virtual Vector & - operator=(const double s) override; + operator=(const Number s) override; /** * Imports all the elements present in the vector's IndexSet from the @@ -126,7 +130,7 @@ namespace LinearAlgebra */ virtual void import( - const ReadWriteVector & V, + const ReadWriteVector & V, VectorOperation::values operation, std::shared_ptr communication_pattern = std::shared_ptr()) override; @@ -135,64 +139,64 @@ namespace LinearAlgebra * Multiply the entire vector by a fixed factor. */ virtual Vector & - operator*=(const double factor) override; + operator*=(const Number factor) override; /** * Divide the entire vector by a fixed factor. */ virtual Vector & - operator/=(const double factor) override; + operator/=(const Number factor) override; /** * Add the vector @p V to the present one. */ virtual Vector & - operator+=(const VectorSpaceVector &V) override; + operator+=(const VectorSpaceVector &V) override; /** * Subtract the vector @p V from the present one. */ virtual Vector & - operator-=(const VectorSpaceVector &V) override; + 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; + virtual Number + 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; + add(const Number 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; + add(const Number 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; + add(const Number a, + const VectorSpaceVector &V, + const Number 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; + sadd(const Number s, + const Number a, + const VectorSpaceVector &V) override; /** * Scale each element of this vector by the corresponding element in the @@ -201,13 +205,13 @@ namespace LinearAlgebra * vectors need to have the same layout. */ virtual void - scale(const VectorSpaceVector &scaling_factors) override; + scale(const VectorSpaceVector &scaling_factors) override; /** * Assignment *this = a*V. */ virtual void - equ(const double a, const VectorSpaceVector &V) override; + equ(const Number a, const VectorSpaceVector &V) override; /** * Return whether the vector contains only elements with value zero. @@ -218,28 +222,28 @@ namespace LinearAlgebra /** * Return the mean value of the element of this vector. */ - virtual double + virtual Number 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 + virtual typename LinearAlgebra::VectorSpaceVector::real_type 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 + virtual typename LinearAlgebra::VectorSpaceVector::real_type 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 + virtual typename LinearAlgebra::VectorSpaceVector::real_type linfty_norm() const override; /** @@ -264,10 +268,10 @@ namespace LinearAlgebra * 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; + virtual Number + add_and_dot(const Number a, + const VectorSpaceVector &V, + const VectorSpaceVector &W) override; /** * This function always returns false and is present only for backward * compatibility. @@ -306,14 +310,14 @@ namespace LinearAlgebra * Return a const reference to the underlying Trilinos * Tpetra::Vector class. */ - const Tpetra::Vector<> & + const Tpetra::Vector & trilinos_vector() const; /** * Return a (modifyable) reference to the underlying Trilinos * Tpetra::Vector class. */ - Tpetra::Vector<> & + Tpetra::Vector & trilinos_vector(); /** @@ -367,7 +371,7 @@ namespace LinearAlgebra /** * Pointer to the actual Tpetra vector object. */ - std::unique_ptr> vector; + std::unique_ptr> vector; /** * IndexSet of the elements of the last imported vector. @@ -383,8 +387,9 @@ namespace LinearAlgebra }; + template inline bool - Vector::has_ghost_elements() const + Vector::has_ghost_elements() const { return false; } @@ -395,8 +400,9 @@ namespace LinearAlgebra /** * Declare dealii::LinearAlgebra::TpetraWrappers::Vector as distributed vector. */ -template <> -struct is_serial_vector : std::false_type +template +struct is_serial_vector> + : std::false_type {}; DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/lac/vector_element_access.h b/include/deal.II/lac/vector_element_access.h index b168260588..22e541033c 100644 --- a/include/deal.II/lac/vector_element_access.h +++ b/include/deal.II/lac/vector_element_access.h @@ -127,13 +127,13 @@ namespace internal template <> inline void - ElementAccess::add( - const double value, - const types::global_dof_index i, - LinearAlgebra::TpetraWrappers::Vector &V) + 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(); + Tpetra::Vector vector = V.trilinos_vector(); TrilinosWrappers::types::int_type trilinos_i = vector.getMap()->getLocalElement( static_cast(i)); @@ -144,20 +144,68 @@ namespace internal // We're going to modify the data on host. vector.modify(); vector_1d(trilinos_i) += value; - vector.sync::device_type::memory_space>(); + vector.sync::device_type::memory_space>(); } template <> inline void - ElementAccess::set( - const double value, - const types::global_dof_index i, - LinearAlgebra::TpetraWrappers::Vector &V) + ElementAccess>::add( + const float 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 void + ElementAccess>::set( + const float value, + const types::global_dof_index i, + LinearAlgebra::TpetraWrappers::Vector &V) { // Extract local indices in the vector. - Tpetra::Vector<> vector = V.trilinos_vector(); + Tpetra::Vector vector = V.trilinos_vector(); TrilinosWrappers::types::int_type trilinos_i = vector.getMap()->getLocalElement( static_cast(i)); @@ -168,18 +216,40 @@ namespace internal // We're going to modify the data on host. vector.modify(); vector_1d(trilinos_i) = value; - vector.sync::device_type::memory_space>(); + vector.sync::device_type::memory_space>(); } + template <> inline double - ElementAccess::get( - const LinearAlgebra::TpetraWrappers::Vector &V, - const types::global_dof_index i) + 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); + } + + + + template <> + inline float + 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(); + Tpetra::Vector vector = V.trilinos_vector(); TrilinosWrappers::types::int_type trilinos_i = vector.getMap()->getLocalElement( static_cast(i)); diff --git a/include/deal.II/multigrid/mg_transfer.h b/include/deal.II/multigrid/mg_transfer.h index e4b4f3db2a..961b0223df 100644 --- a/include/deal.II/multigrid/mg_transfer.h +++ b/include/deal.II/multigrid/mg_transfer.h @@ -132,8 +132,8 @@ namespace internal }; # ifdef DEAL_II_WITH_MPI - template <> - struct MatrixSelector + template + struct MatrixSelector> { using Sparsity = ::dealii::TrilinosWrappers::SparsityPattern; using Matrix = ::dealii::TrilinosWrappers::SparseMatrix; diff --git a/source/fe/mapping_q1_eulerian.cc b/source/fe/mapping_q1_eulerian.cc index 666a7bef6e..3c99fb8c83 100644 --- a/source/fe/mapping_q1_eulerian.cc +++ b/source/fe/mapping_q1_eulerian.cc @@ -77,7 +77,8 @@ MappingQ1Eulerian::get_vertices( Assert(dof_cell->active() == true, ExcInactiveCell()); // now get the values of the shift vectors at the vertices - Vector mapping_values(shiftmap_dof_handler->get_fe().dofs_per_cell); + Vector mapping_values( + shiftmap_dof_handler->get_fe().dofs_per_cell); dof_cell->get_dof_values(*euler_transform_vectors, mapping_values); for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; ++i) diff --git a/source/lac/trilinos_sparse_matrix.cc b/source/lac/trilinos_sparse_matrix.cc index 350006da3c..0996f76053 100644 --- a/source/lac/trilinos_sparse_matrix.cc +++ b/source/lac/trilinos_sparse_matrix.cc @@ -48,28 +48,28 @@ namespace TrilinosWrappers namespace internal { template - double * + typename VectorType::value_type * begin(VectorType &V) { return V.begin(); } template - const double * + const typename VectorType::value_type * begin(const VectorType &V) { return V.begin(); } template - double * + typename VectorType::value_type * end(VectorType &V) { return V.end(); } template - const double * + const typename VectorType::value_type * end(const VectorType &V) { return V.end(); @@ -104,31 +104,31 @@ namespace TrilinosWrappers return V.trilinos_vector()[0] + V.trilinos_vector().MyLength(); } - template <> - double * - begin(LinearAlgebra::TpetraWrappers::Vector &V) + template + Number * + begin(LinearAlgebra::TpetraWrappers::Vector &V) { return V.trilinos_vector().getDataNonConst().get(); } - template <> - const double * - begin(const LinearAlgebra::TpetraWrappers::Vector &V) + template + const Number * + begin(const LinearAlgebra::TpetraWrappers::Vector &V) { return V.trilinos_vector().getData().get(); } - template <> - double * - end(LinearAlgebra::TpetraWrappers::Vector &V) + template + Number * + end(LinearAlgebra::TpetraWrappers::Vector &V) { return V.trilinos_vector().getDataNonConst().get() + V.trilinos_vector().getLocalLength(); } - template <> - const double * - end(const LinearAlgebra::TpetraWrappers::Vector &V) + template + const Number * + end(const LinearAlgebra::TpetraWrappers::Vector &V) { return V.trilinos_vector().getData().get() + V.trilinos_vector().getLocalLength(); @@ -1540,6 +1540,22 @@ namespace TrilinosWrappers + template + void + SparseMatrix::set(const size_type row, + const size_type n_cols, + const size_type *col_indices, + const Number * values, + const bool elide_zero_values) + { + std::vector trilinos_values(n_cols); + std::copy(values, values + n_cols, trilinos_values.begin()); + this->set( + row, n_cols, col_indices, trilinos_values.data(), elide_zero_values); + } + + + void SparseMatrix::set(const size_type row, const size_type n_cols, @@ -3476,78 +3492,98 @@ namespace TrilinosWrappers template void SparseMatrix::vmult(MPI::Vector &, const MPI::Vector &) const; + template void SparseMatrix::vmult(dealii::Vector &, const dealii::Vector &) const; + template void SparseMatrix::vmult( 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; + dealii::LinearAlgebra::TpetraWrappers::Vector &, + const dealii::LinearAlgebra::TpetraWrappers::Vector &) const; + + /*template void + SparseMatrix::vmult( + dealii::LinearAlgebra::TpetraWrappers::Vector &, + const dealii::LinearAlgebra::TpetraWrappers::Vector &) const;*/ template void SparseMatrix::vmult( dealii::LinearAlgebra::EpetraWrappers::Vector &, const dealii::LinearAlgebra::EpetraWrappers::Vector &) const; # endif + template void SparseMatrix::Tvmult(MPI::Vector &, const MPI::Vector &) const; + template void SparseMatrix::Tvmult(dealii::Vector &, const dealii::Vector &) const; + template void SparseMatrix::Tvmult( 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; + dealii::LinearAlgebra::TpetraWrappers::Vector &, + const dealii::LinearAlgebra::TpetraWrappers::Vector &) const; template void SparseMatrix::Tvmult( dealii::LinearAlgebra::EpetraWrappers::Vector &, const dealii::LinearAlgebra::EpetraWrappers::Vector &) const; # endif + template void SparseMatrix::vmult_add(MPI::Vector &, const MPI::Vector &) const; + template void SparseMatrix::vmult_add(dealii::Vector &, const dealii::Vector &) const; + template void SparseMatrix::vmult_add( 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; + dealii::LinearAlgebra::TpetraWrappers::Vector &, + const dealii::LinearAlgebra::TpetraWrappers::Vector &) const; template void SparseMatrix::vmult_add( dealii::LinearAlgebra::EpetraWrappers::Vector &, const dealii::LinearAlgebra::EpetraWrappers::Vector &) const; # endif + template void SparseMatrix::Tvmult_add(MPI::Vector &, const MPI::Vector &) const; + template void SparseMatrix::Tvmult_add(dealii::Vector &, const dealii::Vector &) const; + template void SparseMatrix::Tvmult_add( 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; + dealii::LinearAlgebra::TpetraWrappers::Vector &, + const dealii::LinearAlgebra::TpetraWrappers::Vector &) const; template void SparseMatrix::Tvmult_add( diff --git a/source/lac/trilinos_tpetra_vector.cc b/source/lac/trilinos_tpetra_vector.cc index 5e3c00b53a..88324e6d53 100644 --- a/source/lac/trilinos_tpetra_vector.cc +++ b/source/lac/trilinos_tpetra_vector.cc @@ -40,37 +40,41 @@ namespace LinearAlgebra { namespace TpetraWrappers { - Vector::Vector() - : vector(new Tpetra::Vector<>(Teuchos::RCP>( + template + Vector::Vector() + : vector(new Tpetra::Vector(Teuchos::RCP>( new Tpetra::Map<>(0, 0, Utilities::Trilinos::tpetra_comm_self())))) {} - Vector::Vector(const Vector &V) + template + Vector::Vector(const Vector &V) : Subscriptor() - , vector(new Tpetra::Vector<>(V.trilinos_vector(), Teuchos::Copy)) + , 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<>( + template + 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))))) {} + template void - Vector::reinit(const IndexSet ¶llel_partitioner, - const MPI_Comm &communicator, - const bool omit_zeroing_entries) + 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>( + vector = std_cxx14::make_unique>( Teuchos::rcp(new Tpetra::Map<>(input_map))); else if (omit_zeroing_entries == false) { @@ -80,16 +84,17 @@ namespace LinearAlgebra + template void - Vector::reinit(const VectorSpaceVector &V, - const bool omit_zeroing_entries) + Vector::reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries) { // Check that casting will work. - Assert(dynamic_cast(&V) != nullptr, + Assert(dynamic_cast *>(&V) != nullptr, ExcVectorTypeNotCompatible()); // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast(V); + const Vector &down_V = dynamic_cast &>(V); reinit(down_V.locally_owned_elements(), down_V.get_mpi_communicator(), @@ -98,8 +103,9 @@ namespace LinearAlgebra - Vector & - Vector::operator=(const Vector &V) + template + Vector & + Vector::operator=(const Vector &V) { // Distinguish three cases: // - First case: both vectors have the same layout. @@ -119,8 +125,8 @@ namespace LinearAlgebra Tpetra::REPLACE); } else - vector = - std_cxx14::make_unique>(V.trilinos_vector()); + vector = std_cxx14::make_unique>( + V.trilinos_vector()); } return *this; @@ -128,10 +134,12 @@ namespace LinearAlgebra - Vector & - Vector::operator=(const double s) + template + Vector & + Vector::operator=(const Number s) { - Assert(s == 0., ExcMessage("Only 0 can be assigned to a vector.")); + Assert(s == Number(0.), + ExcMessage("Only 0 can be assigned to a vector.")); vector->putScalar(s); @@ -140,9 +148,10 @@ namespace LinearAlgebra + template void - Vector::import( - const ReadWriteVector & V, + Vector::import( + const ReadWriteVector & V, VectorOperation::values operation, std::shared_ptr communication_pattern) { @@ -177,17 +186,18 @@ namespace LinearAlgebra } Tpetra::Export<> tpetra_export(tpetra_comm_pattern->get_tpetra_export()); - Tpetra::Vector<> source_vector(tpetra_export.getSourceMap()); + Tpetra::Vector source_vector(tpetra_export.getSourceMap()); - source_vector.sync(); - auto x_2d = source_vector.getLocalView(); + source_vector.template sync(); + auto x_2d = source_vector.template getLocalView(); auto x_1d = Kokkos::subview(x_2d, Kokkos::ALL(), 0); - source_vector.modify(); + source_vector.template 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>(); + source_vector.template sync< + typename Tpetra::Vector::device_type::memory_space>(); if (operation == VectorOperation::insert) vector->doExport(source_vector, tpetra_export, Tpetra::REPLACE); else if (operation == VectorOperation::add) @@ -198,8 +208,9 @@ namespace LinearAlgebra - Vector & - Vector::operator*=(const double factor) + template + Vector & + Vector::operator*=(const Number factor) { AssertIsFinite(factor); vector->scale(factor); @@ -209,27 +220,29 @@ namespace LinearAlgebra - Vector & - Vector::operator/=(const double factor) + template + Vector & + Vector::operator/=(const Number factor) { AssertIsFinite(factor); - Assert(factor != 0., ExcZero()); - *this *= 1. / factor; + Assert(factor != Number(0.), ExcZero()); + *this *= Number(1.) / factor; return *this; } - Vector & - Vector::operator+=(const VectorSpaceVector &V) + template + Vector & + Vector::operator+=(const VectorSpaceVector &V) { // Check that casting will work. - Assert(dynamic_cast(&V) != nullptr, + Assert(dynamic_cast *>(&V) != nullptr, ExcVectorTypeNotCompatible()); // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast(V); + 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()))) { @@ -247,8 +260,8 @@ namespace LinearAlgebra data_exchange, Tpetra::ADD);*/ - Tpetra::Vector<> dummy(vector->getMap(), false); - Tpetra::Import<> data_exchange(dummy.getMap(), + Tpetra::Vector dummy(vector->getMap(), false); + Tpetra::Import<> data_exchange(dummy.getMap(), down_V.trilinos_vector().getMap()); dummy.doExport(down_V.trilinos_vector(), @@ -263,8 +276,9 @@ namespace LinearAlgebra - Vector & - Vector::operator-=(const VectorSpaceVector &V) + template + Vector & + Vector::operator-=(const VectorSpaceVector &V) { this->add(-1., V); @@ -273,14 +287,15 @@ namespace LinearAlgebra - double Vector::operator*(const VectorSpaceVector &V) const + template + Number Vector::operator*(const VectorSpaceVector &V) const { // Check that casting will work. - Assert(dynamic_cast(&V) != nullptr, + Assert(dynamic_cast *>(&V) != nullptr, ExcVectorTypeNotCompatible()); // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast(V); + 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()), @@ -291,34 +306,37 @@ namespace LinearAlgebra + template void - Vector::add(const double a) + Vector::add(const Number a) { AssertIsFinite(a); - vector->sync(); - auto vector_2d = vector->getLocalView(); + vector->template sync(); + auto vector_2d = vector->template getLocalView(); auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); - vector->modify(); + vector->template 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>(); + vector->template sync< + typename Tpetra::Vector::device_type::memory_space>(); } + template void - Vector::add(const double a, const VectorSpaceVector &V) + Vector::add(const Number a, const VectorSpaceVector &V) { // Check that casting will work. - Assert(dynamic_cast(&V) != nullptr, + Assert(dynamic_cast *>(&V) != nullptr, ExcVectorTypeNotCompatible()); // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast(V); + const Vector &down_V = dynamic_cast &>(V); AssertIsFinite(a); Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())), ExcDifferentParallelPartitioning()); @@ -328,23 +346,24 @@ namespace LinearAlgebra + template void - Vector::add(const double a, - const VectorSpaceVector &V, - const double b, - const VectorSpaceVector &W) + Vector::add(const Number a, + const VectorSpaceVector &V, + const Number b, + const VectorSpaceVector &W) { // Check that casting will work. - Assert(dynamic_cast(&V) != nullptr, + Assert(dynamic_cast *>(&V) != nullptr, ExcVectorTypeNotCompatible()); // Check that casting will work. - Assert(dynamic_cast(&W) != nullptr, + Assert(dynamic_cast *>(&W) != nullptr, ExcVectorTypeNotCompatible()); // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast(V); + const Vector &down_V = dynamic_cast &>(V); // Downcast W. If fails, throws an exception. - const Vector &down_W = dynamic_cast(W); + 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())), @@ -358,35 +377,37 @@ namespace LinearAlgebra + template void - Vector::sadd(const double s, - const double a, - const VectorSpaceVector &V) + Vector::sadd(const Number s, + const Number a, + const VectorSpaceVector &V) { // Check that casting will work. - Assert(dynamic_cast(&V) != nullptr, + 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); + const Vector &down_V = dynamic_cast &>(V); + Vector tmp(down_V); tmp *= a; *this += tmp; } + template void - Vector::scale(const VectorSpaceVector &scaling_factors) + Vector::scale(const VectorSpaceVector &scaling_factors) { // Check that casting will work. - Assert(dynamic_cast(&scaling_factors) != nullptr, + Assert(dynamic_cast *>(&scaling_factors) != nullptr, ExcVectorTypeNotCompatible()); // Downcast scaling_factors. If fails, throws an exception. - const Vector &down_scaling_factors = - dynamic_cast(scaling_factors); + const Vector &down_scaling_factors = + dynamic_cast &>(scaling_factors); Assert(vector->getMap()->isSameAs( *(down_scaling_factors.trilinos_vector().getMap())), ExcDifferentParallelPartitioning()); @@ -399,15 +420,16 @@ namespace LinearAlgebra + template void - Vector::equ(const double a, const VectorSpaceVector &V) + Vector::equ(const Number a, const VectorSpaceVector &V) { // Check that casting will work. - Assert(dynamic_cast(&V) != nullptr, + Assert(dynamic_cast *>(&V) != nullptr, ExcVectorTypeNotCompatible()); // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast(V); + 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) @@ -421,18 +443,19 @@ namespace LinearAlgebra + template bool - Vector::all_zero() const + 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, + Number * start_ptr = vector->getDataNonConst().get(); + const Number *ptr = start_ptr, *eptr = start_ptr + vector->getLocalLength(); unsigned int flag = 0; while (ptr != eptr) { - if (*ptr != 0) + if (*ptr != Number(0)) { flag = 1; break; @@ -453,42 +476,47 @@ namespace LinearAlgebra - double - Vector::mean_value() const + template + Number + Vector::mean_value() const { return vector->meanValue(); } - double - Vector::l1_norm() const + template + typename LinearAlgebra::VectorSpaceVector::real_type + Vector::l1_norm() const { return vector->norm1(); } - double - Vector::l2_norm() const + template + typename LinearAlgebra::VectorSpaceVector::real_type + Vector::l2_norm() const { return vector->norm2(); } - double - Vector::linfty_norm() const + template + typename LinearAlgebra::VectorSpaceVector::real_type + Vector::linfty_norm() const { return vector->normInf(); } - double - Vector::add_and_dot(const double a, - const VectorSpaceVector &V, - const VectorSpaceVector &W) + template + Number + Vector::add_and_dot(const Number a, + const VectorSpaceVector &V, + const VectorSpaceVector &W) { this->add(a, V); @@ -497,16 +525,18 @@ namespace LinearAlgebra - Vector::size_type - Vector::size() const + template + typename Vector::size_type + Vector::size() const { return vector->getGlobalLength(); } + template MPI_Comm - Vector::get_mpi_communicator() const + Vector::get_mpi_communicator() const { const auto tpetra_comm = dynamic_cast *>( vector->getMap()->getComm().get()); @@ -516,21 +546,17 @@ namespace LinearAlgebra + template ::dealii::IndexSet - Vector::locally_owned_elements() const + 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) { @@ -546,27 +572,30 @@ namespace LinearAlgebra - const Tpetra::Vector<> & - Vector::trilinos_vector() const + template + const Tpetra::Vector & + Vector::trilinos_vector() const { return *vector; } - Tpetra::Vector<> & - Vector::trilinos_vector() + template + Tpetra::Vector & + Vector::trilinos_vector() { return *vector; } + template void - Vector::print(std::ostream & out, - const unsigned int precision, - const bool scientific, - const bool across) const + 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); @@ -581,9 +610,9 @@ namespace LinearAlgebra 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); + vector->template sync(); + auto vector_2d = vector->template getLocalView(); + auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); const size_t local_length = vector->getLocalLength(); if (across) @@ -601,25 +630,33 @@ namespace LinearAlgebra + template std::size_t - Vector::memory_consumption() const + Vector::memory_consumption() const { return sizeof(*this) + vector->getLocalLength() * - (sizeof(double) + sizeof(TrilinosWrappers::types::int_type)); + (sizeof(Number) + sizeof(TrilinosWrappers::types::int_type)); } + template void - Vector::create_tpetra_comm_pattern(const IndexSet &source_index_set, - const MPI_Comm &mpi_comm) + 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); } + + template class Vector; + template class Vector; + template class Vector>; + template class Vector>; + } // namespace TpetraWrappers } // namespace LinearAlgebra diff --git a/tests/trilinos/tpetra_vector_01.cc b/tests/trilinos/tpetra_vector_01.cc index c10ac87bb1..ba39dfa34e 100644 --- a/tests/trilinos/tpetra_vector_01.cc +++ b/tests/trilinos/tpetra_vector_01.cc @@ -26,7 +26,7 @@ // Check LinearAlgebra::TpetraWrappers::Vector assignment and import - +template void test() { @@ -45,10 +45,10 @@ test() } 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); + 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.")); @@ -69,9 +69,9 @@ test() 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); + 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) @@ -204,7 +204,7 @@ test() b.import(read_write_1, VectorOperation::insert); c.import(read_write_1, VectorOperation::insert); - const double val = b * c; + const Number val = b * c; AssertThrow(val == 285., ExcMessage("Problem in operator *.")); } @@ -217,7 +217,8 @@ main(int argc, char **argv) Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); - test(); + test(); + test(); deallog << "OK" << std::endl; diff --git a/tests/trilinos/tpetra_vector_02.cc b/tests/trilinos/tpetra_vector_02.cc index 8e48753530..a775c7ec30 100644 --- a/tests/trilinos/tpetra_vector_02.cc +++ b/tests/trilinos/tpetra_vector_02.cc @@ -26,6 +26,7 @@ // Check LinearAlgebra::TpetraWrappers::Vector add and sadd. +template void test() { @@ -44,12 +45,12 @@ test() } 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); + 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) @@ -58,9 +59,9 @@ test() 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); + 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) @@ -188,14 +189,14 @@ test() AssertThrow(b.l1_norm() == 95., ExcMessage("Problem in l1_norm.")); - const double eps = 1e-6; + const Number 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); + const Number val = a.add_and_dot(2., a, b); AssertThrow(val == 1530., ExcMessage("Problem in add_and_dot")); } @@ -208,7 +209,8 @@ main(int argc, char **argv) Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); - test(); + test(); + test(); deallog << "OK" << std::endl; diff --git a/tests/trilinos/tpetra_vector_03.cc b/tests/trilinos/tpetra_vector_03.cc index 4927aaedb8..650296a491 100644 --- a/tests/trilinos/tpetra_vector_03.cc +++ b/tests/trilinos/tpetra_vector_03.cc @@ -28,6 +28,7 @@ #include "../tests.h" +template void test() { @@ -38,7 +39,8 @@ test() locally_owned.add_range(my_id * 2, my_id * 2 + 2); locally_owned.compress(); - LinearAlgebra::TpetraWrappers::Vector v(locally_owned, MPI_COMM_WORLD); + 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); @@ -47,7 +49,7 @@ test() workaround_set.add_index(my_second_index); workaround_set.add_index(0); workaround_set.compress(); - LinearAlgebra::ReadWriteVector rw_vector(workaround_set); + LinearAlgebra::ReadWriteVector rw_vector(workaround_set); rw_vector(my_first_index) = my_id + 10; rw_vector(my_second_index) = my_id + 100; @@ -97,5 +99,6 @@ main(int argc, char **argv) argc, argv, testing_max_num_threads()); MPILogInitAll log; - test(); + test(); + test(); } diff --git a/tests/trilinos/tpetra_vector_03.mpirun=2.output b/tests/trilinos/tpetra_vector_03.mpirun=2.output index 638319c2f1..e94f192176 100644 --- a/tests/trilinos/tpetra_vector_03.mpirun=2.output +++ b/tests/trilinos/tpetra_vector_03.mpirun=2.output @@ -31,7 +31,66 @@ IndexSet: {[0,2]} [0]: 4.000e+00 [1]: 2.000e+02 [2]: 1.000e+02 +DEAL:0::Tpetra first import add: +2.000e+00 1.000e+02 +IndexSet: {[0,2]} + +[0]: 1.000e+00 +[1]: 1.000e+02 +[2]: 1.000e+01 +DEAL:0::Tpetra second import add: +4.000e+00 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: +1.000e+01 1.010e+02 +IndexSet: {0, 3} + +[0]: 1.000e+00 +[3]: 1.010e+02 +DEAL:1::Tpetra second import add: +3.000e+01 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 DEAL:1::Tpetra first import add: 1.000e+01 1.010e+02 IndexSet: {0, 3} diff --git a/tests/trilinos/tpetra_vector_03.mpirun=4.output b/tests/trilinos/tpetra_vector_03.mpirun=4.output index f70ed13b09..5fbf8dea2d 100644 --- a/tests/trilinos/tpetra_vector_03.mpirun=4.output +++ b/tests/trilinos/tpetra_vector_03.mpirun=4.output @@ -28,6 +28,38 @@ IndexSet: {[0,2]} DEAL:0::ReadWrite import max: IndexSet: {[0,2]} +[0]: 8.000e+00 +[1]: 2.000e+02 +[2]: 1.000e+02 +DEAL:0::Tpetra first import add: +4.000e+00 1.000e+02 +IndexSet: {[0,2]} + +[0]: 1.000e+00 +[1]: 1.000e+02 +[2]: 1.000e+01 +DEAL:0::Tpetra second import add: +8.000e+00 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 @@ -61,6 +93,38 @@ IndexSet: {0, [3,4]} DEAL:1::ReadWrite import max: IndexSet: {0, [3,4]} +[0]: 8.000e+00 +[3]: 2.010e+02 +[4]: 1.010e+02 +DEAL:1::Tpetra first import add: +1.000e+01 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: +3.000e+01 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 @@ -98,8 +162,67 @@ IndexSet: {0, [5,6]} [0]: 8.000e+00 [5]: 2.020e+02 [6]: 1.020e+02 +DEAL:2::Tpetra first import add: +1.100e+01 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: +3.200e+01 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: +1.200e+01 1.030e+02 +IndexSet: {0, 7} + +[0]: 1.000e+00 +[7]: 1.030e+02 +DEAL:3::Tpetra second import add: +3.400e+01 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 DEAL:3::Tpetra first import add: 1.200e+01 1.030e+02 IndexSet: {0, 7} -- 2.39.5