From 23f11135d1b1c5f2f79185aece523e2544e5b2ad Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sat, 17 Nov 2012 21:03:29 +0000 Subject: [PATCH] Skip hashes for comparing constraints in favor of a std::map with user-defined comparator. This makes code cleaner and runs faster since boost's hash value is slow. Fix error in dof_info.templates.h when hp elements are used but only one element is present. Fix error in AVX detection. git-svn-id: https://svn.dealii.org/trunk@27557 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/aclocal.m4 | 5 +- deal.II/configure | 19 +- .../deal.II/matrix_free/dof_info.templates.h | 203 ++++++------------ .../deal.II/matrix_free/helper_functions.h | 76 +++---- .../matrix_free/mapping_info.templates.h | 160 ++++++-------- .../include/deal.II/matrix_free/matrix_free.h | 16 ++ .../matrix_free/matrix_free.templates.h | 140 +++++++----- 7 files changed, 267 insertions(+), 352 deletions(-) diff --git a/deal.II/aclocal.m4 b/deal.II/aclocal.m4 index 9fb7765dff..85da1c5535 100644 --- a/deal.II/aclocal.m4 +++ b/deal.II/aclocal.m4 @@ -2210,7 +2210,7 @@ AC_DEFUN(DEAL_II_DETECT_VECTORIZATION_LEVEL, dnl [ AC_LANG(C++) CXXFLAGS="$CXXFLAGSG" - dnl SSE2 check in debug mode + dnl SSE2 check AC_MSG_CHECKING(whether CPU supports SSE2) AC_TRY_RUN( [ @@ -2243,7 +2243,7 @@ AC_DEFUN(DEAL_II_DETECT_VECTORIZATION_LEVEL, dnl ], [ AC_MSG_RESULT(yes) - dnl AVX check in debug mode + dnl AVX check AC_MSG_CHECKING(whether CPU supports AVX) AC_TRY_RUN( [ @@ -2271,6 +2271,7 @@ AC_DEFUN(DEAL_II_DETECT_VECTORIZATION_LEVEL, dnl if (ptr[i] != 5.0625) return_value = 1; _mm_free (data); + return return_value; } ], [ diff --git a/deal.II/configure b/deal.II/configure index 391072862f..5d843cb4e1 100755 --- a/deal.II/configure +++ b/deal.II/configure @@ -1,5 +1,5 @@ #! /bin/sh -# From configure.in Revision: 27122 . +# From configure.in Revision: 27277 . # Guess values for system-dependent variables and create Makefiles. # Generated by GNU Autoconf 2.69 for deal.II 7.3.pre. # @@ -814,7 +814,6 @@ enable_mpi enable_threads enable_shared enable_parser -enable_mgcompatibility enable_compat_blocker with_boost with_petsc @@ -1475,10 +1474,6 @@ Optional Features: --enable-parser While switched on by default, this option allows to switch off support for the function parser in the contrib directory. - --enable-mgcompatibility - Use preconditioner interface in MGSmootherRelaxation - instead of the new interface using the function - step. Defaults to disabled. --enable-compat-blocker=mapping Block functions that implicitely assume a Q1 mapping @@ -9859,17 +9854,6 @@ fi -# Check whether --enable-mgcompatibility was given. -if test "${enable_mgcompatibility+set}" = set; then : - enableval=$enable_mgcompatibility; if test "x$enableval" = "xyes" ; then - { $as_echo "$as_me:${as_lineno-$LINENO}: result: enable multigrid compatibility mode" >&5 -$as_echo "enable multigrid compatibility mode" >&6; } - -$as_echo "#define DEAL_II_MULTIGRID_COMPATIBILITY 1" >>confdefs.h - - fi -fi - # Check whether --enable-compat-blocker was given. @@ -13583,6 +13567,7 @@ else if (ptr[i] != 5.0625) return_value = 1; _mm_free (data); + return return_value; } _ACEOF diff --git a/deal.II/include/deal.II/matrix_free/dof_info.templates.h b/deal.II/include/deal.II/matrix_free/dof_info.templates.h index 8c680c9646..6f2ee7ba7e 100644 --- a/deal.II/include/deal.II/matrix_free/dof_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/dof_info.templates.h @@ -21,8 +21,17 @@ DEAL_II_NAMESPACE_OPEN namespace internal { -namespace MatrixFreeFunctions -{ + namespace MatrixFreeFunctions + { + + struct ConstraintComparator + { + bool operator()(const std::pair &p1, + const std::pair &p2) const + { + return p1.second < p2.second; + } + }; /** * A struct that takes entries describing @@ -52,135 +61,62 @@ namespace MatrixFreeFunctions unsigned short insert_entries (const std::vector > &entries); - std::vector constraint_pool_data; - std::vector constraint_pool_row_index; - std::vector > pool_locations; - std::vector > constraint_entries; + std::vector > constraint_entries; std::vector constraint_indices; - std::vector one_constraint; - HashValue hashes; + std::pair,unsigned int> next_constraint; + std::map,unsigned int,FPArrayComparator > constraints; }; template ConstraintValues::ConstraintValues () : - hashes (1.) - { - constraint_pool_row_index.push_back (0); - } + constraints(FPArrayComparator(1.)) + {} template unsigned short ConstraintValues:: insert_entries (const std::vector > &entries) { - unsigned int insert_position = deal_II_numbers::invalid_unsigned_int; - - typedef std::vector >::iterator iter; - constraint_entries.resize(entries.size()); - one_constraint.resize(entries.size()); - constraint_indices.resize(entries.size()); - for (unsigned int j=0;j 0) { + constraint_indices.resize(entries.size()); + constraint_entries = entries; + std::sort(constraint_entries.begin(), constraint_entries.end(), + ConstraintComparator()); + for (unsigned int j=0;j test (hash_val, 0); - - // Try to find a constraint in the pool with - // the same hash value. - iter pos = std::lower_bound (pool_locations.begin(), - pool_locations.end(), - test); - - // If constraint has to be added, which will - // be its no. - test.second = constraint_pool_row_index.size()-1; - - // Hash value larger than all the ones - // before. We need to add it. - if (pos == pool_locations.end()) - goto insert; - - // A constraint in the pool with the same hash - // value identified. - else if (pos->first == test.first) - { - bool is_same = true; - while(is_same == true) - { - if(one_constraint.size()!= - (constraint_pool_row_index[pos->second+1]- - constraint_pool_row_index[pos->second])) - // The constraints have different length, and - // hence different. - is_same = false; - else - for (unsigned int q=0; qsecond]+q]- - one_constraint[q])>hashes.scaling) - { - is_same = false; - break; - } - if (is_same == false) - { - // Try if there is another constraint with the - // same hash value. - ++pos; - if (pos != pool_locations.end() && pos->first == test.first) - is_same = true; - else - goto insert; - } - else - { - // The constraint is the same as the - // (pos->second)th in the pool. Add the - // location of constraint in pool - insert_position = pos->second; - break; - } - } - } + // in pool. the initial implementation + // computed a hash value based on the truncated + // array (to given accuracy around 1e-13) in + // order to easily detect different arrays and + // then made a fine-grained check when the + // hash values were equal. this was quite + // lenghty and now we use a std::map with a + // user-defined comparator to compare floating + // point arrays to a tolerance 1e-13. + std::pair, unsigned int, + FPArrayComparator >::iterator, + bool> it = constraints.insert(next_constraint); + unsigned int insert_position = deal_II_numbers::invalid_unsigned_int; + if (it.second == false) + insert_position = it.first->second; else - { - // A new constraint has been identified. It - // needs to be added (at the position - // test->second in pool_locations). - insert: - pool_locations.insert (pos, test); - - // Remember hash value and location of - // constraint. - constraint_pool_data.insert (constraint_pool_data.end(), - one_constraint.begin(), - one_constraint.end()); - constraint_pool_row_index.push_back (constraint_pool_data.size()); - - // Add the location of constraint in pool. - insert_position = test.second; - } + insert_position = next_constraint.second; // we want to store the result as a short // variable, so we have to make sure that the @@ -214,6 +150,9 @@ namespace MatrixFreeFunctions dofs_per_cell (dof_info_in.dofs_per_cell), dofs_per_face (dof_info_in.dofs_per_face), store_plain_indices (dof_info_in.store_plain_indices), + cell_active_fe_index (dof_info_in.cell_active_fe_index), + max_fe_index (dof_info_in.max_fe_index), + fe_index_conversion (dof_info_in.fe_index_conversion), ghost_dofs (dof_info_in.ghost_dofs) {} @@ -229,9 +168,13 @@ namespace MatrixFreeFunctions ghost_dofs.clear(); dofs_per_cell.clear(); dofs_per_face.clear(); + n_components = 0; row_starts_plain_indices.clear(); plain_dof_indices.clear(); store_plain_indices = false; + cell_active_fe_index.clear(); + max_fe_index = 0; + fe_index_conversion.clear(); } @@ -251,7 +194,7 @@ namespace MatrixFreeFunctions const unsigned int n_owned = last_owned - first_owned; std::pair constraint_iterator (0,0); - unsigned int dofs_this_cell = (cell_active_fe_index.size() == 0) ? + unsigned int dofs_this_cell = (cell_active_fe_index.empty()) ? dofs_per_cell[0] : dofs_per_cell[cell_active_fe_index[cell_number]]; for (unsigned int i=0; i * blb = begin_indicators(cell_number); - for (unsigned int j=0; j (new_dof_indices.size(), new_constraint_indicator.size(), 0); - + AssertDimension(dof_indices.size(), new_dof_indices.size()); AssertDimension(constraint_indicator.size(), new_constraint_indicator.size()); @@ -984,7 +912,8 @@ namespace MatrixFreeFunctions for(counter=start_nonboundary*vectorization_length; counter))); out << " Memory dof indices: "; - size_info.print_memory_statistics + size_info.print_memory_statistics (out, MemoryConsumption::memory_consumption (dof_indices)); out << " Memory constraint indicators: "; size_info.print_memory_statistics (out, MemoryConsumption::memory_consumption (constraint_indicator)); out << " Memory plain indices: "; - size_info.print_memory_statistics + size_info.print_memory_statistics (out, MemoryConsumption::memory_consumption (row_starts_plain_indices)+ MemoryConsumption::memory_consumption (plain_dof_indices)); out << " Memory vector partitioner: "; @@ -2169,7 +2100,7 @@ namespace MatrixFreeFunctions } -} // end of namespace MatrixFreeFunctions + } // end of namespace MatrixFreeFunctions } // end of namespace internal DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/include/deal.II/matrix_free/helper_functions.h b/deal.II/include/deal.II/matrix_free/helper_functions.h index 5fccb9d8b6..b9558fa3c0 100644 --- a/deal.II/include/deal.II/matrix_free/helper_functions.h +++ b/deal.II/include/deal.II/matrix_free/helper_functions.h @@ -22,8 +22,6 @@ #include #include -#include - DEAL_II_NAMESPACE_OPEN @@ -100,7 +98,7 @@ namespace MatrixFreeFunctions * to zero. */ void clear(); - + /** * Prints minimum, average, and * maximal memory consumption over the @@ -144,54 +142,38 @@ namespace MatrixFreeFunctions /** * Data type to identify cell type. - */ + */ enum CellType {cartesian=0, affine=1, general=2, undefined=3}; - // ----------------- hash structure -------------------------------- - /** - * A class that is - * used to quickly find out whether two - * vectors of floating point numbers are the - * same without going through all the - * elements: store a hash value for each - * vector. Generate the - * hash value by a sum of all values - * multiplied by random numbers (cast to - * int). Of course, this is not a sure - * criterion and one must manually check for - * equality before this hash is telling - * something useful. However, inequalities are - * easily detected (unless roundoff spoils the - * hash function) - */ - struct HashValue + /** + * A class that is used to compare floating point arrays (e.g. std::vectors, + * Tensor<1,dim>, etc.). The idea of this class is to consider two arrays as + * equal if they are the same within a given tolerance. We use this + * comparator class within an std::map<> of the given arrays. Note that this + * comparison operator does not satisfy all the mathematical properties one + * usually wants to have (consider e.g. the numbers a=0, b=0.1, c=0.2 with + * tolerance 0.15; the operator gives a + struct FPArrayComparator { - // Constructor: sets the size of Number values - // with the typical magnitude that is to be - // expected. - HashValue (const double element_size = 1.); - - // get hash value for a vector of floating - // point numbers (which are assumed to be of - // order of magnitude one). Do this by first - // truncating everything that is smaller than - // the scaling (in order to eliminate noise - // from roundoff errors) and then calling the - // boost hash function - unsigned int operator ()(const std::vector &vec); - - // get hash value for a tensor of rank - // two where the magnitude of the - // entries is given by the parameter - // weight - template - unsigned int operator ()(const Tensor<2,dim,VectorizedArray > - &input, - const bool is_diagonal); - - - const double scaling; + FPArrayComparator (const Number scaling); + + bool operator() (const std::vector &v1, + const std::vector &v2) const; + + template + bool operator ()(const Tensor<1,dim,VectorizedArray > *t1, + const Tensor<1,dim,VectorizedArray > *t2) const; + + template + bool operator ()(const Tensor<2,dim,VectorizedArray > *t1, + const Tensor<2,dim,VectorizedArray > *t2) const; + + Number tolerance; }; // Note: Implementation in matrix_free.templates.h diff --git a/deal.II/include/deal.II/matrix_free/mapping_info.templates.h b/deal.II/include/deal.II/matrix_free/mapping_info.templates.h index 9d1626f45a..f36c058080 100644 --- a/deal.II/include/deal.II/matrix_free/mapping_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/mapping_info.templates.h @@ -202,7 +202,7 @@ namespace MatrixFreeFunctions Utilities::fixed_power(n_q_points_1d[q]); current_data.n_q_points.push_back (n_q_points); - current_data.n_q_points_face.push_back + current_data.n_q_points_face.push_back (Utilities::fixed_power(n_q_points_1d[q])); current_data.quadrature.push_back (Quadrature(quad[my_q][q])); @@ -245,7 +245,7 @@ namespace MatrixFreeFunctions // considered together, this variable holds // the individual info of the last chunk of // cells - CellType cell_t [vectorization_length], + CellType cell_t [vectorization_length], cell_t_prev [vectorization_length]; for (unsigned int j=0; j > hash_collection; - HashValue hash_value (jacobian_size); + FPArrayComparator comparator(jacobian_size); + std::map >*, unsigned int, + FPArrayComparator > cartesians(comparator); + std::map >*, unsigned int, + FPArrayComparator > affines(comparator); // loop over all cells for (unsigned int cell=0; cell (mapping, dummy_fe, current_data.quadrature[fe_index], update_flags_feval)); @@ -341,105 +341,73 @@ namespace MatrixFreeFunctions if (cell_t[j] > most_general_type) most_general_type = cell_t[j]; AssertIndexRange (most_general_type, 3); + unsigned int insert_position = numbers::invalid_unsigned_int; // Cartesian cell with diagonal Jacobian: only // insert the diagonal of the inverse and the - // Jacobian determinant - unsigned int insert_position = numbers::invalid_unsigned_int; - typedef std::vector >::iterator iter; + // Jacobian determinant. We do this by using + // an std::map that collects pointers to all + // Cartesian Jacobians. We need a pointer in + // the std::map because it cannot store data + // based on VectorizedArray (alignment + // issue). We circumvent the problem by + // temporarily filling the next value into the + // cartesian_data field and, in case we did an + // insertion, the data is already in the + // correct place. if (most_general_type == cartesian) { - std::pair >, - VectorizedArray > new_entry; + std::pair >*, + unsigned int> new_entry; + new_entry.second = cartesian_data.size(); + Tensor<1,dim,VectorizedArray > cart; for (unsigned int d=0; d (data.const_jac, true); - std::pair insertion (hash, -insert_position); - iter pos = std::lower_bound (hash_collection.begin(), - hash_collection.end(), - insertion); - - // ok, found a data field with the same - // key. check whether we really hit a - // duplicate, i.e., whether the hash really - // was effective - bool duplicate = true; - if (pos != hash_collection.end() && - pos->first == hash) - { - for (unsigned int d=0; dsecond].first[d][j])> - hash_value.scaling) - duplicate = false; - } - else - duplicate = false; - - // if no duplicate, insert the data - if (duplicate == false) + cart[d] = data.const_jac[d][d]; + cartesian_data.push_back + (std::pair >, + VectorizedArray > + (cart, VectorizedArray())); + new_entry.first = &cartesian_data[new_entry.second].first; + + std::pair >*, + unsigned int, + FPArrayComparator >::iterator, + bool> it = cartesians.insert(new_entry); + if (it.second == false) { - hash_collection.insert (pos, insertion); - cartesian_data.push_back (new_entry); + insert_position = it.first->second; + cartesian_data.resize(new_entry.second); } - // else, remember the position else - insert_position = -pos->second; + insert_position = new_entry.second; } // Constant Jacobian case. same strategy as // before, but with other data fields else if (most_general_type == affine) { - insert_position = affine_data.size(); - - // check whether everything is the same as on - // the previous cell - const unsigned int hash = hash_value.template operator() (data.const_jac, false); - std::pair insertion (hash, -insert_position); - iter pos = std::lower_bound (hash_collection.begin(), - hash_collection.end(), - insertion); - - // ok, found a data field with the same - // key. check whether we really hit a - // duplicate - bool duplicate = true; - if (pos != hash_collection.end() && - pos->first == hash) + std::pair >*, + unsigned int> new_entry; + new_entry.second = affine_data.size(); + affine_data.push_back + (std::pair >, + VectorizedArray > + (data.const_jac, VectorizedArray())); + new_entry.first = &affine_data[new_entry.second].first; + + std::pair >*, + unsigned int, + FPArrayComparator >::iterator, + bool> it = affines.insert(new_entry); + if (it.second == false) { - for (unsigned int d=0; dsecond].first[d] - [e][j])> - hash_value.scaling) - duplicate = false; + insert_position = it.first->second; + affine_data.resize(new_entry.second); } else - duplicate = false; - - if (duplicate == false) - { - hash_collection.insert (pos, insertion); - affine_data.push_back - (std::pair >, - VectorizedArray >(data.const_jac, - make_vectorized_array - (Number(0.)))); - } - else - insert_position = -pos->second; + insert_position = new_entry.second; } // general cell case: first resize the data @@ -598,7 +566,7 @@ namespace MatrixFreeFunctions if (update_flags & update_quadrature_points) { - // eventually we turn to the quadrature points + // eventually we turn to the quadrature points // that we can compress in case we have // Cartesian cells. we also need to reorder // them into arrays of vectorized data types. @@ -636,7 +604,7 @@ namespace MatrixFreeFunctions } } } // end for ( cell < n_macro_cells ) - current_data.rowstart_jacobians.push_back + current_data.rowstart_jacobians.push_back (current_data.jacobians.size()); current_data.rowstart_q_points[n_macro_cells] = current_data.quadrature_points.size(); @@ -711,7 +679,7 @@ namespace MatrixFreeFunctions // and we already have determined that this // cell is either Cartesian or with constant // Jacobian, we have nothing more to do. - if (my_q > 0 && (get_cell_type(cell) == cartesian + if (my_q > 0 && (get_cell_type(cell) == cartesian || get_cell_type(cell) == affine) ) continue; @@ -957,7 +925,7 @@ namespace MatrixFreeFunctions if (general_size_glob > 0) { out << " Memory Jacobian data: "; - size_info.print_memory_statistics + size_info.print_memory_statistics (out, MemoryConsumption::memory_consumption (jacobians) + MemoryConsumption::memory_consumption (JxW_values)); out << " Memory second derivative data: "; @@ -976,7 +944,7 @@ namespace MatrixFreeFunctions if (quad_size_glob > 0) { out << " Memory quadrature points: "; - size_info.print_memory_statistics + size_info.print_memory_statistics (out, MemoryConsumption::memory_consumption (rowstart_q_points) + MemoryConsumption::memory_consumption (quadrature_points)); } @@ -990,10 +958,10 @@ namespace MatrixFreeFunctions const SizeInfo &size_info) const { out << " Cell types: "; - size_info.print_memory_statistics + size_info.print_memory_statistics (out, MemoryConsumption::memory_consumption (cell_type)); out << " Memory transformations compr: "; - size_info.print_memory_statistics + size_info.print_memory_statistics (out, MemoryConsumption::memory_consumption (affine_data) + MemoryConsumption::memory_consumption (cartesian_data)); for (unsigned int j=0; j::get_dof_info (unsigned int dof_index) const +template +inline +unsigned int +MatrixFree::n_constraint_pool_entries() const +{ + return constraint_pool_row_index.size()-1; +} + + + template inline const Number* diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h index c26e375a72..8faeee770d 100644 --- a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h +++ b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h @@ -771,19 +771,38 @@ void MatrixFree::initialize_indices AssertDimension (cell_level_index.size(),size_info.n_macro_cells*vectorization_length); } - // set constraint pool and reorder the indices - constraint_pool_row_index = - constraint_values.constraint_pool_row_index; - constraint_pool_data.resize (constraint_values.constraint_pool_data.size()); - std::copy (constraint_values.constraint_pool_data.begin(), - constraint_values.constraint_pool_data.end(), - constraint_pool_data.begin()); - for (unsigned int no=0; no, unsigned int, + internal::MatrixFreeFunctions::FPArrayComparator >::iterator + it = constraint_values.constraints.begin(), + end = constraint_values.constraints.end(); + std::vector*> + constraints (constraint_values.constraints.size()); + unsigned int length = 0; + for ( ; it != end; ++it) + { + AssertIndexRange(it->second, constraints.size()); + constraints[it->second] = &it->first; + length += it->first.size(); + } + constraint_pool_data.clear(); + constraint_pool_data.reserve (length); + constraint_pool_row_index.reserve(constraint_values.constraints.size()+1); + constraint_pool_row_index.resize(1, 0); + for (unsigned int i=0; ibegin(), + constraints[i]->end()); + constraint_pool_row_index.push_back(constraint_pool_data.size()); } + AssertDimension(constraint_pool_data.size(), length); + for (unsigned int no=0; no::epsilon() * - 1024.) - {} + template + FPArrayComparator::FPArrayComparator (const Number scaling) + : + tolerance (scaling * std::numeric_limits::epsilon() * 1024.) + {} - unsigned int HashValue::operator ()(const std::vector &vec) - { - std::vector mod_vec(vec); - for (unsigned int i=0; i(boost::hash_range (mod_vec.begin(), mod_vec.end())); - } + template + bool + FPArrayComparator::operator() (const std::vector &v1, + const std::vector &v2) const + { + const unsigned int s1 = v1.size(), s2 = v2.size(); + if (s1 < s2) + return true; + else if (s1 > s2) + return false; + else + for (unsigned int i=0; i v2[i] + tolerance) + return false; + return false; + } - template - unsigned int HashValue::operator () - (const Tensor<2,dim,VectorizedArray > &input, - const bool is_diagonal) - { - const unsigned int vectorization_length = - VectorizedArray::n_array_elements; - if (is_diagonal) - { - number mod_tensor [dim][vectorization_length]; - for (unsigned int i=0; i - (boost::hash_range(&mod_tensor[0][0], - &mod_tensor[dim-1][vectorization_length-1]+1)); - } - else - { - number mod_tensor [dim][dim][vectorization_length]; - for (unsigned int i=0; i(boost::hash_range - (&mod_tensor[0][0][0], - &mod_tensor[dim-1][dim-1][vectorization_length-1]+1)); - } - } + template + template + bool + FPArrayComparator:: + operator ()(const Tensor<1,dim,VectorizedArray > *t1, + const Tensor<1,dim,VectorizedArray > *t2) const + { + for (unsigned int d=0; d::n_array_elements;++k) + if ((*t1)[d][k] < (*t2)[d][k] - tolerance) + return true; + else if ((*t1)[d][k] > (*t2)[d][k] + tolerance) + return false; + return false; + } + + + template + template + bool + FPArrayComparator:: + operator ()(const Tensor<2,dim,VectorizedArray > *t1, + const Tensor<2,dim,VectorizedArray > *t2) const + { + for (unsigned int d=0; d::n_array_elements;++k) + if ((*t1)[d][e][k] < (*t2)[d][e][k] - tolerance) + return true; + else if ((*t1)[d][e][k] > (*t2)[d][e][k] + tolerance) + return false; + return false; + } } } -- 2.39.5