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