sort.cpp
Go to the documentation of this file.
1 
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
16 namespace 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 
39 static void MyQuickSortDescending(Real* first, Real* last, int depth);
40 static void InsertionSortDescending(Real* first, const int length,
41  int guard);
42 static Real SortThreeDescending(Real* a, Real* b, Real* c);
43 
44 static void MyQuickSortAscending(Real* first, Real* last, int depth);
45 static void InsertionSortAscending(Real* first, const int length,
46  int guard);
47 
48 
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);
58 
59 }
60 
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 
75 static 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 
100 static 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 
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);
134 
135 }
136 
137 static 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 }
161 static 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 
194 void 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 
228 void 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 
static void test(int n)
Definition: tmtf.cpp:115
Real & element(int, int)
Definition: newmat6.cpp:748
#define MaxDepth
Definition: sort.cpp:37
double Real
Definition: include.h:307
int Nrows() const
Definition: newmat.h:494
#define REPORT
Definition: sort.cpp:22
static void MyQuickSortDescending(Real *first, Real *last, int depth)
Definition: sort.cpp:100
void sort_ascending(GeneralMatrix &GM)
Definition: sort.cpp:125
int Storage() const
Definition: newmat.h:496
static void InsertionSortAscending(Real *first, const int length, int guard)
Definition: sort.cpp:137
Real * Store() const
Definition: newmat.h:497
#define Throw(E)
Definition: myexcept.h:191
The usual rectangular matrix.
Definition: newmat.h:625
Incompatible dimensions exception.
Definition: newmat.h:1998
#define DoSimpleSort
Definition: sort.cpp:36
#define max(a, b)
Definition: kmlExt.cpp:32
FloatVector FloatVector * a
int Ncols() const
Definition: newmat.h:495
Real trace(const BaseMatrix &B)
Definition: newmat.h:2099
Diagonal matrix.
Definition: newmat.h:896
void SortSV(DiagonalMatrix &D, Matrix &U, bool ascending)
Definition: sort.cpp:194
Covergence failure exception.
Definition: newmat.h:1922
static void InsertionSortDescending(Real *first, const int length, int guard)
Definition: sort.cpp:75
void sort_descending(GeneralMatrix &GM)
Definition: sort.cpp:49
The classes for matrices that can contain data are derived from this.
Definition: newmat.h:447
static Real SortThreeDescending(Real *a, Real *b, Real *c)
Definition: sort.cpp:61
static void MyQuickSortAscending(Real *first, Real *last, int depth)
Definition: sort.cpp:161


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:16