diff options
Diffstat (limited to 'gsl-1.9/sort/subsetind_source.c')
-rw-r--r-- | gsl-1.9/sort/subsetind_source.c | 146 |
1 files changed, 146 insertions, 0 deletions
diff --git a/gsl-1.9/sort/subsetind_source.c b/gsl-1.9/sort/subsetind_source.c new file mode 100644 index 0000000..5505fdf --- /dev/null +++ b/gsl-1.9/sort/subsetind_source.c @@ -0,0 +1,146 @@ +/* sort/subsetind_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_index) (size_t * p, 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]; + p[0] = 0; + + /* 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 > src[p[i1 - 1] * stride]) + break; + + p[i1] = p[i1 - 1]; + } + + p[i1] = i; + + xbound = src[p[j-1] * stride]; + } + + return GSL_SUCCESS; +} + + +int +FUNCTION (gsl_sort_vector,smallest_index) (size_t * p, const size_t k, + const TYPE (gsl_vector) * v) +{ + return FUNCTION (gsl_sort, smallest_index) (p, k, v->data, v->stride, v->size); +} + +int +FUNCTION (gsl_sort, largest_index) (size_t * p, 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]; + p[0] = 0; + + /* 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 < src[stride * p[i1 - 1]]) + break; + + p[i1] = p[i1 - 1]; + } + + p[i1] = i; + + xbound = src[stride * p[j-1]]; + } + + return GSL_SUCCESS; +} + + +int +FUNCTION (gsl_sort_vector,largest_index) (size_t * p, const size_t k, + const TYPE (gsl_vector) * v) +{ + return FUNCTION (gsl_sort, largest_index) (p, k, v->data, v->stride, v->size); +} |