moab
TupleList.cpp
Go to the documentation of this file.
00001 #include <string.h>
00002 #include <limits.h>
00003 #include <stdio.h>
00004 #include <stdlib.h>
00005 #include <stdarg.h>
00006 #include <iostream>
00007 
00008 #include "moab/TupleList.hpp"
00009 
00010 namespace moab {
00011 
00012 void fail(const char *fmt, ...)
00013 {
00014   va_list ap;
00015   va_start(ap, fmt);
00016   vfprintf(stderr, fmt, ap);
00017   va_end(ap);
00018   exit(1);
00019 }
00020 
00021 TupleList::buffer::buffer(size_t sz)
00022 {
00023   ptr = NULL;
00024   buffSize = 0;
00025   this->buffer_init_(sz, __FILE__);
00026 }
00027 
00028 TupleList::buffer::buffer()
00029 {
00030   buffSize = 0;
00031   ptr = NULL;
00032 }
00033 
00034 void TupleList::buffer::buffer_init_(size_t sizeIn, const char *file)
00035 {
00036   this->buffSize = sizeIn;
00037   void *res = malloc(this->buffSize);
00038   if (!res && buffSize > 0)
00039     fail("%s: allocation of %d bytes failed\n", file, (int) buffSize);
00040   ptr = (char*) res;
00041 }
00042 
00043 void TupleList::buffer::buffer_reserve_(size_t min, const char *file)
00044 {
00045   if (this->buffSize < min)
00046   {
00047     size_t newSize = this->buffSize;
00048     newSize += newSize / 2 + 1;
00049     if (newSize < min)
00050       newSize = min;
00051     void *res = realloc(ptr, newSize);
00052     if (!res && newSize > 0)
00053       fail("%s: reallocation of %d bytes failed\n", file, newSize);
00054     ptr = (char*) res;
00055     this->buffSize = newSize;
00056   }
00057 }
00058 
00059 void TupleList::buffer::reset()
00060 {
00061   free(ptr);
00062   ptr = NULL;
00063   buffSize = 0;
00064 }
00065 
00066 TupleList::TupleList(uint p_mi, uint p_ml, uint p_mul, uint p_mr, uint p_max)
00067 {
00068   vi = NULL;
00069   vl = NULL;
00070   vul = NULL;
00071   vr = NULL;
00072   initialize(p_mi, p_ml, p_mul, p_mr, p_max);
00073 }
00074 
00075 TupleList::TupleList()
00076 {
00077   vi = NULL;
00078   vl = NULL;
00079   vul = NULL;
00080   vr = NULL;
00081   disableWriteAccess();
00082 }
00083 
00084 // Allocates space for the tuple list in memory according to parameters
00085 void TupleList::initialize(uint p_mi, uint p_ml, uint p_mul, uint p_mr,
00086     uint p_max)
00087 {
00088   this->n = 0;
00089   this->max = p_max;
00090   this->mi = p_mi;
00091   this->ml = p_ml;
00092   this->mul = p_mul;
00093   this->mr = p_mr;
00094   size_t sz;
00095 
00096   if (max * mi > 0)
00097   {
00098     sz = max * mi * sizeof(sint);
00099     void *resi = malloc(sz);
00100     if (!resi && max * mi > 0)
00101       fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
00102     vi = (sint*) resi;
00103   }
00104   else
00105     vi = NULL;
00106   if (max * ml > 0)
00107   {
00108     sz = max * ml * sizeof(slong);
00109     void *resl = malloc(sz);
00110     if (!resl && max * ml > 0)
00111       fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
00112     vl = (slong*) resl;
00113   }
00114   else
00115     vl = NULL;
00116   if (max * mul > 0)
00117   {
00118     sz = max * mul * sizeof(ulong);
00119     void *resu = malloc(sz);
00120     if (!resu && max * mul > 0)
00121       fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
00122     vul = (ulong*) resu;
00123   }
00124   else
00125     vul = NULL;
00126   if (max * mr > 0)
00127   {
00128     sz = max * mr * sizeof(realType);
00129     void *resr = malloc(sz);
00130     if (!resr && max * ml > 0)
00131       fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
00132     vr = (realType*) resr;
00133   }
00134   else
00135     vr = NULL;
00136 
00137   //Begin with write access disabled
00138   this->disableWriteAccess();
00139 
00140   //Set read variables
00141   vi_rd = vi;
00142   vl_rd = vl;
00143   vul_rd = vul;
00144   vr_rd = vr;
00145 }
00146 
00147 // Resizes a tuplelist to the given uint max
00148 ErrorCode TupleList::resize(uint maxIn)
00149 {
00150   this->max = maxIn;
00151   size_t sz;
00152 
00153   if (vi || (max * mi > 0))
00154   {
00155     sz = max * mi * sizeof(sint);
00156     void *resi = realloc(vi, sz);
00157     if (!resi && max * mi > 0)
00158     {
00159       fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
00160       return moab::MB_MEMORY_ALLOCATION_FAILED;
00161     }
00162     vi = (sint*) resi;
00163   }
00164   if (vl || (max * ml > 0))
00165   {
00166     sz = max * ml * sizeof(slong);
00167     void *resl = realloc(vl, sz);
00168     if (!resl && max * ml > 0)
00169     {
00170       fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
00171       return moab::MB_MEMORY_ALLOCATION_FAILED;
00172     }
00173     vl = (slong*) resl;
00174   }
00175   if (vul || (max * mul > 0))
00176   {
00177     sz = max * mul * sizeof(ulong);
00178     void *resu = realloc(vul, sz);
00179     if (!resu && max * mul > 0)
00180     {
00181       fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
00182       return moab::MB_MEMORY_ALLOCATION_FAILED;
00183     }
00184     vul = (ulong*) resu;
00185   }
00186   if (vr || (max * mr > 0))
00187   {
00188     sz = max * mr * sizeof(realType);
00189     void *resr = realloc(vr, sz);
00190     if (!resr && max * mr > 0)
00191     {
00192       fail("%s: allocation of %d bytes failed\n", __FILE__, (int) sz);
00193       return moab::MB_MEMORY_ALLOCATION_FAILED;
00194     }
00195     vr = (realType*) resr;
00196   }
00197 
00198   //Set read variables
00199   vi_rd = vi;
00200   vl_rd = vl;
00201   vul_rd = vul;
00202   vr_rd = vr;
00203 
00204   //Set the write variables if necessary
00205   if (writeEnabled)
00206   {
00207     vi_wr = vi;
00208     vl_wr = vl;
00209     vul_wr = vul;
00210     vr_wr = vr;
00211   }
00212   return moab::MB_SUCCESS;
00213 }
00214 
00215 // Frees the memory used by the tuplelist
00216 void TupleList::reset()
00217 {
00218   //free up the pointers
00219   free(vi);
00220   free(vl);
00221   free(vul);
00222   free(vr);
00223   //Set them all to null
00224   vr = NULL;
00225   vi = NULL;
00226   vul = NULL;
00227   vl = NULL;
00228   //Set the read and write pointers to null
00229   disableWriteAccess();
00230   vi_rd = NULL;
00231   vl_rd = NULL;
00232   vul_rd = NULL;
00233   vr_rd = NULL;
00234 }
00235 
00236 // Increments n; if n>max, increase the size of the tuplelist
00237 void TupleList::reserve()
00238 {
00239   n++;
00240   while (n > max)
00241     resize((max ? max + max / 2 + 1 : 2));
00242   last_sorted = -1;
00243 }
00244 
00245 // Given the value and the position in the field, finds the index of the tuple
00246 // to which the value belongs
00247 int TupleList::find(unsigned int key_num, sint value)
00248 {
00249   ulong uvalue = (ulong) value;
00250   if (!(key_num > mi))
00251   {
00252     // Binary search: only if the tuple_list is sorted
00253     if (last_sorted == (int) key_num)
00254     {
00255       int lb = 0, ub = n, index; // lb=lower bound, ub=upper bound, index=mid
00256       for (; lb <= ub;)
00257       {
00258         index = (lb + ub) / 2;
00259         if (vi[index * mi + key_num] == (long) uvalue)
00260           return index;
00261         else if (vi[index * mi + key_num] > (long) uvalue)
00262           ub = index - 1;
00263         else if (vi[index * mi + key_num] < (long) uvalue)
00264           lb = index + 1;
00265       }
00266     }
00267     else
00268     {
00269       // Sequential search: if tuple_list is not sorted
00270       for (long index = 0; index < n; index++)
00271       {
00272         if (vi[index * mi + key_num] == (long) uvalue)
00273           return index;
00274       }
00275     }
00276   }
00277   return -1; // If the value wasn't present or an invalid key was given
00278 }
00279 
00280 int TupleList::find(unsigned int key_num, slong value)
00281 {
00282   ulong uvalue = (ulong) value;
00283   if (!(key_num > ml))
00284   {
00285     if (last_sorted - mi == key_num)
00286     {
00287       int lb = 0, ub = n, index; // lb=lower bound, ub=upper bound, index=mid
00288       for (; lb <= ub;)
00289       {
00290         index = (lb + ub) / 2;
00291         if (vl[index * ml + key_num] == (long) uvalue)
00292           return index;
00293         else if (vl[index * ml + key_num] > (long) uvalue)
00294           ub = index - 1;
00295         else if (vl[index * ml + key_num] < (long) uvalue)
00296           lb = index + 1;
00297       }
00298     }
00299     else
00300     {
00301       // Sequential search: if tuple_list is not sorted
00302       for (uint index = 0; index < n; index++)
00303       {
00304         if (vl[index * ml + key_num] == (long) uvalue)
00305           return index;
00306       }
00307     }
00308   }
00309   return -1; // If the value wasn't present or an invalid key was given
00310 }
00311 
00312 int TupleList::find(unsigned int key_num, ulong value)
00313 {
00314   if (!(key_num > mul))
00315   {
00316     if (last_sorted - mi - ml == key_num)
00317     {
00318       int lb = 0, ub = n - 1, index; // lb=lower bound, ub=upper bound, index=mid
00319       for (; lb <= ub;)
00320       {
00321         index = (lb + ub) / 2;
00322         if (vul[index * mul + key_num] == value)
00323           return index;
00324         else if (vul[index * mul + key_num] > value)
00325           ub = index - 1;
00326         else if (vul[index * mul + key_num] < value)
00327           lb = index + 1;
00328       }
00329     }
00330     else
00331     {
00332       // Sequential search: if tuple_list is not sorted
00333       for (uint index = 0; index < n; index++)
00334       {
00335         if (vul[index * mul + key_num] == value)
00336           return index;
00337       }
00338     }
00339   }
00340   return -1; // If the value wasn't present or an invalid key was given
00341 }
00342 
00343 int TupleList::find(unsigned int key_num, realType value)
00344 {
00345   if (!key_num > mr)
00346   {
00347     // Sequential search: TupleList cannot be sorted by reals
00348     for (uint index = 0; index < n; index++)
00349     {
00350       if (vr[index * mr + key_num] == value)
00351         return index;
00352     }
00353   }
00354   return -1; // If the value wasn't present or an invalid key was given
00355 }
00356 
00357 sint TupleList::get_sint(unsigned int index, unsigned int m)
00358 {
00359   if (mi > m && n > index)
00360     return vi[index * mi + m];
00361   return 0;
00362 }
00363 
00364 slong TupleList::get_int(unsigned int index, unsigned int m)
00365 {
00366   if (ml > m && n > index)
00367     return vl[index * ml + m];
00368   return 0;
00369 }
00370 
00371 ulong TupleList::get_ulong(unsigned int index, unsigned int m)
00372 {
00373   if (mul > m && n > index)
00374     return vul[index * mul + m];
00375   return 0;
00376 }
00377 
00378 realType TupleList::get_double(unsigned int index, unsigned int m)
00379 {
00380   if (mr > m && n > index)
00381     return vr[index * mr + m];
00382   return 0;
00383 }
00384 
00385 ErrorCode TupleList::get(unsigned int index, const sint *&sp, const slong *&ip,
00386     const ulong *&lp, const realType *&dp)
00387 {
00388   if (index <= n)
00389   {
00390     if (mi)
00391       *&sp = &vi[index * mi];
00392     else
00393       *&sp = NULL;
00394     if (ml)
00395       *&ip = &vl[index * ml];
00396     else
00397       *&ip = NULL;
00398     if (mul)
00399       *&lp = &vul[index * mul];
00400     else
00401       *&lp = NULL;
00402     if (mr)
00403       *&dp = &vr[index * mr];
00404     else
00405       *&dp = NULL;
00406 
00407     return MB_SUCCESS;
00408   }
00409   return MB_FAILURE;
00410 }
00411 
00412 unsigned int TupleList::push_back(sint *sp, slong *ip, ulong *lp, realType *dp)
00413 {
00414   reserve();
00415   if (mi)
00416     memcpy(&vi[mi * (n - 1)], sp, mi * sizeof(sint));
00417   if (ml)
00418     memcpy(&vl[ml * (n - 1)], ip, ml * sizeof(long));
00419   if (mul)
00420     memcpy(&vul[mul * (n - 1)], lp, mul * sizeof(ulong));
00421   if (mr)
00422     memcpy(&vr[mr * (n - 1)], dp, mr * sizeof(realType));
00423 
00424   last_sorted = -1;
00425   return n - 1;
00426 }
00427 
00428 void TupleList::enableWriteAccess()
00429 {
00430   writeEnabled = true;
00431   last_sorted = -1;
00432   vi_wr = vi;
00433   vl_wr = vl;
00434   vul_wr = vul;
00435   vr_wr = vr;
00436 }
00437 
00438 void TupleList::disableWriteAccess()
00439 {
00440   writeEnabled = false;
00441   vi_wr = NULL;
00442   vl_wr = NULL;
00443   vul_wr = NULL;
00444   vr_wr = NULL;
00445 }
00446 
00447 void TupleList::getTupleSize(uint &mi_out, uint &ml_out, uint &mul_out,
00448     uint &mr_out) const
00449 {
00450   mi_out = mi;
00451   ml_out = ml;
00452   mul_out = mul;
00453   mr_out = mr;
00454 }
00455 
00456 uint TupleList::inc_n()
00457 {
00458   //Check for direct write access
00459   if (!writeEnabled)
00460   {
00461     enableWriteAccess();
00462   }
00463   n++;
00464   return n;
00465 }
00466 
00467 void TupleList::set_n(uint n_in)
00468 {
00469   //Check for direct write access;
00470   if (!writeEnabled)
00471   {
00472     enableWriteAccess();
00473   }
00474   n = n_in;
00475 }
00476 
00477 void TupleList::print(const char *name) const
00478 {
00479   std::cout << "Printing Tuple " << name << "===================" << std::endl;
00480   unsigned long i = 0, l = 0, ul = 0, r = 0;
00481   for (uint k = 0; k < n; k++)
00482   {
00483     for (uint j = 0; j < mi; j++)
00484     {
00485       std::cout << vi[i++] << " | ";
00486     }
00487     for (uint j = 0; j < ml; j++)
00488     {
00489       std::cout << vl[l++] << " | ";
00490     }
00491     for (uint j = 0; j < mul; j++)
00492     {
00493       std::cout << vul[ul++] << " | ";
00494     }
00495     for (uint j = 0; j < mr; j++)
00496     {
00497       std::cout << vr[r++] << " | ";
00498     }
00499     std::cout << std::endl;
00500   }
00501   std::cout << "=======================================" << std::endl
00502       << std::endl;
00503 }
00504 
00505 void TupleList::permute(uint *perm, void *work)
00506 {
00507   const unsigned int_size = mi * sizeof(sint), long_size = ml * sizeof(slong),
00508       ulong_size = mul * sizeof(ulong), real_size = mr * sizeof(realType);
00509   if (mi)
00510   {
00511     uint *p = perm, *pe = p + n;
00512     char *sorted = (char *) work;
00513     while (p != pe)
00514       memcpy((void *) sorted, &vi[mi * (*p++)], int_size), sorted += int_size;
00515     memcpy(vi, work, int_size * n);
00516   }
00517   if (ml)
00518   {
00519     uint *p = perm, *pe = p + n;
00520     char *sorted = (char *) work;
00521     while (p != pe)
00522       memcpy((void *) sorted, &vl[ml * (*p++)], long_size), sorted += long_size;
00523     memcpy(vl, work, long_size * n);
00524   }
00525   if (mul)
00526   {
00527     uint *p = perm, *pe = p + n;
00528     char *sorted = (char *) work;
00529     while (p != pe)
00530       memcpy((void *) sorted, &vul[mul * (*p++)], ulong_size), sorted +=
00531           ulong_size;
00532     memcpy(vul, work, ulong_size * n);
00533   }
00534   if (mr)
00535   {
00536     uint *p = perm, *pe = p + n;
00537     char *sorted = (char *) work;
00538     while (p != pe)
00539       memcpy((void *) sorted, &vr[mr * (*p++)], real_size), sorted += real_size;
00540     memcpy(vr, work, real_size * n);
00541   }
00542 }
00543 
00544 # define umax_2(a, b) (((a)>(b)) ? (a):(b))
00545 
00546 ErrorCode TupleList::sort(uint key, TupleList::buffer *buf)
00547 {
00548   const unsigned int_size = mi * sizeof(sint);
00549   const unsigned long_size = ml * sizeof(slong);
00550   const unsigned ulong_size = mul * sizeof(ulong);
00551   const unsigned real_size = mr * sizeof(realType);
00552   const unsigned width = umax_2(umax_2(int_size,long_size),
00553       umax_2(ulong_size,real_size));
00554   const unsigned data_size =
00555       key >= mi ? sizeof(SortData<long> ) : sizeof(SortData<uint> );
00556 
00557   uint work_min = n * umax_2(2*data_size,sizeof(sint)+width);
00558   uint *work;
00559   buf->buffer_reserve(work_min);
00560   work = (uint *) buf->ptr;
00561   if (key < mi)
00562     index_sort((uint *) &vi[key], n, mi, work, (SortData<uint>*) work);
00563   else if (key < mi + ml)
00564     index_sort((long*) &vl[key - mi], n, ml, work, (SortData<long>*) work);
00565   else if (key < mi + ml + mul)
00566     index_sort((ulong*) &vul[key - mi - ml], n, mul, work,
00567         (SortData<ulong>*) work);
00568   else
00569     return MB_NOT_IMPLEMENTED;
00570 
00571   permute(work, work + n);
00572 
00573   if (!writeEnabled)
00574     last_sorted = key;
00575   return MB_SUCCESS;
00576 }
00577 
00578 #undef umax_2
00579 
00580 #define DIGIT_BITS   8
00581 #define DIGIT_VALUES (1<<DIGIT_BITS)
00582 #define DIGIT_MASK   ((Value)(DIGIT_VALUES-1))
00583 #define CEILDIV(a,b) (((a)+(b)-1)/(b))
00584 #define DIGITS       CEILDIV(CHAR_BIT*sizeof(Value),DIGIT_BITS)
00585 #define VALUE_BITS   (DIGIT_BITS*DIGITS)
00586 #define COUNT_SIZE   (DIGITS*DIGIT_VALUES)
00587 
00588 /* used to unroll a tiny loop: */
00589 #define COUNT_DIGIT_01(n,i) \
00590     if(n>i) count[i][val&DIGIT_MASK]++, val>>=DIGIT_BITS
00591 #define COUNT_DIGIT_02(n,i) COUNT_DIGIT_01(n,i); COUNT_DIGIT_01(n,i+ 1)
00592 #define COUNT_DIGIT_04(n,i) COUNT_DIGIT_02(n,i); COUNT_DIGIT_02(n,i+ 2)
00593 #define COUNT_DIGIT_08(n,i) COUNT_DIGIT_04(n,i); COUNT_DIGIT_04(n,i+ 4)
00594 #define COUNT_DIGIT_16(n,i) COUNT_DIGIT_08(n,i); COUNT_DIGIT_08(n,i+ 8)
00595 #define COUNT_DIGIT_32(n,i) COUNT_DIGIT_16(n,i); COUNT_DIGIT_16(n,i+16)
00596 #define COUNT_DIGIT_64(n,i) COUNT_DIGIT_32(n,i); COUNT_DIGIT_32(n,i+32)
00597 
00598 template<class Value>
00599 Value TupleList::radix_count(const Value *A, const Value *end, Index stride,
00600     Index count[DIGITS][DIGIT_VALUES])
00601 {
00602   Value bitorkey = 0;
00603   memset(count, 0, COUNT_SIZE * sizeof(Index));
00604   do
00605   {
00606     Value val = *A;
00607     bitorkey |= val;
00608     COUNT_DIGIT_64(DIGITS, 0);
00609     // above macro expands to:
00610     //if(DIGITS> 0) count[ 0][val&DIGIT_MASK]++, val>>=DIGIT_BITS;
00611     //if(DIGITS> 1) count[ 1][val&DIGIT_MASK]++, val>>=DIGIT_BITS;
00612     //  ...
00613     //if(DIGITS>63) count[63][val&DIGIT_MASK]++, val>>=DIGIT_BITS;
00614 
00615   } while (A += stride, A != end);
00616   return bitorkey;
00617 }
00618 
00619 #undef COUNT_DIGIT_01
00620 #undef COUNT_DIGIT_02
00621 #undef COUNT_DIGIT_04
00622 #undef COUNT_DIGIT_08
00623 #undef COUNT_DIGIT_16
00624 #undef COUNT_DIGIT_32
00625 #undef COUNT_DIGIT_64
00626 
00627 void TupleList::radix_offsets(Index *c)
00628 {
00629   Index sum = 0, t, *ce = c + DIGIT_VALUES;
00630   do
00631     t = *c, *c++ = sum, sum += t;
00632   while (c != ce);
00633 }
00634 
00635 template<class Value>
00636 unsigned TupleList::radix_zeros(Value bitorkey,
00637     Index count[DIGITS][DIGIT_VALUES], unsigned *shift, Index **offsets)
00638 {
00639   unsigned digits = 0, sh = 0;
00640   Index *c = &count[0][0];
00641   do
00642   {
00643     if (bitorkey & DIGIT_MASK)
00644       *shift++ = sh, *offsets++ = c, ++digits, radix_offsets(c);
00645   } while (bitorkey >>= DIGIT_BITS,sh += DIGIT_BITS,c += DIGIT_VALUES,sh
00646       != VALUE_BITS);
00647   return digits;
00648 }
00649 
00650 template<class Value>
00651 void TupleList::radix_index_pass_b(const Value *A, Index n, Index stride,
00652     unsigned sh, Index *off, SortData<Value> *out)
00653 {
00654   Index i = 0;
00655   do
00656   {
00657     Value v = *A;
00658     SortData<Value> *d = &out[off[(v >> sh) & DIGIT_MASK]++];
00659     d->v = v, d->i = i++;
00660   } while (A += stride, i != n);
00661 }
00662 
00663 template<class Value>
00664 void TupleList::radix_index_pass_m(const SortData<Value> *src,
00665     const SortData<Value> *end, unsigned sh, Index *off, SortData<Value> *out)
00666 {
00667   do
00668   {
00669     SortData<Value> *d = &out[off[(src->v >> sh) & DIGIT_MASK]++];
00670     d->v = src->v, d->i = src->i;
00671   } while (++src != end);
00672 }
00673 
00674 template<class Value>
00675 void TupleList::radix_index_pass_e(const SortData<Value> *src,
00676     const SortData<Value> *end, unsigned sh, Index *off, Index *out)
00677 {
00678   do
00679     out[off[(src->v >> sh) & DIGIT_MASK]++] = src->i;
00680   while (++src != end);
00681 }
00682 
00683 template<class Value>
00684 void TupleList::radix_index_pass_be(const Value *A, Index n, Index stride,
00685     unsigned sh, Index *off, Index *out)
00686 {
00687   Index i = 0;
00688   do
00689     out[off[(*A >> sh) & DIGIT_MASK]++] = i++;
00690   while (A += stride, i != n);
00691 }
00692 
00693 template<class Value>
00694 void TupleList::radix_index_sort(const Value *A, Index n, Index stride,
00695     Index *idx, SortData<Value> *work)
00696 {
00697   Index count[DIGITS][DIGIT_VALUES];
00698   Value bitorkey = radix_count(A, A + n * stride, stride, count);
00699   unsigned shift[DIGITS];
00700   Index *offsets[DIGITS];
00701   unsigned digits = radix_zeros(bitorkey, count, shift, offsets);
00702   if (digits == 0)
00703   {
00704     Index i = 0;
00705     do
00706       *idx++ = i++;
00707     while (i != n);
00708   }
00709   else if (digits == 1)
00710   {
00711     radix_index_pass_be(A, n, stride, shift[0], offsets[0], idx);
00712   }
00713   else
00714   {
00715     SortData<Value> *src, *dst;
00716     unsigned d;
00717     if ((digits & 1) == 0)
00718       dst = work, src = dst + n;
00719     else
00720       src = work, dst = src + n;
00721     radix_index_pass_b(A, n, stride, shift[0], offsets[0], src);
00722     for (d = 1; d != digits - 1; ++d)
00723     {
00724       SortData<Value> *t;
00725       radix_index_pass_m(src, src + n, shift[d], offsets[d], dst);
00726       t = src, src = dst, dst = t;
00727     }
00728     radix_index_pass_e(src, src + n, shift[d], offsets[d], idx);
00729   }
00730 }
00731 
00732 template<class Value>
00733 void TupleList::merge_index_sort(const Value *A, const Index An, Index stride,
00734     Index *idx, SortData<Value> *work)
00735 {
00736   SortData<Value> * const buf[2] = { work + An, work };
00737   Index n = An, base = -n, odd = 0, c = 0, b = 1;
00738   Index i = 0;
00739   for (;;)
00740   {
00741     SortData<Value> *p;
00742     if ((c & 1) == 0)
00743     {
00744       base += n, n += (odd & 1), c |= 1, b ^= 1;
00745       while (n > 3)
00746         odd <<= 1, odd |= (n & 1), n >>= 1, c <<= 1, b ^= 1;
00747     }
00748     else
00749       base -= n - (odd & 1), n <<= 1, n -= (odd & 1), odd >>= 1, c >>= 1;
00750     if (c == 0)
00751       break;
00752     p = buf[b] + base;
00753     if (n == 2)
00754     {
00755       Value v[2];
00756       v[0] = *A, A += stride, v[1] = *A, A += stride;
00757       if (v[1] < v[0])
00758         p[0].v = v[1], p[0].i = i + 1, p[1].v = v[0], p[1].i = i;
00759       else
00760         p[0].v = v[0], p[0].i = i, p[1].v = v[1], p[1].i = i + 1;
00761       i += 2;
00762     }
00763     else if (n == 3)
00764     {
00765       Value v[3];
00766       v[0] = *A, A += stride, v[1] = *A, A += stride, v[2] = *A, A += stride;
00767       if (v[1] < v[0])
00768       {
00769         if (v[2] < v[1])
00770           p[0].v = v[2], p[1].v = v[1], p[2].v = v[0], p[0].i = i + 2, p[1].i =
00771               i + 1, p[2].i = i;
00772         else
00773         {
00774           if (v[2] < v[0])
00775             p[0].v = v[1], p[1].v = v[2], p[2].v = v[0], p[0].i = i + 1, p[1].i =
00776                 i + 2, p[2].i = i;
00777           else
00778             p[0].v = v[1], p[1].v = v[0], p[2].v = v[2], p[0].i = i + 1, p[1].i =
00779                 i, p[2].i = i + 2;
00780         }
00781       }
00782       else
00783       {
00784         if (v[2] < v[0])
00785           p[0].v = v[2], p[1].v = v[0], p[2].v = v[1], p[0].i = i + 2, p[1].i =
00786               i, p[2].i = i + 1;
00787         else
00788         {
00789           if (v[2] < v[1])
00790             p[0].v = v[0], p[1].v = v[2], p[2].v = v[1], p[0].i = i, p[1].i = i
00791                 + 2, p[2].i = i + 1;
00792           else
00793             p[0].v = v[0], p[1].v = v[1], p[2].v = v[2], p[0].i = i, p[1].i = i
00794                 + 1, p[2].i = i + 2;
00795         }
00796       }
00797       i += 3;
00798     }
00799     else
00800     {
00801       const Index na = n >> 1, nb = (n + 1) >> 1;
00802       const SortData<Value> *ap = buf[b ^ 1] + base, *ae = ap + na;
00803       SortData<Value> *bp = p + na, *be = bp + nb;
00804       for (;;)
00805       {
00806         if (bp->v < ap->v)
00807         {
00808           *p++ = *bp++;
00809           if (bp != be)
00810             continue;
00811           do
00812             *p++ = *ap++;
00813           while (ap != ae);
00814           break;
00815         }
00816         else
00817         {
00818           *p++ = *ap++;
00819           if (ap == ae)
00820             break;
00821         }
00822       }
00823     }
00824   }
00825   {
00826     const SortData<Value> *p = buf[0], *pe = p + An;
00827     do
00828       *idx++ = (p++)->i;
00829     while (p != pe);
00830   }
00831 }
00832 
00833 template<class Value>
00834 void TupleList::index_sort(const Value *A, Index n, Index stride, Index *idx,
00835     SortData<Value> *work)
00836 {
00837   if (n < DIGIT_VALUES)
00838   {
00839     if (n == 0)
00840       return;
00841     if (n == 1)
00842       *idx = 0;
00843     else
00844       merge_index_sort(A, n, stride, idx, work);
00845   }
00846   else
00847     radix_index_sort(A, n, stride, idx, work);
00848 }
00849 
00850 #undef DIGIT_BITS
00851 #undef DIGIT_VALUES
00852 #undef DIGIT_MASK
00853 #undef CEILDIV
00854 #undef DIGITS
00855 #undef VALUE_BITS
00856 #undef COUNT_SIZE
00857 #undef sort_data_long
00858 
00859 } //namespace
00860 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines