From d990f148c0c6ee8fce5edb8cb374aaaca6c77dd7 Mon Sep 17 00:00:00 2001
From: David Wells <wellsd2@rpi.edu>
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 <algorithm>
#include <chrono>
#include <iostream>
#include <random>
#include <vector>

template <typename T>
void qsort_index (const std::vector<T>     &a,
                  std::vector<std::size_t> &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) && (i<r));
      do
        {
          --j;
        }
      while ((a[ind[j]]<v) && (j>0));

      if (i<j)
        std::swap (ind[i], ind[j]);
      else
        std::swap (ind[i], ind[r]);
    }
  while (i<j);
  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<std::size_t> indices(length);
  std::iota(indices.begin(), indices.end(), std::size_t(0));

  std::uniform_int_distribution<std::size_t> unif(0, length);
  std::default_random_engine re;

  std::vector<std::size_t> 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<std::chrono::milliseconds>(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<std::chrono::milliseconds>(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 <class VectorType>
-  void qsort_index (const VectorType          &a,
-                    std::vector<unsigned int> &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) && (i<r));
-        do
-          {
-            --j;
-          }
-        while ((a(ind[j])<v) && (j>0));
-
-        if (i<j)
-          std::swap (ind[i], ind[j]);
-        else
-          std::swap (ind[i], ind[r]);
-      }
-    while (i<j);
-    qsort_index(a,ind,l,i-1);
-    qsort_index(a,ind,i+1,r);
-  }
-}
-
-
-
 
 template <int dim, class VectorType, int spacedim>
 void GridRefinement::refine (Triangulation<dim,spacedim> &tria,
@@ -515,26 +466,29 @@ GridRefinement::refine_and_coarsen_optimize (Triangulation<dim,spacedim> &tria,
           ExcDimensionMismatch(criteria.size(), tria.n_active_cells()));
   Assert (criteria.is_non_negative (), ExcNegativeCriteria());
 
-  // get an increasing order on
-  // the error indicator
-  std::vector<unsigned int> tmp(criteria.size());
-  for (unsigned int i=0; i<criteria.size(); ++i)
-    tmp[i] = i;
+  // get a decreasing order on the error indicator
+  std::vector<unsigned int> 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<double>::max();
-  unsigned int min_arg = 0;
+  std::size_t min_arg = 0;
 
-  for (unsigned int M = 0; M<criteria.size(); ++M)
+  for (std::size_t M = 0; M<criteria.size(); ++M)
     {
-      expected_error_reduction += (1-std::pow(2.,-1.*order)) * criteria(tmp[M]);
+      expected_error_reduction += (1-std::pow(2.,-1.*order)) * criteria(cell_indices[M]);
 
       const double cost = std::pow(((std::pow(2.,dim)-1)*(1+M)+N),
                                    (double)order/dim) *
@@ -546,7 +500,7 @@ GridRefinement::refine_and_coarsen_optimize (Triangulation<dim,spacedim> &tria,
         }
     }
 
-  refine (tria, criteria, criteria(tmp[min_arg]));
+  refine (tria, criteria, criteria(cell_indices[min_arg]));
 }
 
 
-- 
2.39.5