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

   1 /*
   2  * Implement Heap sort -- direct and indirect sorting
   3  * Based on descriptions in Sedgewick "Algorithms in C"
   4  *
   5  * Copyright (C) 1999  Thomas Walter
   6  *
   7  * 18 February 2000: Modified for GSL by Brian Gough
   8  *
   9  * This is free software; you can redistribute it and/or modify it
  10  * under the terms of the GNU General Public License as published by the
  11  * Free Software Foundation; either version 3, or (at your option) any
  12  * later version.
  13  *
  14  * This source is distributed in the hope that it will be useful, but WITHOUT
  15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  16  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  17  * for more details.
  18  */
  19
  20 static inline void FUNCTION (index, downheap) (size_t * p, const BASE * data, const size_t stride, const size_t N, size_t k);
  21
  22 static inline void
  23 FUNCTION (index, downheap) (size_t * p, const BASE * data, const size_t stride, const size_t N, size_t k)
  24 {
  25   const size_t pki = p[k];
  26
  27   while (k <= N / 2)
  28     {
  29       size_t j = 2 * k;
  30
  31       if (j < N && data[p[j] * stride] < data[p[j + 1] * stride])
  32         {
  33           j++;
  34         }
  35
  36       if (!(data[pki * stride] < data[p[j] * stride])) /* avoid infinite loop if nan */
  37         {
  38           break;
  39         }
  40
  41       p[k] = p[j];
  42
  43       k = j;
  44     }
  45
  46   p[k] = pki;
  47 }
  48
  49 void
  50 FUNCTION (gsl_sort, index) (size_t * p, const BASE * data, const size_t stride, const size_t n)
  51 {
  52   size_t N;
  53   size_t i, k;
  54
  55   if (n == 0)
  56     {
  57       return;   /* No data to sort */
  58     }
  59
  60   /* set permutation to identity */
  61
  62   for (i = 0 ; i < n ; i++)
  63     {
  64       p[i] = i ;
  65     }
  66
  67   /* We have n_data elements, last element is at 'n_data-1', first at
  68      '0' Set N to the last element number. */
  69
  70   N = n - 1;
  71
  72   k = N / 2;
  73   k++;                          /* Compensate the first use of 'k--' */
  74   do
  75     {
  76       k--;
  77       FUNCTION (index, downheap) (p, data, stride, N, k);
  78     }
  79   while (k > 0);
  80
  81   while (N > 0)
  82     {
  83       /* first swap the elements */
  84       size_t tmp = p[0];
  85       p[0] = p[N];
  86       p[N] = tmp;
  87
  88       /* then process the heap */
  89       N--;
  90
  91       FUNCTION (index, downheap) (p, data, stride, N, 0);
  92     }
  93 }
  94
  95 int
  96 FUNCTION (gsl_sort_vector, index) (gsl_permutation * permutation, const TYPE (gsl_vector) * v)
  97 {
  98   if (permutation->size != v->size)
  99     {
 100       GSL_ERROR ("permutation and vector lengths are not equal", GSL_EBADLEN);
 101     }
 102
 103   FUNCTION (gsl_sort, index) (permutation->data, v->data, v->stride, v->size) ;
 104
 105   return GSL_SUCCESS ;
 106 }