htswanalysis/MACS/lib/gsl/gsl-1.11/sort/sortind.c

   1 /*
   2  * Implement Heap sort -- direct and indirect sorting
   3  * Based on descriptions in Sedgewick "Algorithms in C"
   4  * Copyright (C) 1999  Thomas Walter
   5  *
   6  * 18 February 2000: Modified for GSL by Brian Gough
   7  *
   8  * This is free software; you can redistribute it and/or modify it
   9  * under the terms of the GNU General Public License as published by the
  10  * Free Software Foundation; either version 3, or (at your option) any
  11  * later version.
  12  *
  13  * This source is distributed in the hope that it will be useful, but WITHOUT
  14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16  * for more details.
  17  */
  18
  19 #include <config.h>
  20 #include <stdlib.h>
  21 #include <gsl/gsl_heapsort.h>
  22
  23 static inline void downheap (size_t * p, const void *data, const size_t size, const size_t N, size_t k, gsl_comparison_fn_t compare);
  24
  25 #define CMP(data,size,j,k) (compare((const char *)(data) + (size) * (j), (const char *)(data) + (size) * (k)))
  26
  27 static inline void
  28 downheap (size_t * p, const void *data, const size_t size, const size_t N, size_t k, gsl_comparison_fn_t compare)
  29 {
  30   const size_t pki = p[k];
  31
  32   while (k <= N / 2)
  33     {
  34       size_t j = 2 * k;
  35
  36       if (j < N && CMP (data, size, p[j], p[j + 1]) < 0)
  37         {
  38           j++;
  39         }
  40
  41       if (CMP (data, size, pki, p[j]) >= 0)
  42         {
  43           break;
  44         }
  45
  46       p[k] = p[j];
  47
  48       k = j;
  49     }
  50
  51   p[k] = pki;
  52 }
  53
  54 int
  55 gsl_heapsort_index (size_t * p, const void *data, size_t count, size_t size, gsl_comparison_fn_t compare)
  56 {
  57   /* Sort the array in ascending order. This is a true inplace
  58      algorithm with N log N operations. Worst case (an already sorted
  59      array) is something like 20% slower */
  60
  61   size_t i, k, N;
  62
  63   if (count == 0)
  64     {
  65       return GSL_SUCCESS;       /* No data to sort */
  66     }
  67
  68   for (i = 0; i < count; i++)
  69     {
  70       p[i] = i ;                /* set permutation to identity */
  71     }
  72
  73   /* We have n_data elements, last element is at 'n_data-1', first at
  74      '0' Set N to the last element number. */
  75
  76   N = count - 1;
  77
  78   k = N / 2;
  79   k++;                          /* Compensate the first use of 'k--' */
  80   do
  81     {
  82       k--;
  83       downheap (p, data, size, N, k, compare);
  84     }
  85   while (k > 0);
  86
  87   while (N > 0)
  88     {
  89       /* first swap the elements */
  90       size_t tmp = p[0];
  91       p[0] = p[N];
  92       p[N] = tmp;
  93
  94       /* then process the heap */
  95       N--;
  96
  97       downheap (p, data, size, N, 0, compare);
  98     }
  99
 100   return GSL_SUCCESS;
 101 }