From d990f148c0c6ee8fce5edb8cb374aaaca6c77dd7 Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 9 May 2017 21:50:19 -0400 Subject: [PATCH] Get rid of a custom sort implementation. This sorting function is about 10% faster for very large grids (O(10^7) cell error indicators) but otherwise has the same performance as std::sort. This patch replaces this with std::sort and a custom comparator. To see the performance change, try running: #include #include #include #include #include template void qsort_index (const std::vector &a, std::vector &ind, std::size_t l, std::size_t r) { std::size_t i,j; T v; if (r<=l) return; v = a[ind[r]]; i = l-1; j = r; do { do { ++i; } while ((a[ind[i]]>v) && (i0)); if (i 0) { qsort_index(a,ind,l,i-1); } qsort_index(a,ind,i+1,r); } int main() { // switch to length = 10000000 to see timing information constexpr std::size_t length = 20; std::vector indices(length); std::iota(indices.begin(), indices.end(), std::size_t(0)); std::uniform_int_distribution unif(0, length); std::default_random_engine re; std::vector values; for (std::size_t index = 0; index < indices.size(); ++index) { values.push_back(unif(re)); } // for (std::size_t index = 0; index < length; ++index) // { // std::cout << index << ": " // << indices[index] << ": " // << values[indices[index]] << '\n'; // } { auto v = values; auto i = indices; const auto t0 = std::chrono::high_resolution_clock::now(); qsort_index(v, i, 0, length - 1); const auto t1 = std::chrono::high_resolution_clock::now(); std::cout << "milliseconds:" << std::chrono::duration_cast(t1 - t0).count() << std::endl; std::cout << "\nafter qsort_index:\n"; // for (std::size_t index = 0; index < length; ++index) // { // std::cout << index << ": " // << i[index] << ": " // << v[i[index]] << '\n'; // } } { auto v = values; auto i = indices; const auto t0 = std::chrono::high_resolution_clock::now(); std::sort(i.begin(), i.end(), [&](const std::size_t left, const std::size_t right) { return v[left] >= v[right]; }); const auto t1 = std::chrono::high_resolution_clock::now(); std::cout << "milliseconds:" << std::chrono::duration_cast(t1 - t0).count() << std::endl; std::cout << "\nafter std::sort:\n"; // for (std::size_t index = 0; index < length; ++index) // { // std::cout << index << ": " // << i[index] << ": " // << v[i[index]] << '\n'; // } } } --- source/grid/grid_refinement.cc | 74 +++++++--------------------------- 1 file changed, 14 insertions(+), 60 deletions(-) diff --git a/source/grid/grid_refinement.cc b/source/grid/grid_refinement.cc index 870e3de68d..809a6f1d1a 100644 --- a/source/grid/grid_refinement.cc +++ b/source/grid/grid_refinement.cc @@ -132,55 +132,6 @@ namespace } -namespace -{ - /** - * Sorts the vector @p ind as an index vector of @p a in increasing order. - * This implementation of quicksort seems to be faster than the standard - * library version and is needed in @p refine_and_coarsen_optimize. - */ - - template - void qsort_index (const VectorType &a, - std::vector &ind, - int l, - int r) - { - int i,j; - typename VectorType::value_type v; - - if (r<=l) - return; - - v = a(ind[r]); - i = l-1; - j = r; - do - { - do - { - ++i; - } - while ((a(ind[i])>v) && (i0)); - - if (i void GridRefinement::refine (Triangulation &tria, @@ -515,26 +466,29 @@ GridRefinement::refine_and_coarsen_optimize (Triangulation &tria, ExcDimensionMismatch(criteria.size(), tria.n_active_cells())); Assert (criteria.is_non_negative (), ExcNegativeCriteria()); - // get an increasing order on - // the error indicator - std::vector tmp(criteria.size()); - for (unsigned int i=0; i cell_indices(criteria.size()); + std::iota(cell_indices.begin(), cell_indices.end(), 0); - qsort_index (criteria, tmp, 0, criteria.size()-1); + std::sort(cell_indices.begin(), cell_indices.end(), + [&criteria](const unsigned int left, + const unsigned int right) + { + return criteria[left] > criteria[right]; + }); double expected_error_reduction = 0; const double original_error = criteria.l1_norm(); - const unsigned int N = criteria.size(); + const std::size_t N = criteria.size(); // minimize the cost functional discussed in the documentation double min_cost = std::numeric_limits::max(); - unsigned int min_arg = 0; + std::size_t min_arg = 0; - for (unsigned int M = 0; M &tria, } } - refine (tria, criteria, criteria(tmp[min_arg])); + refine (tria, criteria, criteria(cell_indices[min_arg])); } -- 2.39.5