moab
|
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