source: nanovis/trunk/newmat11/sort.cpp @ 4794

Last change on this file since 4794 was 2096, checked in by ldelgass, 9 years ago

Normalize line endings, set eol-style to native on *.cpp, *.h files

  • Property svn:eol-style set to native
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.