515
|
1 // f-sort.cc -*- C++ -*- |
|
2 /* |
|
3 |
1009
|
4 Copyright (C) 1994, 1995 John W. Eaton |
515
|
5 |
|
6 This file is part of Octave. |
|
7 |
|
8 Octave is free software; you can redistribute it and/or modify it |
|
9 under the terms of the GNU General Public License as published by the |
|
10 Free Software Foundation; either version 2, or (at your option) any |
|
11 later version. |
|
12 |
|
13 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
16 for more details. |
|
17 |
|
18 You should have received a copy of the GNU General Public License |
|
19 along with Octave; see the file COPYING. If not, write to the Free |
1315
|
20 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
515
|
21 |
|
22 */ |
|
23 |
|
24 #ifdef HAVE_CONFIG_H |
1192
|
25 #include <config.h> |
515
|
26 #endif |
|
27 |
1352
|
28 #include "defun-dld.h" |
777
|
29 #include "error.h" |
|
30 #include "gripes.h" |
544
|
31 #include "help.h" |
1352
|
32 #include "tree-const.h" |
515
|
33 |
|
34 static void |
|
35 mx_sort (Matrix& m, Matrix& idx, int return_idx) |
|
36 { |
|
37 int nr = m.rows (); |
|
38 int nc = m.columns (); |
|
39 idx.resize (nr, nc); |
|
40 int i, j; |
|
41 |
|
42 if (return_idx) |
|
43 { |
|
44 for (j = 0; j < nc; j++) |
|
45 for (i = 0; i < nr; i++) |
|
46 idx.elem (i, j) = i+1; |
|
47 } |
|
48 |
|
49 for (j = 0; j < nc; j++) |
|
50 { |
|
51 for (int gap = nr/2; gap > 0; gap /= 2) |
|
52 for (i = gap; i < nr; i++) |
|
53 for (int k = i - gap; |
|
54 k >= 0 && m.elem (k, j) > m.elem (k+gap, j); |
|
55 k -= gap) |
|
56 { |
|
57 double tmp = m.elem (k, j); |
|
58 m.elem (k, j) = m.elem (k+gap, j); |
|
59 m.elem (k+gap, j) = tmp; |
|
60 |
|
61 if (return_idx) |
|
62 { |
|
63 double tmp = idx.elem (k, j); |
|
64 idx.elem (k, j) = idx.elem (k+gap, j); |
|
65 idx.elem (k+gap, j) = tmp; |
|
66 } |
|
67 } |
|
68 } |
|
69 } |
|
70 |
|
71 static void |
|
72 mx_sort (RowVector& v, RowVector& idx, int return_idx) |
|
73 { |
|
74 int n = v.capacity (); |
|
75 idx.resize (n); |
|
76 int i; |
|
77 |
|
78 if (return_idx) |
|
79 for (i = 0; i < n; i++) |
|
80 idx.elem (i) = i+1; |
|
81 |
|
82 for (int gap = n/2; gap > 0; gap /= 2) |
|
83 for (i = gap; i < n; i++) |
|
84 for (int k = i - gap; |
|
85 k >= 0 && v.elem (k) > v.elem (k+gap); |
|
86 k -= gap) |
|
87 { |
|
88 double tmp = v.elem (k); |
|
89 v.elem (k) = v.elem (k+gap); |
|
90 v.elem (k+gap) = tmp; |
|
91 |
|
92 if (return_idx) |
|
93 { |
|
94 double tmp = idx.elem (k); |
|
95 idx.elem (k) = idx.elem (k+gap); |
|
96 idx.elem (k+gap) = tmp; |
|
97 } |
|
98 } |
|
99 } |
|
100 |
|
101 static void |
|
102 mx_sort (ComplexMatrix& cm, Matrix& idx, int return_idx) |
|
103 { |
|
104 int nr = cm.rows (); |
|
105 int nc = cm.columns (); |
|
106 idx.resize (nr, nc); |
|
107 int i, j; |
|
108 |
|
109 if (return_idx) |
|
110 { |
|
111 for (j = 0; j < nc; j++) |
|
112 for (i = 0; i < nr; i++) |
|
113 idx.elem (i, j) = i+1; |
|
114 } |
|
115 |
|
116 for (j = 0; j < nc; j++) |
|
117 { |
|
118 for (int gap = nr/2; gap > 0; gap /= 2) |
|
119 for (i = gap; i < nr; i++) |
|
120 for (int k = i - gap; |
|
121 k >= 0 && abs (cm.elem (k, j)) > abs (cm.elem (k+gap, j)); |
|
122 k -= gap) |
|
123 { |
|
124 Complex ctmp = cm.elem (k, j); |
|
125 cm.elem (k, j) = cm.elem (k+gap, j); |
|
126 cm.elem (k+gap, j) = ctmp; |
|
127 |
|
128 if (return_idx) |
|
129 { |
|
130 double tmp = idx.elem (k, j); |
|
131 idx.elem (k, j) = idx.elem (k+gap, j); |
|
132 idx.elem (k+gap, j) = tmp; |
|
133 } |
|
134 } |
|
135 } |
|
136 } |
|
137 |
|
138 static void |
|
139 mx_sort (ComplexRowVector& cv, RowVector& idx, int return_idx) |
|
140 { |
|
141 int n = cv.capacity (); |
|
142 idx.resize (n); |
|
143 int i; |
|
144 |
|
145 if (return_idx) |
|
146 for (i = 0; i < n; i++) |
|
147 idx.elem (i) = i+1; |
|
148 |
|
149 for (int gap = n/2; gap > 0; gap /= 2) |
|
150 for (i = gap; i < n; i++) |
|
151 for (int k = i - gap; |
|
152 k >= 0 && abs (cv.elem (k)) > abs (cv.elem (k+gap)); |
|
153 k -= gap) |
|
154 { |
|
155 Complex tmp = cv.elem (k); |
|
156 cv.elem (k) = cv.elem (k+gap); |
|
157 cv.elem (k+gap) = tmp; |
|
158 |
|
159 if (return_idx) |
|
160 { |
|
161 double tmp = idx.elem (k); |
|
162 idx.elem (k) = idx.elem (k+gap); |
|
163 idx.elem (k+gap) = tmp; |
|
164 } |
|
165 } |
|
166 } |
|
167 |
701
|
168 DEFUN_DLD_BUILTIN ("sort", Fsort, Ssort, 2, 2, |
519
|
169 "[S, I] = sort (X)\n\ |
|
170 \n\ |
|
171 sort the columns of X, optionally return sort index") |
515
|
172 { |
519
|
173 Octave_object retval; |
|
174 |
|
175 int nargin = args.length (); |
515
|
176 |
712
|
177 if (nargin != 1) |
519
|
178 { |
|
179 print_usage ("sort"); |
|
180 return retval; |
|
181 } |
515
|
182 |
|
183 int return_idx = nargout > 1; |
|
184 if (return_idx) |
|
185 retval.resize (2); |
|
186 else |
|
187 retval.resize (1); |
|
188 |
712
|
189 tree_constant arg = args(0); |
620
|
190 |
636
|
191 if (arg.is_real_type ()) |
620
|
192 { |
636
|
193 Matrix m = arg.matrix_value (); |
|
194 |
|
195 if (! error_state) |
620
|
196 { |
636
|
197 if (m.rows () == 1) |
|
198 { |
|
199 int nc = m.columns (); |
|
200 RowVector v (nc); |
|
201 for (int i = 0; i < nc; i++) |
|
202 v.elem (i) = m.elem (0, i); |
|
203 RowVector idx; |
|
204 mx_sort (v, idx, return_idx); |
515
|
205 |
636
|
206 retval(0) = tree_constant (v, 0); |
|
207 if (return_idx) |
|
208 retval(1) = tree_constant (idx, 0); |
|
209 } |
|
210 else |
|
211 { |
1357
|
212 // Sorts m in place, optionally computes index Matrix. |
|
213 |
636
|
214 Matrix idx; |
|
215 mx_sort (m, idx, return_idx); |
515
|
216 |
636
|
217 retval(0) = m; |
|
218 if (return_idx) |
|
219 retval(1) = idx; |
|
220 } |
620
|
221 } |
|
222 } |
636
|
223 else if (arg.is_complex_type ()) |
620
|
224 { |
636
|
225 ComplexMatrix cm = arg.complex_matrix_value (); |
|
226 |
|
227 if (! error_state) |
620
|
228 { |
636
|
229 if (cm.rows () == 1) |
|
230 { |
|
231 int nc = cm.columns (); |
|
232 ComplexRowVector cv (nc); |
|
233 for (int i = 0; i < nc; i++) |
|
234 cv.elem (i) = cm.elem (0, i); |
|
235 RowVector idx; |
|
236 mx_sort (cv, idx, return_idx); |
515
|
237 |
636
|
238 retval(0) = tree_constant (cv, 0); |
|
239 if (return_idx) |
|
240 retval(1) = tree_constant (idx, 0); |
|
241 } |
|
242 else |
|
243 { |
1357
|
244 // Sorts cm in place, optionally computes index Matrix. |
|
245 |
636
|
246 Matrix idx; |
|
247 mx_sort (cm, idx, return_idx); |
515
|
248 |
636
|
249 retval(0) = cm; |
|
250 if (return_idx) |
|
251 retval(1) = idx; |
|
252 } |
620
|
253 } |
|
254 } |
|
255 else |
|
256 { |
636
|
257 gripe_wrong_type_arg ("sort", arg); |
515
|
258 } |
|
259 |
|
260 return retval; |
|
261 } |
|
262 |
|
263 /* |
|
264 ;;; Local Variables: *** |
|
265 ;;; mode: C++ *** |
|
266 ;;; page-delimiter: "^/\\*" *** |
|
267 ;;; End: *** |
|
268 */ |