1 files changed, 106 insertions, 0 deletions
diff --git a/gsl-1.9/sort/sortvecind_source.c b/gsl-1.9/sort/sortvecind_source.c
new file mode 100644
index 0000000..04769ee
--- /dev/null
+++ b/gsl-1.9/sort/sortvecind_source.c
@@ -0,0 +1,106 @@
+/*
+ * Implement Heap sort -- direct and indirect sorting
+ * Based on descriptions in Sedgewick "Algorithms in C"
+ *
+ * Copyright (C) 1999  Thomas Walter
+ *
+ * 18 February 2000: Modified for GSL by Brian Gough
+ *
+ * This is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the
+ * Free Software Foundation; either version 2, or (at your option) any
+ * later version.
+ *
+ * This source is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+ * for more details.
+ */
+
+static inline void FUNCTION (index, downheap) (size_t * p, const BASE * data, const size_t stride, const size_t N, size_t k);
+
+static inline void
+FUNCTION (index, downheap) (size_t * p, const BASE * data, const size_t stride, const size_t N, size_t k)
+{
+  const size_t pki = p[k];
+
+  while (k <= N / 2)
+    {
+      size_t j = 2 * k;
+
+      if (j < N && data[p[j] * stride] < data[p[j + 1] * stride])
+        {
+          j++;
+        }
+
+      if (!(data[pki * stride] < data[p[j] * stride])) /* avoid infinite loop if nan */
+        {
+          break;
+        }
+
+      p[k] = p[j];
+
+      k = j;
+    }
+
+  p[k] = pki;
+}
+
+void
+FUNCTION (gsl_sort, index) (size_t * p, const BASE * data, const size_t stride, const size_t n)
+{
+  size_t N;
+  size_t i, k;
+
+  if (n == 0)
+    {
+      return;   /* No data to sort */
+    }
+
+  /* set permutation to identity */
+
+  for (i = 0 ; i < n ; i++)
+    {
+      p[i] = i ;
+    }
+
+  /* We have n_data elements, last element is at 'n_data-1', first at
+     '0' Set N to the last element number. */
+
+  N = n - 1;
+
+  k = N / 2;
+  k++;                          /* Compensate the first use of 'k--' */
+  do
+    {
+      k--;
+      FUNCTION (index, downheap) (p, data, stride, N, k);
+    }
+  while (k > 0);
+
+  while (N > 0)
+    {
+      /* first swap the elements */
+      size_t tmp = p[0];
+      p[0] = p[N];
+      p[N] = tmp;
+
+      /* then process the heap */
+      N--;
+
+      FUNCTION (index, downheap) (p, data, stride, N, 0);
+    }
+}
+
+int
+FUNCTION (gsl_sort_vector, index) (gsl_permutation * permutation, const TYPE (gsl_vector) * v)
+{
+  if (permutation->size != v->size)
+    {
+      GSL_ERROR ("permutation and vector lengths are not equal", GSL_EBADLEN);
+    }
+
+  FUNCTION (gsl_sort, index) (permutation->data, v->data, v->stride, v->size) ;
+  
+  return GSL_SUCCESS ;
+}