From bba2c791a07ccfee4e0867d196de9cdaf9ce47e4 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 19 Feb 2014 09:56:31 +0000 Subject: [PATCH] Make it possible to allocate FEEvaluation on the heap using an AlignedVector type. It turns out that the revision 32498 must be reverted again: a shared ptr is needed to keep track of the allocated shape info in case it is copied around git-svn-id: https://svn.dealii.org/trunk@32511 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/base/aligned_vector.h | 62 ++-- .../deal.II/matrix_free/fe_evaluation.h | 239 ++++++++++-- tests/matrix_free/copy_feevaluation.cc | 349 ++++++++++++++++++ tests/matrix_free/copy_feevaluation.output | 12 + 4 files changed, 592 insertions(+), 70 deletions(-) create mode 100644 tests/matrix_free/copy_feevaluation.cc create mode 100644 tests/matrix_free/copy_feevaluation.output diff --git a/deal.II/include/deal.II/base/aligned_vector.h b/deal.II/include/deal.II/base/aligned_vector.h index c013786d69..2916bcd0ba 100644 --- a/deal.II/include/deal.II/base/aligned_vector.h +++ b/deal.II/include/deal.II/base/aligned_vector.h @@ -23,7 +23,6 @@ #include #include #include -#include #if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 @@ -76,7 +75,8 @@ public: * Sets the vector size to the given size and initializes all elements with * T(). */ - AlignedVector (const size_type size); + AlignedVector (const size_type size, + const T &init = T()); /** * Destructor. @@ -259,13 +259,10 @@ namespace internal static const std::size_t minimum_parallel_grain_size = 160000/sizeof(T)+1; public: /** - * Constructor. Issues a parallel call if - * there are sufficiently many elements, - * otherwise work in serial. Copies the data - * from source to destination and then calls - * destructor on the source. If the optional - * argument is set to true, the source is left - * untouched instead. + * Constructor. Issues a parallel call if there are sufficiently many + * elements, otherwise work in serial. Copies the data from source to + * destination and then calls destructor on the source. If the optional + * argument is set to true, the source is left untouched instead. */ AlignedVectorMove (T *source_begin, T *source_end, @@ -285,30 +282,26 @@ namespace internal } /** - * This method moves elements from the source - * to the destination given in the constructor - * on a subrange given by two integers. + * This method moves elements from the source to the destination given in + * the constructor on a subrange given by two integers. */ virtual void apply_to_subrange (const std::size_t begin, const std::size_t end) const { - // for classes trivial assignment can use - // memcpy + // for classes trivial assignment can use memcpy if (std_cxx1x::is_trivial::value == true) std::memcpy (destination_+begin, source_+begin, (end-begin)*sizeof(T)); else if (copy_only_ == false) for (std::size_t i=begin; i::value == true && trivial_element) std::memset (destination_+begin, 0, (end-begin)*sizeof(T)); else + // initialize memory and set for (std::size_t i=begin; i::AlignedVector () template < class T > inline -AlignedVector::AlignedVector (const size_type size) +AlignedVector::AlignedVector (const size_type size, + const T &init) : _data (0), _end_data (0), _end_allocated (0) { if (size > 0) - resize (size); + resize (size, init); } @@ -504,11 +492,9 @@ AlignedVector::reserve (const size_type size_alloc) #if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 - // allocate and align along boundaries of the size of - // VectorizedArray, which is 16 bytes for SSE and 32 bytes for - // AVX - T *new_data = static_cast(_mm_malloc (size_actual_allocate, - sizeof(VectorizedArray))); + // allocate and align along 64-byte boundaries (this is enough for all + // levels of vectorization currently supported by deal.II) + T *new_data = static_cast(_mm_malloc (size_actual_allocate, 64)); #else T *new_data = static_cast(malloc (size_actual_allocate)); #endif diff --git a/deal.II/include/deal.II/matrix_free/fe_evaluation.h b/deal.II/include/deal.II/matrix_free/fe_evaluation.h index 92af89a2f7..141ebe579e 100644 --- a/deal.II/include/deal.II/matrix_free/fe_evaluation.h +++ b/deal.II/include/deal.II/matrix_free/fe_evaluation.h @@ -88,11 +88,6 @@ public: * @name 1: General operations */ //@{ - /** - * Destructor. - */ - ~FEEvaluationBase(); - /** * Initializes the operation pointer to the current cell. Unlike the * FEValues::reinit function, where the information related to a particular @@ -535,6 +530,11 @@ protected: const DoFHandler &dof_handler, const unsigned int base_element = 0); + /** + * Copy constructor + */ + FEEvaluationBase (const FEEvaluationBase &other); + /** * A unified function to read from and write into vectors based on the given * template operation. It can perform the operation for @p read_dof_values, @@ -641,10 +641,10 @@ protected: * Stores a reference to the unit cell data, i.e., values, gradients and * Hessians in 1D at the quadrature points that constitute the tensor * product. Also contained in matrix_info, but it simplifies code if we - * store a reference to it. If initialized without MatrixFree object, this - * call actually initializes the evaluation. + * store a reference to it. If the object is initialized without MatrixFree + * object, the constructor creates this data structure. */ - const internal::MatrixFreeFunctions::ShapeInfo *data; + std_cxx1x::shared_ptr > data; /** * A pointer to the Cartesian Jacobian information of the present cell. Only @@ -797,7 +797,6 @@ protected: const unsigned int dofs_per_cell, const unsigned int n_q_points); - /** * Constructor that comes with reduced functionality and works similar as * FEValues. The user has to provide a structure of type EvaluatedGeometry @@ -817,6 +816,11 @@ protected: FEEvaluationAccess (const EvaluatedGeometry &geometry, const DoFHandler &dof_handler, const unsigned int base_element = 0); + + /** + * Copy constructor + */ + FEEvaluationAccess (const FEEvaluationAccess &other); }; @@ -959,6 +963,11 @@ protected: FEEvaluationAccess (const EvaluatedGeometry &geometry, const DoFHandler &dof_handler, const unsigned int base_element = 0); + + /** + * Copy constructor + */ + FEEvaluationAccess (const FEEvaluationAccess &other); }; @@ -1110,6 +1119,11 @@ protected: FEEvaluationAccess (const EvaluatedGeometry &geometry, const DoFHandler &dof_handler, const unsigned int base_element = 0); + + /** + * Copy constructor + */ + FEEvaluationAccess (const FEEvaluationAccess &other); }; @@ -1197,6 +1211,11 @@ public: const DoFHandler &dof_handler, const unsigned int base_element = 0); + /** + * Copy constructor + */ + FEEvaluationGeneral (const FEEvaluationGeneral &other); + /** * Evaluates the function values, the gradients, and the Laplacians of the * FE function given at the DoF values in the input vector at the quadrature @@ -1267,6 +1286,12 @@ protected: * Internally stored variables for the different data fields. */ VectorizedArray my_data_array[n_components*(dofs_per_cell+(dim*dim+2*dim+1)*n_q_points)]; + +private: + /** + * Sets the pointers from the data array to values_dof, etc. + */ + void set_data_pointers(); }; @@ -1347,6 +1372,11 @@ public: const unsigned int fe_no = 0, const unsigned int quad_no = 0); + /** + * Copy constructor + */ + FEEvaluation (const FEEvaluation &other); + /** * Constructor that comes with reduced functionality and works similar as * FEValues. The user has to provide a structure of type EvaluatedGeometry @@ -1521,6 +1551,11 @@ public: const DoFHandler &dof_handler, const unsigned int base_element = 0); + /** + * Copy constructor + */ + FEEvaluationGL (const FEEvaluationGL &other); + /** * Evaluates the function values, the gradients, and the Hessians of the FE * function given at the DoF values in the input vector at the quadrature @@ -1643,6 +1678,11 @@ public: const DoFHandler &dof_handler, const unsigned int base_element = 0); + /** + * Copy constructor + */ + FEEvaluationDGP (const FEEvaluationDGP &other); + /** * Evaluates the function values, the gradients, and the Hessians of the FE * function given at the DoF values in the input vector at the quadrature @@ -1672,6 +1712,34 @@ public: #ifndef DOXYGEN +namespace internal +{ + namespace MatrixFreeFunctions + { + // a small class that gives control over the delete behavior of + // std::shared_ptr: we need to disable it when we initialize a pointer + // from another structure. + template + struct DummyDeleter + { + DummyDeleter (const bool do_delete = false) + : + do_delete(do_delete) + {} + + void operator () (CLASS *pointer) + { + if (do_delete) + delete pointer; + } + + const bool do_delete; + }; + } +} + + + /*----------------------- FEEvaluationBase ----------------------------------*/ template @@ -1695,7 +1763,9 @@ FEEvaluationBase mapping_info (&data_in.get_mapping_info()), data (&data_in.get_shape_info (fe_no_in, quad_no_in, active_fe_index, - active_quad_index)), + active_quad_index), + internal::MatrixFreeFunctions::DummyDeleter + >(false)), cartesian_data (0), jacobian (0), J_value (0), @@ -1786,16 +1856,44 @@ FEEvaluationBase template inline -FEEvaluationBase::~FEEvaluationBase() +FEEvaluationBase +::FEEvaluationBase (const FEEvaluationBase &other) + : + quad_no (other.quad_no), + n_fe_components (other.n_fe_components), + active_fe_index (other.active_fe_index), + active_quad_index (other.active_quad_index), + matrix_info (other.matrix_info), + dof_info (other.dof_info), + mapping_info (other.mapping_info), + data (other.data), + cartesian_data (other.cartesian_data), + jacobian (other.jacobian), + J_value (other.J_value), + quadrature_weights (other.quadrature_weights), + quadrature_points (other.quadrature_points), + jacobian_grad (other.jacobian_grad), + jacobian_grad_upper(other.jacobian_grad_upper), + cell (other.cell), + cell_type (other.cell_type), + cell_data_number (other.cell_data_number), + evaluated_geometry (other.evaluated_geometry), + dof_handler (other.dof_handler) { - // delete memory held by data in case this structure was initialized without - // a MatrixFree object which could hold the data structure. - if (matrix_info == 0) - delete data; + for (unsigned int c=0; c inline void @@ -3460,6 +3558,16 @@ FEEvaluationAccess +template +inline +FEEvaluationAccess +::FEEvaluationAccess (const FEEvaluationAccess &other) + : + FEEvaluationBase (other) +{} + + + /*-------------------- FEEvaluationAccess scalar ----------------------------*/ @@ -3491,6 +3599,16 @@ FEEvaluationAccess +template +inline +FEEvaluationAccess +::FEEvaluationAccess (const FEEvaluationAccess&other) + : + FEEvaluationBase (other) +{} + + + template inline VectorizedArray @@ -3715,6 +3833,16 @@ FEEvaluationAccess +template +inline +FEEvaluationAccess +::FEEvaluationAccess (const FEEvaluationAccess&other) + : + FEEvaluationBase (other) +{} + + + template inline Tensor<2,dim,VectorizedArray > @@ -4043,20 +4171,8 @@ FEEvaluationGeneral : BaseClass (data_in, fe_no, quad_no, dofs_per_cell, n_q_points) { - // set the pointers to the correct position in the data array - for (unsigned int c=0; cvalues_dofs[c] = &my_data_array[c*dofs_per_cell]; - this->values_quad[c] = &my_data_array[n_components*dofs_per_cell+c*n_q_points]; - for (unsigned int d=0; dgradients_quad[c][d] = &my_data_array[n_components*(dofs_per_cell+n_q_points) - + - (c*dim+d)*n_q_points]; - for (unsigned int d=0; d<(dim*dim+dim)/2; ++d) - this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+dofs_per_cell) - + - (c*(dim*dim+dim)+d)*n_q_points]; - } + set_data_pointers(); + #ifdef DEBUG // print error message when the dimensions do not match. Propose a possible // fix @@ -4173,6 +4289,30 @@ FEEvaluationGeneral const unsigned int base_element) : BaseClass (geometry, dof_handler, base_element) +{ + set_data_pointers(); +} + + +template +inline +FEEvaluationGeneral +::FEEvaluationGeneral (const FEEvaluationGeneral &other) + : + BaseClass (other) +{ + set_data_pointers(); +} + + + +template +inline +void +FEEvaluationGeneral +::set_data_pointers() { // set the pointers to the correct position in the data array for (unsigned int c=0; c +template +inline +FEEvaluation +::FEEvaluation (const FEEvaluation &other) + : + BaseClass (other) +{ + compute_even_odd_factors(); +} + + + template inline @@ -5957,7 +6110,7 @@ FEEvaluation if (fe_degree > 1 || n_q_points_1d > 3) internal::apply_tensor_product_evenodd, direction, dof_to_quad, add, 0> - (shape_val_evenodd, in, out); + (&shape_val_evenodd[0], in, out); else internal::apply_tensor_product_values, direction, dof_to_quad, add> @@ -5978,7 +6131,7 @@ FEEvaluation if (fe_degree > 1 || n_q_points_1d > 3) internal::apply_tensor_product_evenodd, direction, dof_to_quad, add, 1> - (shape_gra_evenodd, in, out); + (&shape_gra_evenodd[0], in, out); else internal::apply_tensor_product_gradients, direction, dof_to_quad, add> @@ -6002,7 +6155,7 @@ FEEvaluation if (fe_degree > 1 || n_q_points_1d > 3) internal::apply_tensor_product_evenodd, direction, dof_to_quad, add, 2> - (shape_hes_evenodd, in, out); + (&shape_hes_evenodd[0], in, out); else internal::apply_tensor_product_hessians, direction, dof_to_quad, add> @@ -6069,6 +6222,16 @@ FEEvaluationGL +template +inline +FEEvaluationGL +::FEEvaluationGL (const FEEvaluationGL &other) + : + BaseClass (other) +{} + + + template inline void @@ -6293,6 +6456,18 @@ FEEvaluationDGP {} + +template +inline +FEEvaluationDGP +::FEEvaluationDGP (const FEEvaluationDGP &other) + : + BaseClass (other) +{} + + + template inline diff --git a/tests/matrix_free/copy_feevaluation.cc b/tests/matrix_free/copy_feevaluation.cc new file mode 100644 index 0000000000..cb0d126b22 --- /dev/null +++ b/tests/matrix_free/copy_feevaluation.cc @@ -0,0 +1,349 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// same test as matrix_vector_stokes_noflux, but allocating FEEvaluation on +// the heap (using AlignedVector) instead of allocating it on the stack. Tests +// also copy constructors of FEEvaluation. + +#include "../tests.h" + +std::ofstream logfile("output"); + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + + +template +class MatrixFreeTest +{ +public: + typedef typename DoFHandler::active_cell_iterator CellIterator; + typedef double Number; + + MatrixFreeTest(const MatrixFree &data_in): + data (data_in) + {}; + + void + local_apply (const MatrixFree &data, + VectorType &dst, + const VectorType &src, + const std::pair &cell_range) const + { + typedef VectorizedArray vector_t; + // allocate FEEvaluation. This test will test proper alignment + AlignedVector > velocity + (1, FEEvaluation(data, 0)); + AlignedVector > pressure + (1, FEEvaluation(data, 1)); + + for (unsigned int cell=cell_range.first; cell sym_grad_u = + velocity[0].get_symmetric_gradient (q); + vector_t pres = pressure[0].get_value(q); + vector_t div = -velocity[0].get_divergence(q); + pressure[0].submit_value (div, q); + + // subtract p * I + for (unsigned int d=0; d::local_apply, + this, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void test () +{ + Triangulation triangulation; + GridGenerator::hyper_shell (triangulation, Point(), + 0.5, 1., 96, true); + static HyperShellBoundary boundary; + triangulation.set_boundary (0, boundary); + triangulation.set_boundary (1, boundary); + triangulation.begin_active()->set_refine_flag(); + triangulation.last()->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + triangulation.refine_global (3-dim); + triangulation.last()->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + + MappingQ mapping (3); + FE_Q fe_u_scal (fe_degree+1); + FESystem fe_u (fe_u_scal,dim); + FE_Q fe_p (fe_degree); + FESystem fe (fe_u_scal, dim, fe_p, 1); + DoFHandler dof_handler_u (triangulation); + DoFHandler dof_handler_p (triangulation); + DoFHandler dof_handler (triangulation); + + MatrixFree mf_data; + + ConstraintMatrix constraints, constraints_u, constraints_p; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + BlockVector mf_solution; + + dof_handler.distribute_dofs (fe); + dof_handler_u.distribute_dofs (fe_u); + dof_handler_p.distribute_dofs (fe_p); + std::vector stokes_sub_blocks (dim+1,0); + stokes_sub_blocks[dim] = 1; + DoFRenumbering::component_wise (dof_handler, stokes_sub_blocks); + + std::set no_normal_flux_boundaries; + no_normal_flux_boundaries.insert (0); + no_normal_flux_boundaries.insert (1); + DoFTools::make_hanging_node_constraints (dof_handler, + constraints); + VectorTools::compute_no_normal_flux_constraints (dof_handler, 0, + no_normal_flux_boundaries, + constraints, mapping); + constraints.close (); + DoFTools::make_hanging_node_constraints (dof_handler_u, + constraints_u); + VectorTools::compute_no_normal_flux_constraints (dof_handler_u, 0, + no_normal_flux_boundaries, + constraints_u, mapping); + constraints_u.close (); + DoFTools::make_hanging_node_constraints (dof_handler_p, + constraints_p); + constraints_p.close (); + + std::vector dofs_per_block (2); + DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, + stokes_sub_blocks); + + //std::cout << "Number of active cells: " + // << triangulation.n_active_cells() + // << std::endl + // << "Number of degrees of freedom: " + // << dof_handler.n_dofs() + // << " (" << n_u << '+' << n_p << ')' + // << std::endl; + + { + BlockCompressedSimpleSparsityPattern csp (2,2); + + for (unsigned int d=0; d<2; ++d) + for (unsigned int e=0; e<2; ++e) + csp.block(d,e).reinit (dofs_per_block[d], dofs_per_block[e]); + + csp.collect_sizes(); + + DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false); + sparsity_pattern.copy_from (csp); + } + + system_matrix.reinit (sparsity_pattern); + + // this is from step-22 + { + QGauss quadrature_formula(fe_degree+2); + + FEValues fe_values (mapping, fe, quadrature_formula, + update_values | + update_JxW_values | + update_gradients); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + std::vector > phi_grads_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + local_matrix = 0; + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (local_matrix, + local_dof_indices, + system_matrix); + } + } + + + solution.reinit (2); + for (unsigned int d=0; d<2; ++d) + solution.block(d).reinit (dofs_per_block[d]); + solution.collect_sizes (); + + system_rhs.reinit (solution); + mf_solution.reinit (solution); + + // fill system_rhs with random numbers + for (unsigned int j=0; j*> dofs; + dofs.push_back(&dof_handler_u); + dofs.push_back(&dof_handler_p); + std::vector constraints; + constraints.push_back (&constraints_u); + constraints.push_back (&constraints_p); + QGauss<1> quad(fe_degree+2); + // no parallelism + mf_data.reinit (mapping, dofs, constraints, quad, + typename MatrixFree::AdditionalData + (MPI_COMM_WORLD, + MatrixFree::AdditionalData::none)); + } + + system_matrix.vmult (solution, system_rhs); + + MatrixFreeTest > mf (mf_data); + mf.vmult (mf_solution, system_rhs); + + // Verification + mf_solution -= solution; + const double error = mf_solution.linfty_norm(); + const double relative = solution.linfty_norm(); + deallog << "Verification fe degree " << fe_degree << ": " + << error/relative << std::endl << std::endl; +} + + + +int main () +{ + deallog.attach(logfile); + deallog.depth_console(0); + + deallog << std::setprecision (3); + + { + deallog << std::endl << "Test with doubles" << std::endl << std::endl; + deallog.threshold_double(1.e-12); + deallog.push("2d"); + test<2,1>(); + test<2,2>(); + test<2,3>(); + deallog.pop(); + deallog.push("3d"); + test<3,1>(); + deallog.pop(); + } +} + diff --git a/tests/matrix_free/copy_feevaluation.output b/tests/matrix_free/copy_feevaluation.output new file mode 100644 index 0000000000..0434059632 --- /dev/null +++ b/tests/matrix_free/copy_feevaluation.output @@ -0,0 +1,12 @@ + +DEAL:: +DEAL::Test with doubles +DEAL:: +DEAL:2d::Verification fe degree 1: 0 +DEAL:2d:: +DEAL:2d::Verification fe degree 2: 0 +DEAL:2d:: +DEAL:2d::Verification fe degree 3: 0 +DEAL:2d:: +DEAL:3d::Verification fe degree 1: 0 +DEAL:3d:: -- 2.39.5