00001
00002
00003
00004
00005 #ifndef WANT_MATH
00006 #define WANT_MATH
00007 #endif
00008
00009 #include "include.h"
00010
00011 #include "newmatap.h"
00012
00013 #ifdef use_namespace
00014 namespace NEWMAT {
00015 #endif
00016
00017 #ifdef DO_REPORT
00018 #define REPORT { static ExeCounter ExeCount(__LINE__,13); ++ExeCount; }
00019 #else
00020 #define REPORT {}
00021 #endif
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #define DoSimpleSort 17 // when to switch to insert sort
00035 #define MaxDepth 50 // maximum recursion depth
00036
00037 static void MyQuickSortDescending(Real* first, Real* last, int depth);
00038 static void InsertionSortDescending(Real* first, const int length,
00039 int guard);
00040 static Real SortThreeDescending(Real* a, Real* b, Real* c);
00041
00042 static void MyQuickSortAscending(Real* first, Real* last, int depth);
00043 static void InsertionSortAscending(Real* first, const int length,
00044 int guard);
00045
00046
00047 void SortDescending(GeneralMatrix& GM)
00048 {
00049 REPORT
00050 Tracer et("QuickSortDescending");
00051
00052 Real* data = GM.Store(); int max = GM.Storage();
00053
00054 if (max > DoSimpleSort) MyQuickSortDescending(data, data + max - 1, 0);
00055 InsertionSortDescending(data, max, DoSimpleSort);
00056
00057 }
00058
00059 static Real SortThreeDescending(Real* a, Real* b, Real* c)
00060 {
00061
00062 if (*a >= *b)
00063 {
00064 if (*b >= *c) { REPORT return *b; }
00065 else if (*a >= *c) { REPORT Real x = *c; *c = *b; *b = x; return x; }
00066 else { REPORT Real x = *a; *a = *c; *c = *b; *b = x; return x; }
00067 }
00068 else if (*c >= *b) { REPORT Real x = *c; *c = *a; *a = x; return *b; }
00069 else if (*a >= *c) { REPORT Real x = *a; *a = *b; *b = x; return x; }
00070 else { REPORT Real x = *c; *c = *a; *a = *b; *b = x; return x; }
00071 }
00072
00073 static void InsertionSortDescending(Real* first, const int length,
00074 int guard)
00075
00076
00077 {
00078 REPORT
00079 if (length <= 1) return;
00080
00081
00082 Real* f = first; Real v = *f; Real* h = f;
00083 if (guard > length) { REPORT guard = length; }
00084 int i = guard - 1;
00085 while (i--) if (v < *(++f)) { v = *f; h = f; }
00086 *h = *first; *first = v;
00087
00088
00089 i = length - 1; f = first;
00090 while (i--)
00091 {
00092 Real* g = f++; h = f; v = *h;
00093 while (*g < v) *h-- = *g--;
00094 *h = v;
00095 }
00096 }
00097
00098 static void MyQuickSortDescending(Real* first, Real* last, int depth)
00099 {
00100 REPORT
00101 for (;;)
00102 {
00103 const int length = last - first + 1;
00104 if (length < DoSimpleSort) { REPORT return; }
00105 if (depth++ > MaxDepth)
00106 Throw(ConvergenceException("QuickSortDescending fails: "));
00107 Real* centre = first + length/2;
00108 const Real test = SortThreeDescending(first, centre, last);
00109 Real* f = first; Real* l = last;
00110 for (;;)
00111 {
00112 while (*(++f) > test) {}
00113 while (*(--l) < test) {}
00114 if (l <= f) break;
00115 const Real temp = *f; *f = *l; *l = temp;
00116 }
00117 if (f > centre)
00118 { REPORT MyQuickSortDescending(l+1, last, depth); last = f-1; }
00119 else { REPORT MyQuickSortDescending(first, f-1, depth); first = l+1; }
00120 }
00121 }
00122
00123 void SortAscending(GeneralMatrix& GM)
00124 {
00125 REPORT
00126 Tracer et("QuickSortAscending");
00127
00128 Real* data = GM.Store(); int max = GM.Storage();
00129
00130 if (max > DoSimpleSort) MyQuickSortAscending(data, data + max - 1, 0);
00131 InsertionSortAscending(data, max, DoSimpleSort);
00132
00133 }
00134
00135 static void InsertionSortAscending(Real* first, const int length,
00136 int guard)
00137
00138
00139 {
00140 REPORT
00141 if (length <= 1) return;
00142
00143
00144 Real* f = first; Real v = *f; Real* h = f;
00145 if (guard > length) { REPORT guard = length; }
00146 int i = guard - 1;
00147 while (i--) if (v > *(++f)) { v = *f; h = f; }
00148 *h = *first; *first = v;
00149
00150
00151 i = length - 1; f = first;
00152 while (i--)
00153 {
00154 Real* g = f++; h = f; v = *h;
00155 while (*g > v) *h-- = *g--;
00156 *h = v;
00157 }
00158 }
00159 static void MyQuickSortAscending(Real* first, Real* last, int depth)
00160 {
00161 REPORT
00162 for (;;)
00163 {
00164 const int length = last - first + 1;
00165 if (length < DoSimpleSort) { REPORT return; }
00166 if (depth++ > MaxDepth)
00167 Throw(ConvergenceException("QuickSortAscending fails: "));
00168 Real* centre = first + length/2;
00169 const Real test = SortThreeDescending(last, centre, first);
00170 Real* f = first; Real* l = last;
00171 for (;;)
00172 {
00173 while (*(++f) < test) {}
00174 while (*(--l) > test) {}
00175 if (l <= f) break;
00176 const Real temp = *f; *f = *l; *l = temp;
00177 }
00178 if (f > centre)
00179 { REPORT MyQuickSortAscending(l+1, last, depth); last = f-1; }
00180 else { REPORT MyQuickSortAscending(first, f-1, depth); first = l+1; }
00181 }
00182 }
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 void SortSV(DiagonalMatrix& D, Matrix& U, bool ascending)
00193 {
00194 REPORT
00195 Tracer trace("SortSV_DU");
00196 int m = U.Nrows(); int n = U.Ncols();
00197 if (n != D.Nrows()) Throw(IncompatibleDimensionsException(D,U));
00198 Real* u = U.Store();
00199 for (int i=0; i<n; i++)
00200 {
00201 int k = i; Real p = D.element(i);
00202 if (ascending)
00203 {
00204 for (int j=i+1; j<n; j++)
00205 { if (D.element(j) < p) { k = j; p = D.element(j); } }
00206 }
00207 else
00208 {
00209 for (int j=i+1; j<n; j++)
00210 { if (D.element(j) > p) { k = j; p = D.element(j); } }
00211 }
00212 if (k != i)
00213 {
00214 D.element(k) = D.element(i); D.element(i) = p; int j = m;
00215 Real* uji = u + i; Real* ujk = u + k;
00216 if (j) for(;;)
00217 {
00218 p = *uji; *uji = *ujk; *ujk = p;
00219 if (!(--j)) break;
00220 uji += n; ujk += n;
00221 }
00222 }
00223 }
00224 }
00225
00226 void SortSV(DiagonalMatrix& D, Matrix& U, Matrix& V, bool ascending)
00227 {
00228 REPORT
00229 Tracer trace("SortSV_DUV");
00230 int mu = U.Nrows(); int mv = V.Nrows(); int n = D.Nrows();
00231 if (n != U.Ncols()) Throw(IncompatibleDimensionsException(D,U));
00232 if (n != V.Ncols()) Throw(IncompatibleDimensionsException(D,V));
00233 Real* u = U.Store(); Real* v = V.Store();
00234 for (int i=0; i<n; i++)
00235 {
00236 int k = i; Real p = D.element(i);
00237 if (ascending)
00238 {
00239 for (int j=i+1; j<n; j++)
00240 { if (D.element(j) < p) { k = j; p = D.element(j); } }
00241 }
00242 else
00243 {
00244 for (int j=i+1; j<n; j++)
00245 { if (D.element(j) > p) { k = j; p = D.element(j); } }
00246 }
00247 if (k != i)
00248 {
00249 D.element(k) = D.element(i); D.element(i) = p;
00250 Real* uji = u + i; Real* ujk = u + k;
00251 Real* vji = v + i; Real* vjk = v + k;
00252 int j = mu;
00253 if (j) for(;;)
00254 {
00255 p = *uji; *uji = *ujk; *ujk = p; if (!(--j)) break;
00256 uji += n; ujk += n;
00257 }
00258 j = mv;
00259 if (j) for(;;)
00260 {
00261 p = *vji; *vji = *vjk; *vjk = p; if (!(--j)) break;
00262 vji += n; vjk += n;
00263 }
00264 }
00265 }
00266 }
00267
00268
00269
00270
00271 #ifdef use_namespace
00272 }
00273 #endif
00274