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