source: ntrip/branches/BNC_2.12/newmat/sort.cpp

Last change on this file was 2013, checked in by mervart, 13 years ago

* empty log message *

File size: 7.4 KB
Line 
1/// \ingroup newmat
2///@{
3
4/// \file sort.cpp
5/// Sorting functions.
6
7// Copyright (C) 1991,2,3,4: R B Davies
8
9#define WANT_MATH
10
11#include "include.h"
12
13#include "newmatap.h"
14
15#ifdef use_namespace
16namespace NEWMAT {
17#endif
18
19#ifdef DO_REPORT
20#define REPORT { static ExeCounter ExeCount(__LINE__,13); ++ExeCount; }
21#else
22#define REPORT {}
23#endif
24
25/******************************** Quick sort ********************************/
26
27// Quicksort.
28// Essentially the method described in Sedgewick s algorithms in C++
29// My version is still partially recursive, unlike Segewick s, but the
30// smallest segment of each split is used in the recursion, so it should
31// not overlead the stack.
32
33// If the process does not seems to be converging an exception is thrown.
34
35
36#define DoSimpleSort 17 // when to switch to insert sort
37#define MaxDepth 50 // maximum recursion depth
38
39static void MyQuickSortDescending(Real* first, Real* last, int depth);
40static void InsertionSortDescending(Real* first, const int length,
41 int guard);
42static Real SortThreeDescending(Real* a, Real* b, Real* c);
43
44static void MyQuickSortAscending(Real* first, Real* last, int depth);
45static void InsertionSortAscending(Real* first, const int length,
46 int guard);
47
48
49void sort_descending(GeneralMatrix& GM)
50{
51 REPORT
52 Tracer et("sort_descending");
53
54 Real* data = GM.Store(); int max = GM.Storage();
55
56 if (max > DoSimpleSort) MyQuickSortDescending(data, data + max - 1, 0);
57 InsertionSortDescending(data, max, DoSimpleSort);
58
59}
60
61static Real SortThreeDescending(Real* a, Real* b, Real* c)
62{
63 // sort *a, *b, *c; return *b; optimise for already sorted
64 if (*a >= *b)
65 {
66 if (*b >= *c) { REPORT return *b; }
67 else if (*a >= *c) { REPORT Real x = *c; *c = *b; *b = x; return x; }
68 else { REPORT Real x = *a; *a = *c; *c = *b; *b = x; return x; }
69 }
70 else if (*c >= *b) { REPORT Real x = *c; *c = *a; *a = x; return *b; }
71 else if (*a >= *c) { REPORT Real x = *a; *a = *b; *b = x; return x; }
72 else { REPORT Real x = *c; *c = *a; *a = *b; *b = x; return x; }
73}
74
75static void InsertionSortDescending(Real* first, const int length,
76 int guard)
77// guard gives the length of the sequence to scan to find first
78// element (eg = length)
79{
80 REPORT
81 if (length <= 1) return;
82
83 // scan for first element
84 Real* f = first; Real v = *f; Real* h = f;
85 if (guard > length) { REPORT guard = length; }
86 int i = guard - 1;
87 while (i--) if (v < *(++f)) { v = *f; h = f; }
88 *h = *first; *first = v;
89
90 // do the sort
91 i = length - 1; f = first;
92 while (i--)
93 {
94 Real* g = f++; h = f; v = *h;
95 while (*g < v) *h-- = *g--;
96 *h = v;
97 }
98}
99
100static void MyQuickSortDescending(Real* first, Real* last, int depth)
101{
102 REPORT
103 for (;;)
104 {
105 const int length = last - first + 1;
106 if (length < DoSimpleSort) { REPORT return; }
107 if (depth++ > MaxDepth)
108 Throw(ConvergenceException("QuickSortDescending fails: "));
109 Real* centre = first + length/2;
110 const Real test = SortThreeDescending(first, centre, last);
111 Real* f = first; Real* l = last;
112 for (;;)
113 {
114 while (*(++f) > test) {}
115 while (*(--l) < test) {}
116 if (l <= f) break;
117 const Real temp = *f; *f = *l; *l = temp;
118 }
119 if (f > centre)
120 { REPORT MyQuickSortDescending(l+1, last, depth); last = f-1; }
121 else { REPORT MyQuickSortDescending(first, f-1, depth); first = l+1; }
122 }
123}
124
125void sort_ascending(GeneralMatrix& GM)
126{
127 REPORT
128 Tracer et("sort_ascending");
129
130 Real* data = GM.Store(); int max = GM.Storage();
131
132 if (max > DoSimpleSort) MyQuickSortAscending(data, data + max - 1, 0);
133 InsertionSortAscending(data, max, DoSimpleSort);
134
135}
136
137static void InsertionSortAscending(Real* first, const int length,
138 int guard)
139// guard gives the length of the sequence to scan to find first
140// element (eg guard = length)
141{
142 REPORT
143 if (length <= 1) return;
144
145 // scan for first element
146 Real* f = first; Real v = *f; Real* h = f;
147 if (guard > length) { REPORT guard = length; }
148 int i = guard - 1;
149 while (i--) if (v > *(++f)) { v = *f; h = f; }
150 *h = *first; *first = v;
151
152 // do the sort
153 i = length - 1; f = first;
154 while (i--)
155 {
156 Real* g = f++; h = f; v = *h;
157 while (*g > v) *h-- = *g--;
158 *h = v;
159 }
160}
161static void MyQuickSortAscending(Real* first, Real* last, int depth)
162{
163 REPORT
164 for (;;)
165 {
166 const int length = last - first + 1;
167 if (length < DoSimpleSort) { REPORT return; }
168 if (depth++ > MaxDepth)
169 Throw(ConvergenceException("QuickSortAscending fails: "));
170 Real* centre = first + length/2;
171 const Real test = SortThreeDescending(last, centre, first);
172 Real* f = first; Real* l = last;
173 for (;;)
174 {
175 while (*(++f) < test) {}
176 while (*(--l) > test) {}
177 if (l <= f) break;
178 const Real temp = *f; *f = *l; *l = temp;
179 }
180 if (f > centre)
181 { REPORT MyQuickSortAscending(l+1, last, depth); last = f-1; }
182 else { REPORT MyQuickSortAscending(first, f-1, depth); first = l+1; }
183 }
184}
185
186//********* sort diagonal matrix & rearrange matrix columns ****************
187
188// used by SVD
189
190// these are for sorting singular values - should be updated with faster
191// sorts that handle exchange of columns better
192// however time is probably not significant compared with SVD time
193
194void SortSV(DiagonalMatrix& D, Matrix& U, bool ascending)
195{
196 REPORT
197 Tracer trace("SortSV_DU");
198 int m = U.Nrows(); int n = U.Ncols();
199 if (n != D.Nrows()) Throw(IncompatibleDimensionsException(D,U));
200 Real* u = U.Store();
201 for (int i=0; i<n; i++)
202 {
203 int k = i; Real p = D.element(i);
204 if (ascending)
205 {
206 for (int j=i+1; j<n; j++)
207 { if (D.element(j) < p) { k = j; p = D.element(j); } }
208 }
209 else
210 {
211 for (int j=i+1; j<n; j++)
212 { if (D.element(j) > p) { k = j; p = D.element(j); } }
213 }
214 if (k != i)
215 {
216 D.element(k) = D.element(i); D.element(i) = p; int j = m;
217 Real* uji = u + i; Real* ujk = u + k;
218 if (j) for(;;)
219 {
220 p = *uji; *uji = *ujk; *ujk = p;
221 if (!(--j)) break;
222 uji += n; ujk += n;
223 }
224 }
225 }
226}
227
228void SortSV(DiagonalMatrix& D, Matrix& U, Matrix& V, bool ascending)
229{
230 REPORT
231 Tracer trace("SortSV_DUV");
232 int mu = U.Nrows(); int mv = V.Nrows(); int n = D.Nrows();
233 if (n != U.Ncols()) Throw(IncompatibleDimensionsException(D,U));
234 if (n != V.Ncols()) Throw(IncompatibleDimensionsException(D,V));
235 Real* u = U.Store(); Real* v = V.Store();
236 for (int i=0; i<n; i++)
237 {
238 int k = i; Real p = D.element(i);
239 if (ascending)
240 {
241 for (int j=i+1; j<n; j++)
242 { if (D.element(j) < p) { k = j; p = D.element(j); } }
243 }
244 else
245 {
246 for (int j=i+1; j<n; j++)
247 { if (D.element(j) > p) { k = j; p = D.element(j); } }
248 }
249 if (k != i)
250 {
251 D.element(k) = D.element(i); D.element(i) = p;
252 Real* uji = u + i; Real* ujk = u + k;
253 Real* vji = v + i; Real* vjk = v + k;
254 int j = mu;
255 if (j) for(;;)
256 {
257 p = *uji; *uji = *ujk; *ujk = p; if (!(--j)) break;
258 uji += n; ujk += n;
259 }
260 j = mv;
261 if (j) for(;;)
262 {
263 p = *vji; *vji = *vjk; *vjk = p; if (!(--j)) break;
264 vji += n; vjk += n;
265 }
266 }
267 }
268}
269
270
271
272
273#ifdef use_namespace
274}
275#endif
276
277///@}
Note: See TracBrowser for help on using the repository browser.