[gnuastro-commits] master b85eb68: Sorting of indexs of multiple arrays in multiple threads

[Mohammad Akhlaghi]

branch: master
commit b85eb681ea1752c310b020057ba3cdfcbb69b8b8
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Sorting of indexs of multiple arrays in multiple threads
    
    Until now we only had the `gal_qsort_index_float_increasing' function that
    would use a global variable to keep the array to sort the indexs. This was
    good when were were only dealing with a single array (and possibly many
    arrays of indexs, on multiple threads). However, when each thread needs to
    sort indexs that are defined in a separate array, that function would fail
    to give a reasonable result.
    
    With this commit, we now have the new `gal_qsort_index_multi_d' and
    `gal_qsort_index_multi_i' functions for such situations. For such
    situations, it is necessary to keep a pointer to the values array beside
    each index. A structure (`gal_qsort_index_multi') is used to do this.
    
    Also, to follow a logical naming convension, the old
    `gal_qsort_index_float_increasing' and `gal_qsort_index_float_decreasing'
    functions have been renamed to `gal_qsort_index_single_d' and
    `gal_qsort_index_single_i'.
    
    The whole purpose for this was that I wanted to run `gal_label_oversegment'
    on separate datasets, each on a separate thread. Therefore,
    `gal_label_oversegment' will now look into the flags of the `indexs'
    array. It will only attempt to sort the indexs if they are not sorted. In
    the multi-array and multi-threaded scenario above, indexs can be sorted
    outside `gal_label_oversegment' (with the new `gal_qsort_index_single_d')
    and the `indexs->flag' can be set. Afterwards `gal_label_oversegment' will
    work cleanly on every thread/array separately.
---
 NEWS                 |   7 ++-
 doc/gnuastro.texi    | 135 +++++++++++++++++++++++++++++++++++--------------
 lib/arithmetic.c     |  20 ++++----
 lib/gnuastro/qsort.h |  79 +++++++++++++++++++----------
 lib/label.c          |  17 ++++---
 lib/qsort.c          | 140 ++++++++++++++++++++++++++++++++++-----------------
 lib/statistics.c     |  40 +++++++--------
 7 files changed, 291 insertions(+), 147 deletions(-)

diff --git a/NEWS b/NEWS
index ba8fbfb..560f7af 100644
--- a/NEWS
+++ b/NEWS
@@ -91,7 +91,10 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     gal_pdf_suffix_is_pdf: Returns 1 if given suffix is PDF.
     gal_pdf_write: Writes a dataset into an PDF file.
     gal_pointer_allocate_mmap: Allocate space in a file, not in RAM.
-    gal_qsort_index_float_increasing: Sort indexs in increasing value order.
+    gal_qsort_index_single_d: Sort indexs of single array in decreasing order.
+    gal_qsort_index_single_i: Sort indexs of single array in decreasing order.
+    gal_qsort_index_multi_d: Sort indexs in multiple arrays (on different 
threads).
+    gal_qsort_index_multi_i: Sort indexs in multiple arrays (on different 
threads).
     gal_tiff_name_is_tiff: check if name contains a TIFF suffix.
     gal_tiff_suffix_is_tiff: check if suffix is a TIFF suffix.
     gal_tiff_dir_string_read: convert a string to a TIFF directory number.
@@ -168,6 +171,8 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     gal_pointer_allocate: replaces `gal_data_malloc_array' and
         `gal_data_calloc_array', through an argument you can ask for the
         allocated memory to be cleared or not.
+    gal_qsort_TYPE_i: new name for gal_qsort_TYPE_increasing.
+    gal_qsort_TYPE_d: new name for gal_qsort_TYPE_decreasing.
     gal_statistics_is_sorted: can now also update the bit-flags regarding
         the sorted nature of the input (to optimize future calls to the
         function).
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index a974d3b..e9d236f 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -24353,8 +24353,8 @@ the main parameters of this two-layer tessellation and 
help in benefiting
 from it.
 
 @deftp {Type (C @code{struct})} gal_tile_two_layer_params
-This is the general structure to keep all the necessary parameters for a
-two-layer tessellation.
+The general structure to keep all the necessary parameters for a two-layer
+tessellation.
 
 @example
 struct gal_tile_two_layer_params
@@ -24743,22 +24743,65 @@ both polygons have to be sorted in an anti-clock-wise 
manner.
 @node Qsort functions, Permutations, Polygons, Gnuastro library
 @subsection Qsort functions (@file{qsort.h})
 
-The C programming language comes with the @code{qsort} (Quick sort)
-function. @code{qsort} is a generic function which allows you to sort any
-kind of data structure (not just a single number). To define greater and
-smaller (for sorting), @code{qsort} needs another function, even for simple
-numerical types. To facilitate numerical sorting for Gnuastro's
-programs/libraries, Gnuastro defines a function for each type's increasing
-and decreasing function. You can pass these functions to @code{qsort} when
-your array has the respective type (see @ref{Numeric data types}).
+When sorting a dataset is necessary, the C programming language provides
+the @code{qsort} (Quick sort) function. @code{qsort} is a generic function
+which allows you to sort any kind of data structure (not just a single
+array of numbers). To define ``greater'' and ``smaller'' (for sorting),
address@hidden needs another function, even for simple numerical types. The
+functions introduced in this section are to passed onto @code{qsort}.
+
+The first class of functions below (with @code{TYPE} in their names) can be
+used for sorting a simple numeric array. Just replace @code{TYPE} with the
+dataset's numeric datatype. The second set of functions can be used to sort
+indexs (leave the actual numbers untouched). To use the second set of
+functions, a global variable or structure are also necessary as described
+below.
 
address@hidden {Global variable} {gal_qsort_index_arr}
address@hidden {Global variable} {gal_qsort_index_single}
address@hidden Thread-safety
address@hidden Multi-threaded operation
 Pointer to a floating point array (@code{float *}) to use as a reference in
address@hidden, see the explanation there for
-more.
address@hidden or @code{gal_qsort_index_single_i}, see the
+explanation of these functions for more. Note that if more than one array
+is to be sorted in a multi-threaded operation, the two functions will not
+work as expected. Therefore, if the operation is multi-threaded, but all
+the threads just sort the indexs based on a single array, this global
+variable can safely be used.
 @end deffn
 
address@hidden int gal_qsort_index_float_decreasing (const void @code{*a}, 
const void @code{*b})
address@hidden {Type (C @code{struct})} gal_qsort_index_multi
+Structure to get the sorted indexs of multiple datasets on multiple threads
+with @code{gal_qsort_index_multi_d} or @code{gal_qsort_index_multi_i}. Note
+that the @code{values} array will not be changed by these functions, it is
+only read. Therefore all the @code{values} elements in the (to be sorted)
+array of @code{gal_qsort_index_multi} must point to the same place.
+
address@hidden
+struct gal_qsort_index_multi
address@hidden
+  float *values;         /* Array of values (same in all).      */
+  size_t  index;         /* Index of each element to be sorted. */
address@hidden;
address@hidden example
address@hidden deftp
+
address@hidden int gal_qsort_TYPE_d (const void @code{*a}, const void @code{*b})
+When passed to @code{qsort}, this function will sort a @code{TYPE} array in
+decreasing order (first element will be the largest). Please replace
address@hidden (in the function name) with one of the @ref{Numeric data
+types}, for example @code{gal_qsort_int32_d}, or
address@hidden
address@hidden deftypefun
+
address@hidden int gal_qsort_TYPE_i (const void @code{*a}, const void @code{*b})
+When passed to @code{qsort}, this function will sort a @code{TYPE} array in
+increasing order (first element will be the smallest). Please replace
address@hidden (in the function name) with one of the @ref{Numeric data
+types}, for example @code{gal_qsort_int32_i}, or
address@hidden
address@hidden deftypefun
+
address@hidden int gal_qsort_index_single_d (const void @code{*a}, const void 
@code{*b})
 When passed to @code{qsort}, this function will sort a @code{size_t} array
 based on decreasing values in the @code{gal_qsort_index_arr} single
 precision floating point array. The floating point array will not be
@@ -24776,7 +24819,7 @@ main (void)
   size_t address@hidden, 1, 2, address@hidden;
   float address@hidden,0.2,1.8,address@hidden;
   gal_qsort_index_arr=f;
-  qsort(s, 4, sizeof(size_t), gal_qsort_index_float_decreasing);
+  qsort(s, 4, sizeof(size_t), gal_qsort_index_single_d);
   printf("%zu, %zu, %zu, %zu\n", s[0], s[1], s[2], s[3]);
   return EXIT_SUCCESS;
 @}
@@ -24786,27 +24829,28 @@ main (void)
 The output will be: @code{2, 0, 1, 3}.
 @end deftypefun
 
address@hidden int gal_qsort_index_float_increasing (const void @code{*a}, 
const void @code{*b})
-Similar to @code{gal_qsort_index_float_decreasing}, but will sort the
-indexes such that the values of @code{gal_qsort_index_arr} are in
-increasing order.
address@hidden int gal_qsort_index_single_i (const void @code{*a}, const void 
@code{*b})
+Similar to @code{gal_qsort_index_single_d}, but will sort the indexes such
+that the values of @code{gal_qsort_index_arr} are in increasing order.
 @end deftypefun
 
address@hidden int gal_qsort_index_multi_d (const void @code{*a}, const void 
@code{*b})
+When passed to @code{qsort} with an array of @code{gal_qsort_index_multi},
+this function will sort the array based on the values of the given
+indexs. The sorting will be ordered according to the @code{values} pointer
+of @code{gal_qsort_index_multi}. Note that @code{values} must point to the
+same place in all the structures of the @code{gal_qsort_index_multi} array.
 
address@hidden int gal_qsort_TYPE_increasing (const void @code{*a}, const void 
@code{*b})
-When passed to @code{qsort}, this function will sort an @code{TYPE} array
-in increasing order (first element will be the smallest). Please replace
address@hidden (in the function name) with one of the @ref{Numeric data
-types}, for example @code{gal_qsort_int32_increasing}, or
address@hidden
+This function is only useful when the the indexs of multiple arrays on
+multiple threads are to be sorted. If your program is single threaded, or
+all the indexs belong to a single array (sorting different sub-sets of
+indexs in a single array on multiple threads), it is recommended to use
address@hidden
 @end deftypefun
 
address@hidden int gal_qsort_TYPE_decreasing (const void @code{*a}, const void 
@code{*b})
-When passed to @code{qsort}, this function will sort an @code{TYPE} array
-in decreasing order (first element will be the largest). Please replace
address@hidden (in the function name) with one of the @ref{Numeric data
-types}, for example @code{gal_qsort_int32_decreasing}, or
address@hidden
address@hidden int gal_qsort_index_multi_i (const void @code{*a}, const void 
address@hidden)
+Similar to @code{gal_qsort_index_multi_d}, but the result will be sorted in
+increasing order (first element will have the smallest value).
 @end deftypefun
 
 
@@ -25573,13 +25617,26 @@ The @code{values} dataset must have a 32-bit floating 
point type
 read by this function. @code{indexs} must contain the indexs of the
 elements/pixels that will be over-segmented by this function and have a
 @code{GAL_TYPE_SIZE_T} type, see the description of
address@hidden, above. After this function, @code{indexs} will be
-sorted by the values of the respective pixel in @code{values}. The final
-labels will be written in the respective positions of @code{labels}, which
-must have a @code{GAL_TYPE_INT32} type and be the same size as
address@hidden All local minima (maxima), when @code{min0_max1} is @code{1}
-(@code{0}), are considered rivers (watersheds) and given a label of
address@hidden (see above).
address@hidden, above. The final labels will be written in the
+respective positions of @code{labels}, which must have a
address@hidden type and be the same size as @code{values}.
+
+When @code{indexs} is already sorted, this function will ignore
address@hidden To judge if the dataset is sorted or not (by the values
+the indexs correspond to in @code{values}, not the actual indexs), this
+function will look into the bits of @code{indexs->flag}, for the respective
+bit flags, see @ref{Generic data container}. If @code{indexs} is not
+already sorted, this function will sort it according to the values of the
+respective pixel in @code{values}. The increasing/decreasing order will be
+determined by @code{min0_max1}. Note that if this function is called on
+multiple threads @emph{and} @code{values} points to a different array on
+each thread, this function will not return a reasonable result. In this
+case, please sort @code{indexs} prior to calling this function (see
address@hidden in @ref{Qsort functions}).
+
+When @code{indexs} is decreasing (increasing), or @code{min0_max1} is
address@hidden (@code{0}), local minima (maxima), are considered rivers
+(watersheds) and given a label of @code{GAL_LABEL_RIVER} (see above).
 
 Note that rivers/watersheds will also be formed on the edges of the labeled
 regions or when the labeled pixels touch a blank pixel. Therefore this
@@ -25588,7 +25645,7 @@ efficient, it is thus recommended to use 
@code{gal_blank_present} (with
 @code{updateflag=1}) prior to calling this function (see @ref{Library blank
 values}. Once the flag has been set, no other function (including this one)
 that needs special behavior for blank pixels will have to parse the dataset
-any more.
+to see if it has blank values any more.
 
 If you are sure your dataset doesn't have blank values (by the design of
 your software), to avoid an extra parsing of the dataset and improve
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 2d2a908..6f31f51 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -949,34 +949,34 @@ arithmetic_multioperand(int operator, int flags, 
gal_data_t *list)
   switch(list->type)
     {
     case GAL_TYPE_UINT8:
-      MULTIOPERAND_TYPE_SET(uint8_t,   gal_qsort_uint8_increasing);
+      MULTIOPERAND_TYPE_SET(uint8_t,   gal_qsort_uint8_i);
       break;
     case GAL_TYPE_INT8:
-      MULTIOPERAND_TYPE_SET(int8_t,    gal_qsort_int8_increasing);
+      MULTIOPERAND_TYPE_SET(int8_t,    gal_qsort_int8_i);
       break;
     case GAL_TYPE_UINT16:
-      MULTIOPERAND_TYPE_SET(uint16_t,  gal_qsort_uint16_increasing);
+      MULTIOPERAND_TYPE_SET(uint16_t,  gal_qsort_uint16_i);
       break;
     case GAL_TYPE_INT16:
-      MULTIOPERAND_TYPE_SET(int16_t,   gal_qsort_int16_increasing);
+      MULTIOPERAND_TYPE_SET(int16_t,   gal_qsort_int16_i);
       break;
     case GAL_TYPE_UINT32:
-      MULTIOPERAND_TYPE_SET(uint32_t,  gal_qsort_uint32_increasing);
+      MULTIOPERAND_TYPE_SET(uint32_t,  gal_qsort_uint32_i);
       break;
     case GAL_TYPE_INT32:
-      MULTIOPERAND_TYPE_SET(int32_t,   gal_qsort_int32_increasing);
+      MULTIOPERAND_TYPE_SET(int32_t,   gal_qsort_int32_i);
       break;
     case GAL_TYPE_UINT64:
-      MULTIOPERAND_TYPE_SET(uint64_t,  gal_qsort_uint64_increasing);
+      MULTIOPERAND_TYPE_SET(uint64_t,  gal_qsort_uint64_i);
       break;
     case GAL_TYPE_INT64:
-      MULTIOPERAND_TYPE_SET(int64_t,   gal_qsort_int64_increasing);
+      MULTIOPERAND_TYPE_SET(int64_t,   gal_qsort_int64_i);
       break;
     case GAL_TYPE_FLOAT32:
-      MULTIOPERAND_TYPE_SET(float,     gal_qsort_float32_increasing);
+      MULTIOPERAND_TYPE_SET(float,     gal_qsort_float32_i);
       break;
     case GAL_TYPE_FLOAT64:
-      MULTIOPERAND_TYPE_SET(double,    gal_qsort_float64_increasing);
+      MULTIOPERAND_TYPE_SET(double,    gal_qsort_float64_i);
       break;
     default:
       error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
diff --git a/lib/gnuastro/qsort.h b/lib/gnuastro/qsort.h
index 2b20876..c6e8b76 100644
--- a/lib/gnuastro/qsort.h
+++ b/lib/gnuastro/qsort.h
@@ -1,5 +1,5 @@
 /*********************************************************************
-forqsort -- Functions used by qsort to sort an array.
+qsort -- Functions used by qsort to sort an array.
 This is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -48,77 +48,104 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 
-/* Pointer used to sort the indexs of an array based on their flux
-   (value in this array). */
-extern float *gal_qsort_index_arr;
 
+
+/*****************************************************************/
+/**********        Sorting of actual array        ****************/
+/*****************************************************************/
 int
-gal_qsort_index_float_decreasing(const void * a, const void * b);
+gal_qsort_uint8_d(const void *a, const void *b);
 
 int
-gal_qsort_index_float_increasing(const void * a, const void * b);
+gal_qsort_uint8_i(const void *a, const void *b);
 
+int
+gal_qsort_int8_d(const void *a, const void *b);
 
+int
+gal_qsort_int8_i(const void *a, const void *b);
 
 int
-gal_qsort_uint8_decreasing(const void *a, const void *b);
+gal_qsort_uint16_d(const void *a, const void *b);
 
 int
-gal_qsort_uint8_increasing(const void *a, const void *b);
+gal_qsort_uint16_i(const void *a, const void *b);
 
 int
-gal_qsort_int8_decreasing(const void *a, const void *b);
+gal_qsort_int16_d(const void *a, const void *b);
 
 int
-gal_qsort_int8_increasing(const void *a, const void *b);
+gal_qsort_int16_i(const void *a, const void *b);
 
 int
-gal_qsort_uint16_decreasing(const void *a, const void *b);
+gal_qsort_uint32_d(const void *a, const void *b);
 
 int
-gal_qsort_uint16_increasing(const void *a, const void *b);
+gal_qsort_uint32_i(const void *a, const void *b);
 
 int
-gal_qsort_int16_decreasing(const void *a, const void *b);
+gal_qsort_int32_d(const void *a, const void *b);
 
 int
-gal_qsort_int16_increasing(const void *a, const void *b);
+gal_qsort_int32_i(const void *a, const void *b);
 
 int
-gal_qsort_uint32_decreasing(const void *a, const void *b);
+gal_qsort_uint64_d(const void *a, const void *b);
 
 int
-gal_qsort_uint32_increasing(const void *a, const void *b);
+gal_qsort_uint64_i(const void *a, const void *b);
 
 int
-gal_qsort_int32_decreasing(const void *a, const void *b);
+gal_qsort_int64_d(const void *a, const void *b);
 
 int
-gal_qsort_int32_increasing(const void *a, const void *b);
+gal_qsort_int64_i(const void *a, const void *b);
 
 int
-gal_qsort_uint64_decreasing(const void *a, const void *b);
+gal_qsort_float32_d(const void *a, const void *b);
 
 int
-gal_qsort_uint64_increasing(const void *a, const void *b);
+gal_qsort_float32_i(const void *a, const void *b);
 
 int
-gal_qsort_int64_decreasing(const void *a, const void *b);
+gal_qsort_float64_d(const void *a, const void *b);
 
 int
-gal_qsort_int64_increasing(const void *a, const void *b);
+gal_qsort_float64_i(const void *a, const void *b);
+
+
+
+
+
+/*****************************************************************/
+/***************          Sorting indexs        ******************/
+/*****************************************************************/
+/* Pointer used to sort the indexs of an array based on their flux (value
+   in this array). Note: when EACH THREAD USES A DIFFERENT ARRAY, this is
+   not thread-safe . */
+extern float *gal_qsort_index_single;
+
+
+/* When each thread is working on a different array, we'll need to keep the
+   pointer to the array in question for every index. */
+struct gal_qsort_index_multi
+{
+  float *values; /* Array of values (pointer, so original is not touched). */
+                 /* This should be identical in all elements.              */
+  size_t  index; /* Index of each element to be sorted.                    */
+};
 
 int
-gal_qsort_float32_decreasing(const void *a, const void *b);
+gal_qsort_index_single_d(const void *a, const void *b);
 
 int
-gal_qsort_float32_increasing(const void *a, const void *b);
+gal_qsort_index_single_i(const void *a, const void *b);
 
 int
-gal_qsort_float64_decreasing(const void *a, const void *b);
+gal_qsort_index_multi_d(const void *a, const void *b);
 
 int
-gal_qsort_float64_increasing(const void *a, const void *b);
+gal_qsort_index_multi_i(const void *a, const void *b);
 
 
 
diff --git a/lib/label.c b/lib/label.c
index ca37348..4fd8e8d 100644
--- a/lib/label.c
+++ b/lib/label.c
@@ -236,13 +236,18 @@ gal_label_oversegment(gal_data_t *values, gal_data_t 
*indexs,
   if(indexs->size==0) return 0;
 
 
-  /* Sort the given indexs based on their flux (`gal_qsort_index_arr' is
+  /* If the indexs aren't already sorted (by the value they correspond to),
+     sort them given indexs based on their flux (`gal_qsort_index_arr' is
      defined as static in `gnuastro/qsort.h') */
-  gal_qsort_index_arr=values->array;
-  qsort(indexs->array, indexs->size, sizeof(size_t),
-        min0_max1
-        ? gal_qsort_index_float_decreasing
-        : gal_qsort_index_float_increasing );
+  if( !( (indexs->flag & GAL_DATA_FLAG_SORT_CH)
+        && ( indexs->flag
+             & (GAL_DATA_FLAG_SORTED_I
+                | GAL_DATA_FLAG_SORTED_D) ) ) )
+    {
+      gal_qsort_index_single=values->array;
+      qsort(indexs->array, indexs->size, sizeof(size_t),
+            min0_max1 ? gal_qsort_index_single_d : gal_qsort_index_single_i);
+    }
 
 
   /* Initialize the region we want to over-segment. */
diff --git a/lib/qsort.c b/lib/qsort.c
index 258eddb..59e5a59 100644
--- a/lib/qsort.c
+++ b/lib/qsort.c
@@ -1,5 +1,5 @@
 /*********************************************************************
-forqsort -- Functions used by qsort to sort an array.
+qsort -- Functions used by qsort to sort an array.
 This is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -29,133 +29,112 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/qsort.h>
 
-/* Initialize the array for sorting indexs to NULL. */
-float *gal_qsort_index_arr=NULL;
-
-int
-gal_qsort_index_float_decreasing(const void * a, const void * b)
-{
-  float ta=gal_qsort_index_arr[ *(size_t *)a ];
-  float tb=gal_qsort_index_arr[ *(size_t *)b ];
-  return (tb > ta) - (tb < ta);
-}
-
-int
-gal_qsort_index_float_increasing(const void * a, const void * b)
-{
-  float ta=gal_qsort_index_arr[ *(size_t *)a ];
-  float tb=gal_qsort_index_arr[ *(size_t *)b ];
-  return (ta > tb) - (ta < tb);
-}
-
-
-
-
-
-
 
 
 
 
+/*****************************************************************/
+/**********        Sorting of actual array        ****************/
+/*****************************************************************/
 int
-gal_qsort_uint8_decreasing(const void *a, const void *b)
+gal_qsort_uint8_d(const void *a, const void *b)
 {
   return ( *(uint8_t *)b - *(uint8_t *)a );
 }
 
 int
-gal_qsort_uint8_increasing(const void *a, const void *b)
+gal_qsort_uint8_i(const void *a, const void *b)
 {
   return ( *(uint8_t *)a - *(uint8_t *)b );
 }
 
 int
-gal_qsort_int8_decreasing(const void *a, const void *b)
+gal_qsort_int8_d(const void *a, const void *b)
 {
   return ( *(int8_t *)b - *(int8_t *)a );
 }
 
 int
-gal_qsort_int8_increasing(const void *a, const void *b)
+gal_qsort_int8_i(const void *a, const void *b)
 {
   return ( *(int8_t *)a - *(int8_t *)b );
 }
 
 int
-gal_qsort_uint16_decreasing(const void *a, const void *b)
+gal_qsort_uint16_d(const void *a, const void *b)
 {
   return ( *(uint16_t *)b - *(uint16_t *)a );
 }
 
 int
-gal_qsort_uint16_increasing(const void *a, const void *b)
+gal_qsort_uint16_i(const void *a, const void *b)
 {
   return ( *(uint16_t *)a - *(uint16_t *)b );
 }
 
 int
-gal_qsort_int16_decreasing(const void *a, const void *b)
+gal_qsort_int16_d(const void *a, const void *b)
 {
   return ( *(int16_t *)b - *(int16_t *)a );
 }
 
 int
-gal_qsort_int16_increasing(const void *a, const void *b)
+gal_qsort_int16_i(const void *a, const void *b)
 {
   return ( *(int16_t *)a - *(int16_t *)b );
 }
 
 int
-gal_qsort_uint32_decreasing(const void *a, const void *b)
+gal_qsort_uint32_d(const void *a, const void *b)
 {
   return ( *(uint32_t *)b - *(uint32_t *)a );
 }
 
 int
-gal_qsort_uint32_increasing(const void *a, const void *b)
+gal_qsort_uint32_i(const void *a, const void *b)
 {
   return ( *(uint32_t *)a - *(uint32_t *)b );
 }
 
 int
-gal_qsort_int32_decreasing(const void *a, const void *b)
+gal_qsort_int32_d(const void *a, const void *b)
 {
   return ( *(int32_t *)b - *(int32_t *)a );
 }
 
 int
-gal_qsort_int32_increasing(const void *a, const void *b)
+gal_qsort_int32_i(const void *a, const void *b)
 {
   return ( *(int32_t *)a - *(int32_t *)b );
 }
 
 int
-gal_qsort_uint64_decreasing(const void *a, const void *b)
+gal_qsort_uint64_d(const void *a, const void *b)
 {
   return ( *(uint64_t *)b - *(uint64_t *)a );
 }
 
 int
-gal_qsort_uint64_increasing(const void *a, const void *b)
+gal_qsort_uint64_i(const void *a, const void *b)
 {
   return ( *(uint64_t *)a - *(uint64_t *)b );
 }
 
 
 int
-gal_qsort_int64_decreasing(const void *a, const void *b)
+gal_qsort_int64_d(const void *a, const void *b)
 {
   return ( *(int64_t *)b - *(int64_t *)a );
 }
 
 int
-gal_qsort_int64_increasing(const void *a, const void *b)
+gal_qsort_int64_i(const void *a, const void *b)
 {
   return ( *(int64_t *)a - *(int64_t *)b );
 }
 
 int
-gal_qsort_float32_decreasing(const void *a, const void *b)
+gal_qsort_float32_d(const void *a, const void *b)
 {
   float ta=*(float*)a;
   float tb=*(float*)b;
@@ -163,7 +142,7 @@ gal_qsort_float32_decreasing(const void *a, const void *b)
 }
 
 int
-gal_qsort_float32_increasing(const void *a, const void *b)
+gal_qsort_float32_i(const void *a, const void *b)
 {
   float ta=*(float*)a;
   float tb=*(float*)b;
@@ -171,7 +150,7 @@ gal_qsort_float32_increasing(const void *a, const void *b)
 }
 
 int
-gal_qsort_float64_decreasing(const void *a, const void *b)
+gal_qsort_float64_d(const void *a, const void *b)
 {
   double ta=*(double*)a;
   double tb=*(double*)b;
@@ -179,9 +158,80 @@ gal_qsort_float64_decreasing(const void *a, const void *b)
 }
 
 int
-gal_qsort_float64_increasing(const void *a, const void *b)
+gal_qsort_float64_i(const void *a, const void *b)
 {
   double ta=*(double*)a;
   double tb=*(double*)b;
   return (ta > tb) - (ta < tb);
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*****************************************************************/
+/***************          Sorting indexs        ******************/
+/*****************************************************************/
+/* Initialize the array for sorting indexs to NULL. */
+float *gal_qsort_index_single=NULL;
+
+int
+gal_qsort_index_single_d(const void *a, const void *b)
+{
+  float ta=gal_qsort_index_single[ *(size_t *)a ];
+  float tb=gal_qsort_index_single[ *(size_t *)b ];
+  return (tb > ta) - (tb < ta);
+}
+
+int
+gal_qsort_index_single_i(const void *a, const void *b)
+{
+  float ta=gal_qsort_index_single[ *(size_t *)a ];
+  float tb=gal_qsort_index_single[ *(size_t *)b ];
+  return (ta > tb) - (ta < tb);
+}
+
+int
+gal_qsort_index_multi_d(const void *a, const void *b)
+{
+  /* Define the structures. */
+  struct gal_qsort_index_multi *A = (struct gal_qsort_index_multi *)a;
+  struct gal_qsort_index_multi *B = (struct gal_qsort_index_multi *)b;
+
+  /* For easy reading. */
+  float ta=A->values[ A->index ];
+  float tb=B->values[ B->index ];
+
+  /* Return the result. */
+  return (tb > ta) - (tb < ta);
+}
+
+int
+gal_qsort_index_multi_i(const void *a, const void *b)
+{
+  /* Define the structures. */
+  struct gal_qsort_index_multi *A = (struct gal_qsort_index_multi *)a;
+  struct gal_qsort_index_multi *B = (struct gal_qsort_index_multi *)b;
+
+  /* For easy reading. */
+  float ta=A->values[ A->index ];
+  float tb=B->values[ B->index ];
+
+  /* Return the result. */
+  return (ta > tb) - (ta < tb);
+}
diff --git a/lib/statistics.c b/lib/statistics.c
index ea86156..b476f7e 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -1264,25 +1264,25 @@ gal_statistics_sort_increasing(gal_data_t *input)
     switch(input->type)
       {
       case GAL_TYPE_UINT8:
-        STATISTICS_SORT(gal_qsort_uint8_increasing);    break;
+        STATISTICS_SORT(gal_qsort_uint8_i);    break;
       case GAL_TYPE_INT8:
-        STATISTICS_SORT(gal_qsort_int8_increasing);     break;
+        STATISTICS_SORT(gal_qsort_int8_i);     break;
       case GAL_TYPE_UINT16:
-        STATISTICS_SORT(gal_qsort_uint16_increasing);   break;
+        STATISTICS_SORT(gal_qsort_uint16_i);   break;
       case GAL_TYPE_INT16:
-        STATISTICS_SORT(gal_qsort_int16_increasing);    break;
+        STATISTICS_SORT(gal_qsort_int16_i);    break;
       case GAL_TYPE_UINT32:
-        STATISTICS_SORT(gal_qsort_uint32_increasing);   break;
+        STATISTICS_SORT(gal_qsort_uint32_i);   break;
       case GAL_TYPE_INT32:
-        STATISTICS_SORT(gal_qsort_int32_increasing);    break;
+        STATISTICS_SORT(gal_qsort_int32_i);    break;
       case GAL_TYPE_UINT64:
-        STATISTICS_SORT(gal_qsort_uint64_increasing);   break;
+        STATISTICS_SORT(gal_qsort_uint64_i);   break;
       case GAL_TYPE_INT64:
-        STATISTICS_SORT(gal_qsort_int64_increasing);    break;
+        STATISTICS_SORT(gal_qsort_int64_i);    break;
       case GAL_TYPE_FLOAT32:
-        STATISTICS_SORT(gal_qsort_float32_increasing);  break;
+        STATISTICS_SORT(gal_qsort_float32_i);  break;
       case GAL_TYPE_FLOAT64:
-        STATISTICS_SORT(gal_qsort_float64_increasing);  break;
+        STATISTICS_SORT(gal_qsort_float64_i);  break;
       default:
         error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
               __func__, input->type);
@@ -1307,25 +1307,25 @@ gal_statistics_sort_decreasing(gal_data_t *input)
     switch(input->type)
       {
       case GAL_TYPE_UINT8:
-        STATISTICS_SORT(gal_qsort_uint8_decreasing);    break;
+        STATISTICS_SORT(gal_qsort_uint8_d);    break;
       case GAL_TYPE_INT8:
-        STATISTICS_SORT(gal_qsort_int8_decreasing);     break;
+        STATISTICS_SORT(gal_qsort_int8_d);     break;
       case GAL_TYPE_UINT16:
-        STATISTICS_SORT(gal_qsort_uint16_decreasing);   break;
+        STATISTICS_SORT(gal_qsort_uint16_d);   break;
       case GAL_TYPE_INT16:
-        STATISTICS_SORT(gal_qsort_int16_decreasing);    break;
+        STATISTICS_SORT(gal_qsort_int16_d);    break;
       case GAL_TYPE_UINT32:
-        STATISTICS_SORT(gal_qsort_uint32_decreasing);   break;
+        STATISTICS_SORT(gal_qsort_uint32_d);   break;
       case GAL_TYPE_INT32:
-        STATISTICS_SORT(gal_qsort_int32_decreasing);    break;
+        STATISTICS_SORT(gal_qsort_int32_d);    break;
       case GAL_TYPE_UINT64:
-        STATISTICS_SORT(gal_qsort_uint64_decreasing);   break;
+        STATISTICS_SORT(gal_qsort_uint64_d);   break;
       case GAL_TYPE_INT64:
-        STATISTICS_SORT(gal_qsort_int64_decreasing);    break;
+        STATISTICS_SORT(gal_qsort_int64_d);    break;
       case GAL_TYPE_FLOAT32:
-        STATISTICS_SORT(gal_qsort_float32_decreasing);  break;
+        STATISTICS_SORT(gal_qsort_float32_d);  break;
       case GAL_TYPE_FLOAT64:
-        STATISTICS_SORT(gal_qsort_float64_decreasing);  break;
+        STATISTICS_SORT(gal_qsort_float64_d);  break;
       default:
         error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
               __func__, input->type);