new function sorted_indices to replace gsl_sort

lib/vector_ops.cpp

@@ -8,6 +8,8 @@
 #include <cmath>
 #include <string>
 #include <vector>
+#include <algorithm>
+#include <numeric>
 
 using namespace std;
 
@@ -69,3 +71,13 @@ bool triv::is_equidist( double *step, const vector<double>& V )
             return false;
     return true;
 }
+
+//! Compute array of sorted indices.
+
+vector<size_t> triv::sorted_indices(vector<double> const& V)
+{
+    vector<size_t> I(V.size());
+    iota(begin(I), end(I), static_cast<size_t>(0));
+    sort( begin(I), end(I), [&](size_t a, size_t b) { return V[a] < V[b]; } );
+    return I;
+}

lib/vector_ops.hpp

@@ -14,4 +14,5 @@ namespace triv {
                  double tolabs=1e-100, double tolrel=1e-10 );
     bool is_ascending( const vector<double>& V );
     bool is_equidist( double *step, const vector<double>& V );
+    vector<size_t> sorted_indices(vector<double> const& V);
 }