1 /* heap.c
2 $Id $
3
4 ----------------------------------------------------------------------
5 This file is part of SPGL1 (Spectral Projected Gradient for L1).
6
7 Copyright (C) 2007 Ewout van den Berg and Michael P. Friedlander,
8 Department of Computer Science, University of British Columbia, Canada.
9 All rights reserved. E-mail: <{ewout78,mpf}@cs.ubc.ca>.
10
11 SPGL1 is free software; you can redistribute it and/or modify it
12 under the terms of the GNU Lesser General Public License as
13 published by the Free Software Foundation; either version 2.1 of the
14 License, or (at your option) any later version.
15
16 SPGL1 is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General
19 Public License for more details.
20
21 You should have received a copy of the GNU Lesser General Public
22 License along with SPGL1; if not, write to the Free Software
23 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
24 USA
25 ----------------------------------------------------------------------
26 */
27
28 #include <assert.h>
29
30 #include "heap.h"
31
32
33 /* ======================================================================= */
34 /* H E L P E R F U N C T I O N S */
35 /* ======================================================================= */
36
37
38 /*
39 \brief Perform the "sift" operation for the heap-sort algorithm.
40
41 A heap is a collection of items arranged in a binary tree. Each
42 child node is smaller than or equal to its parent. If x[k] is the
43 parent, than its children are x[2k+1] and x[2k+2].
44
45 This routine promotes ("sifts up") children that are larger than
46 their parents. Thus, the largest element of the heap is the root node.
47
48 \param[in] root The root index from which to start sifting.
49 \param[in] lastChild The last child (largest node index) in the sift operation.
50 \param[in,out] x The array to be sifted.
51 */
heap_sift(int root,int lastChild,double x[])52 void heap_sift( int root, int lastChild, double x[] )
53 {
54 int child;
55
56 for (; (child = (root * 2) + 1) <= lastChild; root = child) {
57
58 if (child < lastChild)
59 if ( x[child] < x[child+1] )
60 child++;
61
62 if ( x[child] <= x[root] )
63 break;
64
65 swap_double( x[root], x[child] );
66 }
67 }
68
69
70 /*!
71 \brief Perform the "sift" operation for the heap-sort algorithm.
72
73 A heap is a collection of items arranged in a binary tree. Each
74 child node is smaller than or equal to its parent. If x[k] is the
75 parent, than its children are x[2k+1] and x[2k+2].
76
77 This routine promotes ("sifts up") children that are larger than
78 their parents. Thus, the largest element of the heap is the root node.
79
80 Elements in y are associated with those in x and are reordered accordingly.
81
82 \param[in] root The root index from which to start sifting.
83 \param[in] lastChild The last child (largest node index) in the sift operation.
84 \param[in,out] x The array to be sifted.
85 \param[in,out] y The array to be sifted accordingly.
86 */
heap_sift_2(int root,int lastChild,double x[],double y[])87 void heap_sift_2( int root, int lastChild, double x[], double y[] )
88 {
89 int child;
90
91 for (; (child = (root * 2) + 1) <= lastChild; root = child) {
92
93 if (child < lastChild)
94 if ( x[child] < x[child+1] )
95 child++;
96
97 if ( x[child] <= x[root] )
98 break;
99
100 swap_double( x[root], x[child] );
101 swap_double( y[root], y[child] );
102 }
103 }
104
105
106 /*!
107 \brief Discard the largest element and contract the heap.
108
109 On entry, the numElems of the heap are stored in x[0],...,x[numElems-1],
110 and the biggest element is x[0]. The following operations are performed:
111 -# Swap the first and last elements of the heap
112 -# Shorten the length of the heap by one.
113 -# Restore the heap property to the contracted heap.
114 This effectively makes x[0] the next largest element
115 in the list.
116
117 \param[in] numElems The number of elements in the current heap.
118 \param[in,out] x The array to be modified.
119
120 \return The number of elements in the heap after it has been contracted.
121 */
heap_del_max(int numElems,double x[])122 int heap_del_max(int numElems, double x[])
123 {
124 int lastChild = numElems - 1;
125
126 assert(numElems > 0);
127
128 /* Swap the largest element with the lastChild. */
129 swap_double(x[0], x[lastChild]);
130
131 /* Contract the heap size, thereby discarding the largest element. */
132 lastChild--;
133
134 /* Restore the heap property of the contracted heap. */
135 heap_sift(0, lastChild, x);
136
137 return numElems - 1;
138 }
139
140
141 /*!
142 \brief Discard the largest element of x and contract the heaps.
143
144 On entry, the numElems of the heap are stored in x[0],...,x[numElems-1],
145 and the largest element is x[0]. The following operations are performed:
146 -# Swap the first and last elements of both heaps
147 -# Shorten the length of the heaps by one.
148 -# Restore the heap property to the contracted heap x.
149 This effectively makes x[0] the next largest element
150 in the list.
151
152 \param[in] numElems The number of elements in the current heap.
153 \param[in,out] x The array to be modified.
154 \param[in,out] y The array to be modified accordingly
155
156 \return The number of elements in each heap after they have been contracted.
157 */
heap_del_max_2(int numElems,double x[],double y[])158 int heap_del_max_2( int numElems, double x[], double y[] )
159 {
160 int lastChild = numElems - 1;
161
162 assert(numElems > 0);
163
164 /* Swap the largest element with the lastChild. */
165 swap_double( x[0], x[lastChild] );
166 swap_double( y[0], y[lastChild] );
167
168 /* Contract the heap size, thereby discarding the largest element. */
169 lastChild--;
170
171 /* Restore the heap property of the contracted heap. */
172 heap_sift_2( 0, lastChild, x, y );
173
174 return numElems - 1;
175 }
176
177
178 /*!
179
180 \brief Build a heap by adding one element at a time.
181
182 \param[in] n The length of x and ix.
183 \param[in,out] x The array to be heapified.
184
185 */
heap_build(int n,double x[])186 void heap_build( int n, double x[] )
187 {
188 int i;
189
190 for (i = n/2; i >= 0; i--) heap_sift( i, n-1, x );
191 }
192
193
194 /*!
195
196 \brief Build a heap by adding one element at a time.
197
198 \param[in] n The length of x and ix.
199 \param[in,out] x The array to be heapified.
200 \param[in,out] y The array to be reordered in sync. with x.
201
202 */
heap_build_2(int n,double x[],double y[])203 void heap_build_2( int n, double x[], double y[] )
204 {
205 int i;
206
207 for (i = n/2; i >= 0; i--) heap_sift_2( i, n-1, x, y );
208 }
209