1 files changed, 146 insertions, 0 deletions
diff --git a/gsl-1.9/sort/subset_source.c b/gsl-1.9/sort/subset_source.c
new file mode 100644
index 0000000..8d02426
--- /dev/null
+++ b/gsl-1.9/sort/subset_source.c
@@ -0,0 +1,146 @@
+/* sort/subset_source.c  
+ * 
+ * Copyright (C) 1999,2000,2001  Thomas Walter, 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.
+ */
+
+/* find the k-th smallest elements of the vector data, in ascending order */
+
+int
+FUNCTION (gsl_sort, smallest) (BASE * dest, const size_t k,
+                               const BASE * src, const size_t stride,
+                               const size_t n)
+{
+  size_t i, j;
+  BASE xbound;
+
+  if (k > n)
+    {
+      GSL_ERROR ("subset length k exceeds vector length n", GSL_EINVAL);
+    }
+
+  if (k == 0 || n == 0)
+    {
+      return GSL_SUCCESS;
+    }
+
+  /* take the first element */
+
+  j = 1;
+  xbound = src[0 * stride];
+  dest[0] = xbound;
+
+  /* examine the remaining elements */
+
+  for (i = 1; i < n; i++)
+    {
+      size_t i1;
+
+      BASE xi = src[i * stride];
+
+      if (j < k)
+        {
+          j++;
+        }
+      else if (xi >= xbound)
+        {
+          continue;
+        }
+
+      for (i1 = j - 1; i1 > 0 ; i1--)
+        {
+          if (xi > dest[i1 - 1])
+            break;
+
+          dest[i1] = dest[i1 - 1];
+        }
+
+      dest[i1] = xi;
+
+      xbound = dest[j-1];
+    }
+
+  return GSL_SUCCESS;
+}
+
+
+int
+FUNCTION (gsl_sort_vector,smallest) (BASE * dest, const size_t k, 
+                                     const TYPE (gsl_vector) * v)
+{
+  return FUNCTION (gsl_sort, smallest) (dest, k, v->data, v->stride, v->size);
+}
+
+int
+FUNCTION (gsl_sort, largest) (BASE * dest, const size_t k,
+                              const BASE * src, const size_t stride,
+                              const size_t n)
+{
+  size_t i, j;
+  BASE xbound;
+
+  if (k > n)
+    {
+      GSL_ERROR ("subset length k exceeds vector length n", GSL_EINVAL);
+    }
+
+  if (k == 0 || n == 0)
+    {
+      return GSL_SUCCESS;
+    }
+
+  /* take the first element */
+
+  j = 1;
+  xbound = src[0 * stride];
+  dest[0] = xbound;
+
+  /* examine the remaining elements */
+
+  for (i = 1; i < n; i++)
+    {
+      size_t i1;
+
+      BASE xi = src[i * stride];
+
+      if (j < k)
+        {
+          j++;
+        }
+      else if (xi <= xbound)
+        {
+          continue;
+        }
+
+      for (i1 = j - 1; i1 > 0 ; i1--)
+        {
+          if (xi < dest[i1 - 1])
+            break;
+
+          dest[i1] = dest[i1 - 1];
+        }
+
+      dest[i1] = xi;
+
+      xbound = dest[j-1];
+    }
+
+  return GSL_SUCCESS;
+}
+
+
+int
+FUNCTION (gsl_sort_vector,largest) (BASE * dest, const size_t k, 
+                                    const TYPE (gsl_vector) * v)
+{
+  return FUNCTION (gsl_sort, largest) (dest, k, v->data, v->stride, v->size);
+}