Tekkotsu Homepage
Demos
Overview
Downloads
Dev. Resources
Reference
Credits

sort.cpp

Go to the documentation of this file.
00001 //$$ sort.cpp                            Sorting
00002 
00003 // Copyright (C) 1991,2,3,4: R B Davies
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 /******************************** Quick sort ********************************/
00024 
00025 // Quicksort.
00026 // Essentially the method described in Sedgewick s algorithms in C++
00027 // My version is still partially recursive, unlike Segewick s, but the
00028 // smallest segment of each split is used in the recursion, so it should
00029 // not overlead the stack.
00030 
00031 // If the process does not seems to be converging an exception is thrown.
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    // sort *a, *b, *c; return *b; optimise for already sorted
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 // guard gives the length of the sequence to scan to find first
00076 // element (eg = length)
00077 {
00078    REPORT
00079    if (length <= 1) return;
00080 
00081    // scan for first element
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    // do the sort
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 // guard gives the length of the sequence to scan to find first
00138 // element (eg guard = length)
00139 {
00140    REPORT
00141    if (length <= 1) return;
00142 
00143    // scan for first element
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    // do the sort
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 //********* sort diagonal matrix & rearrange matrix columns ****************
00185 
00186 // used by SVD
00187 
00188 // these are for sorting singular values - should be updated with faster
00189 // sorts that handle exchange of columns better
00190 // however time is probably not significant compared with SVD time
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 

newmat11b
Generated Mon May 9 04:54:18 2016 by Doxygen 1.6.3