1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
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 3, 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 ;
}
|
