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


kni
Author(s): Martin Günther
autogenerated on Thu Aug 27 2015 13:40:07