Get rid of a custom sort implementation.

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';
    //   }
  }
}