From 8ba44d267d2aef3eef1e83fb9ce7d84bcc266049 Mon Sep 17 00:00:00 2001 From: Daniel Jodlbauer Date: Tue, 5 Nov 2019 17:40:23 +0100 Subject: [PATCH] Add SIMD type to template list of MF-Operators --- .../changes/minor/20191106DanielJodlbauer | 3 + include/deal.II/matrix_free/operators.h | 516 +++++++++++------- 2 files changed, 321 insertions(+), 198 deletions(-) create mode 100644 doc/news/changes/minor/20191106DanielJodlbauer diff --git a/doc/news/changes/minor/20191106DanielJodlbauer b/doc/news/changes/minor/20191106DanielJodlbauer new file mode 100644 index 0000000000..a9cdb8cb98 --- /dev/null +++ b/doc/news/changes/minor/20191106DanielJodlbauer @@ -0,0 +1,3 @@ +New: template argument for type of vectorized array (vector width) in matrix-free operators. +
+(Daniel Jodlbauer, 2019/11/06) diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 1906817978..006b0606c6 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -183,7 +183,9 @@ namespace MatrixFreeOperators * @author Denis Davydov, Daniel Arndt, Martin Kronbichler, 2016, 2017 */ template > + typename VectorType = LinearAlgebra::distributed::Vector, + typename VectorizedArrayType = + VectorizedArray> class Base : public Subscriptor { public: @@ -231,7 +233,8 @@ namespace MatrixFreeOperators * selection, defining a diagonal block. */ void - initialize(std::shared_ptr> data, + initialize(std::shared_ptr< + const MatrixFree> data, const std::vector &selected_row_blocks = std::vector(), const std::vector &selected_column_blocks = @@ -251,7 +254,8 @@ namespace MatrixFreeOperators * empty, all components are selected. */ void - initialize(std::shared_ptr> data, + initialize(std::shared_ptr< + const MatrixFree> data, const MGConstrainedDoFs & mg_constrained_dofs, const unsigned int level, const std::vector &selected_row_blocks = @@ -272,7 +276,8 @@ namespace MatrixFreeOperators * empty, all components are selected. */ void - initialize(std::shared_ptr> data_, + initialize(std::shared_ptr< + const MatrixFree> data_, const std::vector &mg_constrained_dofs, const unsigned int level, const std::vector & selected_row_blocks = @@ -359,7 +364,7 @@ namespace MatrixFreeOperators /** * Get read access to the MatrixFree object stored with this operator. */ - std::shared_ptr> + std::shared_ptr> get_matrix_free() const; /** @@ -424,7 +429,8 @@ namespace MatrixFreeOperators /** * MatrixFree object to be used with this operator. */ - std::shared_ptr> data; + std::shared_ptr> + data; /** * A shared pointer to a diagonal matrix that stores the @@ -723,19 +729,23 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d = fe_degree + 1, int n_components = 1, - typename VectorType = LinearAlgebra::distributed::Vector> - class MassOperator : public Base + typename VectorType = LinearAlgebra::distributed::Vector, + typename VectorizedArrayType = + VectorizedArray> + class MassOperator : public Base { public: /** * Number alias. */ - using value_type = typename Base::value_type; + using value_type = + typename Base::value_type; /** * size_type needed for preconditioner classes. */ - using size_type = typename Base::size_type; + using size_type = + typename Base::size_type; /** * Constructor. @@ -763,10 +773,10 @@ namespace MatrixFreeOperators */ void local_apply_cell( - const MatrixFree & data, - VectorType & dst, - const VectorType & src, - const std::pair &cell_range) const; + const MatrixFree &data, + VectorType & dst, + const VectorType & src, + const std::pair & cell_range) const; }; @@ -787,19 +797,23 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d = fe_degree + 1, int n_components = 1, - typename VectorType = LinearAlgebra::distributed::Vector> - class LaplaceOperator : public Base + typename VectorType = LinearAlgebra::distributed::Vector, + typename VectorizedArrayType = + VectorizedArray> + class LaplaceOperator : public Base { public: /** * Number alias. */ - using value_type = typename Base::value_type; + using value_type = + typename Base::value_type; /** * size_type needed for preconditioner classes. */ - using size_type = typename Base::size_type; + using size_type = + typename Base::size_type; /** * Constructor. @@ -862,8 +876,8 @@ namespace MatrixFreeOperators * will delete the table. */ void - set_coefficient(const std::shared_ptr>> - &scalar_coefficient); + set_coefficient( + const std::shared_ptr> &scalar_coefficient); virtual void clear(); @@ -874,7 +888,7 @@ namespace MatrixFreeOperators * The function will throw an error if coefficients are not previously set * by set_coefficient() function. */ - std::shared_ptr>> + std::shared_ptr> get_coefficient(); private: @@ -891,18 +905,18 @@ namespace MatrixFreeOperators */ void local_apply_cell( - const MatrixFree & data, - VectorType & dst, - const VectorType & src, - const std::pair &cell_range) const; + const MatrixFree &data, + VectorType & dst, + const VectorType & src, + const std::pair & cell_range) const; /** * Apply diagonal part of the Laplace operator on a cell. */ void local_diagonal_cell( - const MatrixFree &data, - VectorType & dst, + const MatrixFree &data, + VectorType & dst, const VectorType &, const std::pair &cell_range) const; @@ -918,7 +932,7 @@ namespace MatrixFreeOperators /** * User-provided heterogeneity coefficient. */ - std::shared_ptr>> scalar_coefficient; + std::shared_ptr> scalar_coefficient; }; @@ -1055,20 +1069,21 @@ namespace MatrixFreeOperators //----------------- Base operator ----------------------------- - template - Base::Base() + template + Base::Base() : Subscriptor() , have_interface_matrices(false) {} - template - typename Base::size_type - Base::m() const + template + typename Base::size_type + Base::m() const { Assert(data.get() != nullptr, ExcNotInitialized()); - typename Base::size_type total_size = 0; + typename Base::size_type total_size = + 0; for (unsigned int i = 0; i < selected_rows.size(); ++i) total_size += data->get_vector_partitioner(selected_rows[i])->size(); return total_size; @@ -1076,12 +1091,13 @@ namespace MatrixFreeOperators - template - typename Base::size_type - Base::n() const + template + typename Base::size_type + Base::n() const { Assert(data.get() != nullptr, ExcNotInitialized()); - typename Base::size_type total_size = 0; + typename Base::size_type total_size = + 0; for (unsigned int i = 0; i < selected_columns.size(); ++i) total_size += data->get_vector_partitioner(selected_columns[i])->size(); return total_size; @@ -1089,9 +1105,9 @@ namespace MatrixFreeOperators - template + template void - Base::clear() + Base::clear() { data.reset(); inverse_diagonal_entries.reset(); @@ -1099,10 +1115,10 @@ namespace MatrixFreeOperators - template - typename Base::value_type - Base::el(const unsigned int row, - const unsigned int col) const + template + typename Base::value_type + Base::el(const unsigned int row, + const unsigned int col) const { (void)col; Assert(row == col, ExcNotImplemented()); @@ -1114,9 +1130,10 @@ namespace MatrixFreeOperators - template + template void - Base::initialize_dof_vector(VectorType &vec) const + Base::initialize_dof_vector( + VectorType &vec) const { Assert(data.get() != nullptr, ExcNotInitialized()); AssertDimension(BlockHelper::n_blocks(vec), selected_rows.size()); @@ -1138,11 +1155,11 @@ namespace MatrixFreeOperators - template + template void - Base::initialize( - std::shared_ptr< - const MatrixFree::value_type>> data_, + Base::initialize( + std::shared_ptr> + data_, const std::vector &given_row_selection, const std::vector &given_column_selection) { @@ -1191,11 +1208,11 @@ namespace MatrixFreeOperators - template + template void - Base::initialize( - std::shared_ptr< - const MatrixFree::value_type>> data_, + Base::initialize( + std::shared_ptr> + data_, const MGConstrainedDoFs & mg_constrained_dofs, const unsigned int level, const std::vector &given_row_selection) @@ -1207,13 +1224,14 @@ namespace MatrixFreeOperators - template + template void - Base::initialize( - std::shared_ptr> data_, - const std::vector & mg_constrained_dofs, - const unsigned int level, - const std::vector & given_row_selection) + Base::initialize( + std::shared_ptr> + data_, + const std::vector &mg_constrained_dofs, + const unsigned int level, + const std::vector & given_row_selection) { AssertThrow(level != numbers::invalid_unsigned_int, ExcMessage("level is not set")); @@ -1276,9 +1294,10 @@ namespace MatrixFreeOperators - template + template void - Base::set_constrained_entries_to_one(VectorType &dst) const + Base::set_constrained_entries_to_one( + VectorType &dst) const { for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j) { @@ -1294,43 +1313,49 @@ namespace MatrixFreeOperators - template + template void - Base::vmult(VectorType &dst, const VectorType &src) const + Base::vmult(VectorType & dst, + const VectorType &src) const { - using Number = typename Base::value_type; - dst = Number(0.); + using Number = + typename Base::value_type; + dst = Number(0.); vmult_add(dst, src); } - template + template void - Base::vmult_add(VectorType &dst, const VectorType &src) const + Base::vmult_add( + VectorType & dst, + const VectorType &src) const { mult_add(dst, src, false); } - template + template void - Base::Tvmult_add(VectorType & dst, - const VectorType &src) const + Base::Tvmult_add( + VectorType & dst, + const VectorType &src) const { mult_add(dst, src, true); } - template + template void - Base::adjust_ghost_range_if_necessary( + Base::adjust_ghost_range_if_necessary( const VectorType &src, const bool is_row) const { - using Number = typename Base::value_type; + using Number = + typename Base::value_type; for (unsigned int i = 0; i < BlockHelper::n_blocks(src); ++i) { const unsigned int mf_component = @@ -1361,12 +1386,14 @@ namespace MatrixFreeOperators - template + template void - Base::preprocess_constraints(VectorType & dst, - const VectorType &src) const + Base::preprocess_constraints( + VectorType & dst, + const VectorType &src) const { - using Number = typename Base::value_type; + using Number = + typename Base::value_type; adjust_ghost_range_if_necessary(src, false); adjust_ghost_range_if_necessary(dst, true); @@ -1389,11 +1416,12 @@ namespace MatrixFreeOperators - template + template void - Base::mult_add(VectorType & dst, - const VectorType &src, - const bool transpose) const + Base::mult_add( + VectorType & dst, + const VectorType &src, + const bool transpose) const { AssertDimension(dst.size(), src.size()); AssertDimension(BlockHelper::n_blocks(dst), BlockHelper::n_blocks(src)); @@ -1408,10 +1436,11 @@ namespace MatrixFreeOperators - template + template void - Base::postprocess_constraints(VectorType & dst, - const VectorType &src) const + Base::postprocess_constraints( + VectorType & dst, + const VectorType &src) const { for (unsigned int j = 0; j < BlockHelper::n_blocks(dst); ++j) { @@ -1441,12 +1470,14 @@ namespace MatrixFreeOperators - template + template void - Base::vmult_interface_down(VectorType & dst, - const VectorType &src) const + Base::vmult_interface_down( + VectorType & dst, + const VectorType &src) const { - using Number = typename Base::value_type; + using Number = + typename Base::value_type; AssertDimension(dst.size(), src.size()); adjust_ghost_range_if_necessary(src, false); adjust_ghost_range_if_necessary(dst, true); @@ -1493,12 +1524,14 @@ namespace MatrixFreeOperators - template + template void - Base::vmult_interface_up(VectorType & dst, - const VectorType &src) const + Base::vmult_interface_up( + VectorType & dst, + const VectorType &src) const { - using Number = typename Base::value_type; + using Number = + typename Base::value_type; AssertDimension(dst.size(), src.size()); adjust_ghost_range_if_necessary(src, false); adjust_ghost_range_if_necessary(dst, true); @@ -1532,20 +1565,23 @@ namespace MatrixFreeOperators - template + template void - Base::Tvmult(VectorType &dst, const VectorType &src) const + Base::Tvmult( + VectorType & dst, + const VectorType &src) const { - using Number = typename Base::value_type; - dst = Number(0.); + using Number = + typename Base::value_type; + dst = Number(0.); Tvmult_add(dst, src); } - template + template std::size_t - Base::memory_consumption() const + Base::memory_consumption() const { return inverse_diagonal_entries.get() != nullptr ? inverse_diagonal_entries->memory_consumption() : @@ -1554,19 +1590,22 @@ namespace MatrixFreeOperators - template - std::shared_ptr< - const MatrixFree::value_type>> - Base::get_matrix_free() const + template + std::shared_ptr::value_type, + VectorizedArrayType>> + Base::get_matrix_free() const { return data; } - template + template const std::shared_ptr> & - Base::get_matrix_diagonal_inverse() const + Base::get_matrix_diagonal_inverse() + const { Assert(inverse_diagonal_entries.get() != nullptr && inverse_diagonal_entries->m() > 0, @@ -1576,9 +1615,9 @@ namespace MatrixFreeOperators - template + template const std::shared_ptr> & - Base::get_matrix_diagonal() const + Base::get_matrix_diagonal() const { Assert(diagonal_entries.get() != nullptr && diagonal_entries->m() > 0, ExcNotInitialized()); @@ -1587,22 +1626,24 @@ namespace MatrixFreeOperators - template + template void - Base::Tapply_add(VectorType & dst, - const VectorType &src) const + Base::Tapply_add( + VectorType & dst, + const VectorType &src) const { apply_add(dst, src); } - template + template void - Base::precondition_Jacobi( - VectorType & dst, - const VectorType & src, - const typename Base::value_type omega) const + Base::precondition_Jacobi( + VectorType & dst, + const VectorType & src, + const typename Base::value_type omega) + const { Assert(inverse_diagonal_entries.get() && inverse_diagonal_entries->m() > 0, ExcNotInitialized()); @@ -1699,10 +1740,15 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> - MassOperator:: - MassOperator() - : Base() + typename VectorType, + typename VectorizedArrayType> + MassOperator::MassOperator() + : Base() {} @@ -1711,13 +1757,20 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> + typename VectorType, + typename VectorizedArrayType> void - MassOperator:: - compute_diagonal() + MassOperator::compute_diagonal() { - using Number = typename Base::value_type; - Assert((Base::data.get() != nullptr), ExcNotInitialized()); + using Number = + typename Base::value_type; + Assert((Base::data.get() != nullptr), + ExcNotInitialized()); this->inverse_diagonal_entries.reset(new DiagonalMatrix()); this->diagonal_entries.reset(new DiagonalMatrix()); @@ -1747,15 +1800,19 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> + typename VectorType, + typename VectorizedArrayType> void - MassOperator:: - apply_add(VectorType &dst, const VectorType &src) const + MassOperator::apply_add(VectorType & dst, + const VectorType &src) const { - Base::data->cell_loop(&MassOperator::local_apply_cell, - this, - dst, - src); + Base::data->cell_loop( + &MassOperator::local_apply_cell, this, dst, src); } @@ -1764,18 +1821,33 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> + typename VectorType, + typename VectorizedArrayType> void - MassOperator:: + MassOperator:: local_apply_cell( - const MatrixFree::value_type> &data, - VectorType & dst, - const VectorType & src, + const MatrixFree< + dim, + typename Base::value_type, + VectorizedArrayType> & data, + VectorType & dst, + const VectorType & src, const std::pair &cell_range) const { - using Number = typename Base::value_type; - FEEvaluation phi( - data, this->selected_rows[0]); + using Number = + typename Base::value_type; + FEEvaluation + phi(data, this->selected_rows[0]); for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) { phi.reinit(cell); @@ -1795,10 +1867,15 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> - LaplaceOperator:: - LaplaceOperator() - : Base() + typename VectorType, + typename VectorizedArrayType> + LaplaceOperator::LaplaceOperator() + : Base() {} @@ -1807,12 +1884,17 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> + typename VectorType, + typename VectorizedArrayType> void - LaplaceOperator:: - clear() + LaplaceOperator::clear() { - Base::clear(); + Base::clear(); scalar_coefficient.reset(); } @@ -1822,13 +1904,17 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> + typename VectorType, + typename VectorizedArrayType> void - LaplaceOperator:: + LaplaceOperator:: set_coefficient( - const std::shared_ptr< - Table<2, VectorizedArray::value_type>>> - &scalar_coefficient_) + const std::shared_ptr> &scalar_coefficient_) { scalar_coefficient = scalar_coefficient_; } @@ -1839,16 +1925,15 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> - std::shared_ptr< - Table<2, - VectorizedArray::value_type>>> - LaplaceOperator:: - get_coefficient() + typename VectorType, + typename VectorizedArrayType> + std::shared_ptr> + LaplaceOperator::get_coefficient() { Assert(scalar_coefficient.get(), ExcNotInitialized()); return scalar_coefficient; @@ -1860,13 +1945,20 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> + typename VectorType, + typename VectorizedArrayType> void - LaplaceOperator:: - compute_diagonal() + LaplaceOperator::compute_diagonal() { - using Number = typename Base::value_type; - Assert((Base::data.get() != nullptr), ExcNotInitialized()); + using Number = + typename Base::value_type; + Assert((Base::data.get() != nullptr), + ExcNotInitialized()); this->inverse_diagonal_entries.reset(new DiagonalMatrix()); this->diagonal_entries.reset(new DiagonalMatrix()); @@ -1902,15 +1994,19 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> + typename VectorType, + typename VectorizedArrayType> void - LaplaceOperator:: - apply_add(VectorType &dst, const VectorType &src) const + LaplaceOperator::apply_add(VectorType & dst, + const VectorType &src) const { - Base::data->cell_loop(&LaplaceOperator::local_apply_cell, - this, - dst, - src); + Base::data->cell_loop( + &LaplaceOperator::local_apply_cell, this, dst, src); } namespace Implementation @@ -1933,16 +2029,23 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> + typename VectorType, + typename VectorizedArrayType> void - LaplaceOperator:: + LaplaceOperator:: do_operation_on_cell( - FEEvaluation::value_type> &phi, - const unsigned int cell) const + FEEvaluation< + dim, + fe_degree, + n_q_points_1d, + n_components, + typename Base::value_type> &phi, + const unsigned int cell) const { phi.evaluate(false, true, false); if (scalar_coefficient.get()) @@ -1972,16 +2075,26 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> + typename VectorType, + typename VectorizedArrayType> void - LaplaceOperator:: + LaplaceOperator:: local_apply_cell( - const MatrixFree::value_type> &data, - VectorType & dst, - const VectorType & src, + const MatrixFree< + dim, + typename Base::value_type, + VectorizedArrayType> & data, + VectorType & dst, + const VectorType & src, const std::pair &cell_range) const { - using Number = typename Base::value_type; + using Number = + typename Base::value_type; FEEvaluation phi( data, this->selected_rows[0]); for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) @@ -1998,30 +2111,37 @@ namespace MatrixFreeOperators int fe_degree, int n_q_points_1d, int n_components, - typename VectorType> + typename VectorType, + typename VectorizedArrayType> void - LaplaceOperator:: + LaplaceOperator:: local_diagonal_cell( - const MatrixFree::value_type> &data, - VectorType & dst, + const MatrixFree< + dim, + typename Base::value_type, + VectorizedArrayType> &data, + VectorType & dst, const VectorType &, const std::pair &cell_range) const { - using Number = typename Base::value_type; - using vectorized_value_type = - typename MatrixFree::value_type>:: - vectorized_value_type; + using Number = + typename Base::value_type; FEEvaluation phi( data, this->selected_rows[0]); for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) { phi.reinit(cell); - vectorized_value_type local_diagonal_vector[phi.static_dofs_per_cell]; + VectorizedArrayType local_diagonal_vector[phi.static_dofs_per_cell]; for (unsigned int i = 0; i < phi.dofs_per_component; ++i) { for (unsigned int j = 0; j < phi.dofs_per_component; ++j) - phi.begin_dof_values()[j] = vectorized_value_type(); + phi.begin_dof_values()[j] = VectorizedArrayType(); phi.begin_dof_values()[i] = 1.; do_operation_on_cell(phi, cell); local_diagonal_vector[i] = phi.begin_dof_values()[i]; -- 2.39.5