Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
TupleList.cpp
Go to the documentation of this file.
1 #include <cstring>
2 #include <climits>
3 #include <cstdio>
4 #include <cstdlib>
5 #include <cstdarg>
6 #include <iostream>
7 #include <fstream>
8 
9 #include "moab/TupleList.hpp"
10 
11 namespace moab
12 {
13 
14 void fail( const char* fmt, ... )
15 {
16  va_list ap;
17  va_start( ap, fmt );
18  vfprintf( stderr, fmt, ap );
19  va_end( ap );
20  exit( 1 );
21 }
22 
24 {
25  ptr = NULL;
26  buffSize = 0;
27  this->buffer_init_( sz, __FILE__ );
28 }
29 
31 {
32  buffSize = 0;
33  ptr = NULL;
34 }
35 
36 void TupleList::buffer::buffer_init_( size_t sizeIn, const char* file )
37 {
38  this->buffSize = sizeIn;
39  void* res = malloc( this->buffSize );
40  if( !res && buffSize > 0 ) fail( "%s: allocation of %d bytes failed\n", file, (int)buffSize );
41  ptr = (char*)res;
42 }
43 
44 void TupleList::buffer::buffer_reserve_( size_t min, const char* file )
45 {
46  if( this->buffSize < min )
47  {
48  size_t newSize = this->buffSize;
49  newSize += newSize / 2 + 1;
50  if( newSize < min ) newSize = min;
51  void* res = realloc( ptr, newSize );
52  if( !res && newSize > 0 ) fail( "%s: reallocation of %d bytes failed\n", file, newSize );
53  ptr = (char*)res;
54  this->buffSize = newSize;
55  }
56 }
57 
59 {
60  free( ptr );
61  ptr = NULL;
62  buffSize = 0;
63 }
64 
65 TupleList::TupleList( uint p_mi, uint p_ml, uint p_mul, uint p_mr, uint p_max )
66  : vi( NULL ), vl( NULL ), vul( NULL ), vr( NULL ), last_sorted( -1 )
67 {
68  initialize( p_mi, p_ml, p_mul, p_mr, p_max );
69 }
70 
72  : vi_rd( NULL ), vl_rd( NULL ), vul_rd( NULL ), vr_rd( NULL ), mi( 0 ), ml( 0 ), mul( 0 ), mr( 0 ), n( 0 ),
73  max( 0 ), vi( NULL ), vl( NULL ), vul( NULL ), vr( NULL ), last_sorted( -1 )
74 {
76 }
77 
78 // Allocates space for the tuple list in memory according to parameters
79 void TupleList::initialize( uint p_mi, uint p_ml, uint p_mul, uint p_mr, uint p_max )
80 {
81  this->n = 0;
82  this->max = p_max;
83  this->mi = p_mi;
84  this->ml = p_ml;
85  this->mul = p_mul;
86  this->mr = p_mr;
87  size_t sz;
88 
89  if( max * mi > 0 )
90  {
91  sz = max * mi * sizeof( sint );
92  void* resi = malloc( sz );
93  if( !resi && max * mi > 0 ) fail( "%s: allocation of %d bytes failed\n", __FILE__, (int)sz );
94  vi = (sint*)resi;
95  }
96  else
97  vi = NULL;
98  if( max * ml > 0 )
99  {
100  sz = max * ml * sizeof( slong );
101  void* resl = malloc( sz );
102  if( !resl && max * ml > 0 ) fail( "%s: allocation of %d bytes failed\n", __FILE__, (int)sz );
103  vl = (slong*)resl;
104  }
105  else
106  vl = NULL;
107  if( max * mul > 0 )
108  {
109  sz = max * mul * sizeof( Ulong );
110  void* resu = malloc( sz );
111  if( !resu && max * mul > 0 ) fail( "%s: allocation of %d bytes failed\n", __FILE__, (int)sz );
112  vul = (Ulong*)resu;
113  }
114  else
115  vul = NULL;
116  if( max * mr > 0 )
117  {
118  sz = max * mr * sizeof( realType );
119  void* resr = malloc( sz );
120  if( !resr && max * ml > 0 ) fail( "%s: allocation of %d bytes failed\n", __FILE__, (int)sz );
121  vr = (realType*)resr;
122  }
123  else
124  vr = NULL;
125 
126  // Begin with write access disabled
127  this->disableWriteAccess();
128 
129  // Set read variables
130  vi_rd = vi;
131  vl_rd = vl;
132  vul_rd = vul;
133  vr_rd = vr;
134 }
135 
136 // Resizes a tuplelist to the given uint max
138 {
139  this->max = maxIn;
140  size_t sz;
141 
142  if( vi || ( max * mi > 0 ) )
143  {
144  sz = max * mi * sizeof( sint );
145  void* resi = realloc( vi, sz );
146  if( !resi && max * mi > 0 )
147  {
148  fail( "%s: allocation of %d bytes failed\n", __FILE__, (int)sz );
150  }
151  vi = (sint*)resi;
152  }
153  if( vl || ( max * ml > 0 ) )
154  {
155  sz = max * ml * sizeof( slong );
156  void* resl = realloc( vl, sz );
157  if( !resl && max * ml > 0 )
158  {
159  fail( "%s: allocation of %d bytes failed\n", __FILE__, (int)sz );
161  }
162  vl = (slong*)resl;
163  }
164  if( vul || ( max * mul > 0 ) )
165  {
166  sz = max * mul * sizeof( Ulong );
167  void* resu = realloc( vul, sz );
168  if( !resu && max * mul > 0 )
169  {
170  fail( "%s: allocation of %d bytes failed\n", __FILE__, (int)sz );
172  }
173  vul = (Ulong*)resu;
174  }
175  if( vr || ( max * mr > 0 ) )
176  {
177  sz = max * mr * sizeof( realType );
178  void* resr = realloc( vr, sz );
179  if( !resr && max * mr > 0 )
180  {
181  fail( "%s: allocation of %d bytes failed\n", __FILE__, (int)sz );
183  }
184  vr = (realType*)resr;
185  }
186 
187  // Set read variables
188  vi_rd = vi;
189  vl_rd = vl;
190  vul_rd = vul;
191  vr_rd = vr;
192 
193  // Set the write variables if necessary
194  if( writeEnabled )
195  {
196  vi_wr = vi;
197  vl_wr = vl;
198  vul_wr = vul;
199  vr_wr = vr;
200  }
201  return moab::MB_SUCCESS;
202 }
203 
204 // Frees the memory used by the tuplelist
206 {
207  // free up the pointers
208  free( vi );
209  free( vl );
210  free( vul );
211  free( vr );
212  // Set them all to null
213  vr = NULL;
214  vi = NULL;
215  vul = NULL;
216  vl = NULL;
217  // Set the read and write pointers to null
219  vi_rd = NULL;
220  vl_rd = NULL;
221  vul_rd = NULL;
222  vr_rd = NULL;
223 }
224 
225 // Increments n; if n>max, increase the size of the tuplelist
227 {
228  n++;
229  while( n > max )
230  resize( ( max ? max + max / 2 + 1 : 2 ) );
231  last_sorted = -1;
232 }
233 
234 // Given the value and the position in the field, finds the index of the tuple
235 // to which the value belongs
236 int TupleList::find( unsigned int key_num, sint value )
237 {
238  // we are passing an int, no issue, leave it at long
239  long uvalue = (long)value;
240  if( !( key_num > mi ) )
241  {
242  // Binary search: only if the tuple_list is sorted
243  if( last_sorted == (int)key_num )
244  {
245  int lb = 0, ub = n, index; // lb=lower bound, ub=upper bound, index=mid
246  for( ; lb <= ub; )
247  {
248  index = ( lb + ub ) / 2;
249  if( vi[index * mi + key_num] == uvalue )
250  return index;
251  else if( vi[index * mi + key_num] > uvalue )
252  ub = index - 1;
253  else if( vi[index * mi + key_num] < uvalue )
254  lb = index + 1;
255  }
256  }
257  else
258  {
259  // Sequential search: if tuple_list is not sorted
260  for( uint index = 0; index < n; index++ )
261  {
262  if( vi[index * mi + key_num] == uvalue ) return index;
263  }
264  }
265  }
266  return -1; // If the value wasn't present or an invalid key was given
267 }
268 
269 int TupleList::find( unsigned int key_num, slong value )
270 {
271  long uvalue = (long)value;
272  if( !( key_num > ml ) )
273  {
274  if( last_sorted - mi == key_num )
275  {
276  int lb = 0, ub = n, index; // lb=lower bound, ub=upper bound, index=mid
277  for( ; lb <= ub; )
278  {
279  index = ( lb + ub ) / 2;
280  if( vl[index * ml + key_num] == uvalue )
281  return index;
282  else if( vl[index * ml + key_num] > uvalue )
283  ub = index - 1;
284  else if( vl[index * ml + key_num] < uvalue )
285  lb = index + 1;
286  }
287  }
288  else
289  {
290  // Sequential search: if tuple_list is not sorted
291  for( uint index = 0; index < n; index++ )
292  {
293  if( vl[index * ml + key_num] == uvalue ) return index;
294  }
295  }
296  }
297  return -1; // If the value wasn't present or an invalid key was given
298 }
299 
300 int TupleList::find( unsigned int key_num, Ulong value )
301 {
302  if( !( key_num > mul ) )
303  {
304  if( last_sorted - mi - ml == key_num )
305  {
306  int lb = 0, ub = n - 1, index; // lb=lower bound, ub=upper bound, index=mid
307  for( ; lb <= ub; )
308  {
309  index = ( lb + ub ) / 2;
310  if( vul[index * mul + key_num] == value )
311  return index;
312  else if( vul[index * mul + key_num] > value )
313  ub = index - 1;
314  else if( vul[index * mul + key_num] < value )
315  lb = index + 1;
316  }
317  }
318  else
319  {
320  // Sequential search: if tuple_list is not sorted
321  for( uint index = 0; index < n; index++ )
322  {
323  if( vul[index * mul + key_num] == value ) return index;
324  }
325  }
326  }
327  return -1; // If the value wasn't present or an invalid key was given
328 }
329 
330 int TupleList::find( unsigned int key_num, realType value )
331 {
332  if( !( key_num > mr ) )
333  {
334  // Sequential search: TupleList cannot be sorted by reals
335  for( uint index = 0; index < n; index++ )
336  {
337  if( vr[index * mr + key_num] == value ) return index;
338  }
339  }
340  return -1; // If the value wasn't present or an invalid key was given
341 }
342 
343 sint TupleList::get_sint( unsigned int index, unsigned int m )
344 {
345  if( mi > m && n > index ) return vi[index * mi + m];
346  return 0;
347 }
348 
349 slong TupleList::get_int( unsigned int index, unsigned int m )
350 {
351  if( ml > m && n > index ) return vl[index * ml + m];
352  return 0;
353 }
354 
355 Ulong TupleList::get_ulong( unsigned int index, unsigned int m )
356 {
357  if( mul > m && n > index ) return vul[index * mul + m];
358  return 0;
359 }
360 
361 realType TupleList::get_double( unsigned int index, unsigned int m )
362 {
363  if( mr > m && n > index ) return vr[index * mr + m];
364  return 0;
365 }
366 
367 ErrorCode TupleList::get( unsigned int index, const sint*& sp, const slong*& ip, const Ulong*& lp, const realType*& dp )
368 {
369  if( index <= n )
370  {
371  if( mi )
372  *&sp = &vi[index * mi];
373  else
374  *&sp = NULL;
375  if( ml )
376  *&ip = &vl[index * ml];
377  else
378  *&ip = NULL;
379  if( mul )
380  *&lp = &vul[index * mul];
381  else
382  *&lp = NULL;
383  if( mr )
384  *&dp = &vr[index * mr];
385  else
386  *&dp = NULL;
387 
388  return MB_SUCCESS;
389  }
390  return MB_FAILURE;
391 }
392 
393 unsigned int TupleList::push_back( sint* sp, slong* ip, Ulong* lp, realType* dp )
394 {
395  reserve();
396  if( mi ) memcpy( &vi[mi * ( n - 1 )], sp, mi * sizeof( sint ) );
397  if( ml ) memcpy( &vl[ml * ( n - 1 )], ip, ml * sizeof( long ) );
398  if( mul ) memcpy( &vul[mul * ( n - 1 )], lp, mul * sizeof( Ulong ) );
399  if( mr ) memcpy( &vr[mr * ( n - 1 )], dp, mr * sizeof( realType ) );
400 
401  last_sorted = -1;
402  return n - 1;
403 }
404 
406 {
407  writeEnabled = true;
408  last_sorted = -1;
409  vi_wr = vi;
410  vl_wr = vl;
411  vul_wr = vul;
412  vr_wr = vr;
413 }
414 
416 {
417  writeEnabled = false;
418  vi_wr = NULL;
419  vl_wr = NULL;
420  vul_wr = NULL;
421  vr_wr = NULL;
422 }
423 
424 void TupleList::getTupleSize( uint& mi_out, uint& ml_out, uint& mul_out, uint& mr_out ) const
425 {
426  mi_out = mi;
427  ml_out = ml;
428  mul_out = mul;
429  mr_out = mr;
430 }
431 
433 {
434  // Check for direct write access
435  if( !writeEnabled )
436  {
438  }
439  n++;
440  return n;
441 }
442 
443 void TupleList::set_n( uint n_in )
444 {
445  // Check for direct write access;
446  if( !writeEnabled )
447  {
449  }
450  n = n_in;
451 }
452 
453 void TupleList::print( const char* name ) const
454 {
455  std::cout << "Printing Tuple " << name << "===================" << std::endl;
456  unsigned long i = 0, l = 0, ul = 0, r = 0;
457  for( uint k = 0; k < n; k++ )
458  {
459  for( uint j = 0; j < mi; j++ )
460  {
461  std::cout << vi[i++] << " | ";
462  }
463  for( uint j = 0; j < ml; j++ )
464  {
465  std::cout << vl[l++] << " | ";
466  }
467  for( uint j = 0; j < mul; j++ )
468  {
469  std::cout << vul[ul++] << " | ";
470  }
471  for( uint j = 0; j < mr; j++ )
472  {
473  std::cout << vr[r++] << " | ";
474  }
475  std::cout << std::endl;
476  }
477  std::cout << "=======================================" << std::endl << std::endl;
478 }
479 void TupleList::print_to_file( const char* filename ) const
480 {
481  std::ofstream ofs;
482  ofs.open( filename, std::ofstream::out | std::ofstream::app );
483 
484  ofs << "Printing Tuple " << filename << "===================" << std::endl;
485  unsigned long i = 0, l = 0, ul = 0, r = 0;
486  for( uint k = 0; k < n; k++ )
487  {
488  for( uint j = 0; j < mi; j++ )
489  {
490  ofs << vi[i++] << " | ";
491  }
492  for( uint j = 0; j < ml; j++ )
493  {
494  ofs << vl[l++] << " | ";
495  }
496  for( uint j = 0; j < mul; j++ )
497  {
498  ofs << vul[ul++] << " | ";
499  }
500  for( uint j = 0; j < mr; j++ )
501  {
502  ofs << vr[r++] << " | ";
503  }
504  ofs << std::endl;
505  }
506  ofs << "=======================================" << std::endl << std::endl;
507 
508  ofs.close();
509 }
510 void TupleList::permute( uint* perm, void* work )
511 {
512  const unsigned int_size = mi * sizeof( sint ), long_size = ml * sizeof( slong ), Ulong_size = mul * sizeof( Ulong ),
513  real_size = mr * sizeof( realType );
514  if( mi )
515  {
516  uint *p = perm, *pe = p + n;
517  char* sorted = (char*)work;
518  while( p != pe )
519  memcpy( (void*)sorted, &vi[mi * ( *p++ )], int_size ), sorted += int_size;
520  memcpy( vi, work, int_size * n );
521  }
522  if( ml )
523  {
524  uint *p = perm, *pe = p + n;
525  char* sorted = (char*)work;
526  while( p != pe )
527  memcpy( (void*)sorted, &vl[ml * ( *p++ )], long_size ), sorted += long_size;
528  memcpy( vl, work, long_size * n );
529  }
530  if( mul )
531  {
532  uint *p = perm, *pe = p + n;
533  char* sorted = (char*)work;
534  while( p != pe )
535  memcpy( (void*)sorted, &vul[mul * ( *p++ )], Ulong_size ), sorted += Ulong_size;
536  memcpy( vul, work, Ulong_size * n );
537  }
538  if( mr )
539  {
540  uint *p = perm, *pe = p + n;
541  char* sorted = (char*)work;
542  while( p != pe )
543  memcpy( (void*)sorted, &vr[mr * ( *p++ )], real_size ), sorted += real_size;
544  memcpy( vr, work, real_size * n );
545  }
546 }
547 
548 #define umax_2( a, b ) ( ( ( a ) > ( b ) ) ? ( a ) : ( b ) )
549 
551 {
552  const unsigned int_size = mi * sizeof( sint );
553  const unsigned long_size = ml * sizeof( slong );
554  const unsigned Ulong_size = mul * sizeof( Ulong );
555  const unsigned real_size = mr * sizeof( realType );
556  const unsigned width = umax_2( umax_2( int_size, long_size ), umax_2( Ulong_size, real_size ) );
557  unsigned data_size = key >= mi ? sizeof( SortData< long > ) : sizeof( SortData< uint > );
558 #if defined( WIN32 ) || defined( _WIN32 )
559  if( key >= mi + ml ) data_size = sizeof( SortData< Ulong > );
560 #endif
561 
562  uint work_min = n * umax_2( 2 * data_size, sizeof( sint ) + width );
563  uint* work;
564  buf->buffer_reserve( work_min );
565  work = (uint*)buf->ptr;
566  if( key < mi )
567  index_sort( (uint*)&vi[key], n, mi, work, (SortData< uint >*)work );
568  else if( key < mi + ml )
569  index_sort( (long*)&vl[key - mi], n, ml, work, (SortData< long >*)work );
570  else if( key < mi + ml + mul )
571  index_sort( (Ulong*)&vul[key - mi - ml], n, mul, work, (SortData< Ulong >*)work );
572  else
573  return MB_NOT_IMPLEMENTED;
574 
575  permute( work, work + n );
576 
577  if( !writeEnabled ) last_sorted = key;
578  return MB_SUCCESS;
579 }
580 
581 #undef umax_2
582 
583 #define DIGIT_BITS 8
584 #define DIGIT_VALUES ( 1 << DIGIT_BITS )
585 #define DIGIT_MASK ( (Value)( DIGIT_VALUES - 1 ) )
586 #define CEILDIV( a, b ) ( ( ( a ) + (b)-1 ) / ( b ) )
587 #define DIGITS CEILDIV( CHAR_BIT * sizeof( Value ), DIGIT_BITS )
588 #define VALUE_BITS ( DIGIT_BITS * DIGITS )
589 #define COUNT_SIZE ( DIGITS * DIGIT_VALUES )
590 
591 /* used to unroll a tiny loop: */
592 #define COUNT_DIGIT_01( n, i ) \
593  if( ( n ) > ( i ) ) count[i][val & DIGIT_MASK]++, val >>= DIGIT_BITS
594 #define COUNT_DIGIT_02( n, i ) \
595  COUNT_DIGIT_01( n, i ); \
596  COUNT_DIGIT_01( n, ( i ) + 1 )
597 #define COUNT_DIGIT_04( n, i ) \
598  COUNT_DIGIT_02( n, i ); \
599  COUNT_DIGIT_02( n, ( i ) + 2 )
600 #define COUNT_DIGIT_08( n, i ) \
601  COUNT_DIGIT_04( n, i ); \
602  COUNT_DIGIT_04( n, ( i ) + 4 )
603 #define COUNT_DIGIT_16( n, i ) \
604  COUNT_DIGIT_08( n, i ); \
605  COUNT_DIGIT_08( n, ( i ) + 8 )
606 #define COUNT_DIGIT_32( n, i ) \
607  COUNT_DIGIT_16( n, i ); \
608  COUNT_DIGIT_16( n, ( i ) + 16 )
609 #define COUNT_DIGIT_64( n, i ) \
610  COUNT_DIGIT_32( n, i ); \
611  COUNT_DIGIT_32( n, ( i ) + 32 )
612 
613 template < class Value >
614 Value TupleList::radix_count( const Value* A, const Value* end, Index stride, Index count[DIGITS][DIGIT_VALUES] )
615 {
616  Value bitorkey = 0;
617  memset( count, 0, COUNT_SIZE * sizeof( Index ) );
618  do
619  {
620  Value val = *A;
621  bitorkey |= val;
622  COUNT_DIGIT_64( DIGITS, 0 );
623  // above macro expands to:
624  // if(DIGITS> 0) count[ 0][val&DIGIT_MASK]++, val>>=DIGIT_BITS;
625  // if(DIGITS> 1) count[ 1][val&DIGIT_MASK]++, val>>=DIGIT_BITS;
626  // ...
627  // if(DIGITS>63) count[63][val&DIGIT_MASK]++, val>>=DIGIT_BITS;
628 
629  } while( A += stride, A != end );
630  return bitorkey;
631 }
632 
633 #undef COUNT_DIGIT_01
634 #undef COUNT_DIGIT_02
635 #undef COUNT_DIGIT_04
636 #undef COUNT_DIGIT_08
637 #undef COUNT_DIGIT_16
638 #undef COUNT_DIGIT_32
639 #undef COUNT_DIGIT_64
640 
642 {
643  Index sum = 0, t, *ce = c + DIGIT_VALUES;
644  do
645  t = *c, *c++ = sum, sum += t;
646  while( c != ce );
647 }
648 
649 template < class Value >
650 unsigned TupleList::radix_zeros( Value bitorkey, Index count[DIGITS][DIGIT_VALUES], unsigned* shift, Index** offsets )
651 {
652  unsigned digits = 0, sh = 0;
653  Index* c = &count[0][0];
654  do
655  {
656  if( bitorkey & DIGIT_MASK ) *shift++ = sh, *offsets++ = c, ++digits, radix_offsets( c );
657  } while( bitorkey >>= DIGIT_BITS, sh += DIGIT_BITS, c += DIGIT_VALUES, sh != VALUE_BITS );
658  return digits;
659 }
660 
661 template < class Value >
662 void TupleList::radix_index_pass_b( const Value* A,
663  Index n,
664  Index stride,
665  unsigned sh,
666  Index* off,
667  SortData< Value >* out )
668 {
669  Index i = 0;
670  do
671  {
672  Value v = *A;
673  SortData< Value >* d = &out[off[( v >> sh ) & DIGIT_MASK]++];
674  d->v = v, d->i = i++;
675  } while( A += stride, i != n );
676 }
677 
678 template < class Value >
680  const SortData< Value >* end,
681  unsigned sh,
682  Index* off,
683  SortData< Value >* out )
684 {
685  do
686  {
687  SortData< Value >* d = &out[off[( src->v >> sh ) & DIGIT_MASK]++];
688  d->v = src->v, d->i = src->i;
689  } while( ++src != end );
690 }
691 
692 template < class Value >
694  const SortData< Value >* end,
695  unsigned sh,
696  Index* off,
697  Index* out )
698 {
699  do
700  out[off[( src->v >> sh ) & DIGIT_MASK]++] = src->i;
701  while( ++src != end );
702 }
703 
704 template < class Value >
705 void TupleList::radix_index_pass_be( const Value* A, Index n, Index stride, unsigned sh, Index* off, Index* out )
706 {
707  Index i = 0;
708  do
709  out[off[( *A >> sh ) & DIGIT_MASK]++] = i++;
710  while( A += stride, i != n );
711 }
712 
713 template < class Value >
714 void TupleList::radix_index_sort( const Value* A, Index n, Index stride, Index* idx, SortData< Value >* work )
715 {
716  Index count[DIGITS][DIGIT_VALUES];
717  Value bitorkey = radix_count( A, A + n * stride, stride, count );
718  unsigned shift[DIGITS];
719  Index* offsets[DIGITS];
720  unsigned digits = radix_zeros( bitorkey, count, shift, offsets );
721  if( digits == 0 )
722  {
723  Index i = 0;
724  do
725  *idx++ = i++;
726  while( i != n );
727  }
728  else if( digits == 1 )
729  {
730  radix_index_pass_be( A, n, stride, shift[0], offsets[0], idx );
731  }
732  else
733  {
734  SortData< Value >*src, *dst;
735  unsigned d;
736  if( ( digits & 1 ) == 0 )
737  dst = work, src = dst + n;
738  else
739  src = work, dst = src + n;
740  radix_index_pass_b( A, n, stride, shift[0], offsets[0], src );
741  for( d = 1; d != digits - 1; ++d )
742  {
744  radix_index_pass_m( src, src + n, shift[d], offsets[d], dst );
745  t = src, src = dst, dst = t;
746  }
747  radix_index_pass_e( src, src + n, shift[d], offsets[d], idx );
748  }
749 }
750 
751 template < class Value >
752 void TupleList::merge_index_sort( const Value* A, const Index An, Index stride, Index* idx, SortData< Value >* work )
753 {
754  SortData< Value >* const buf[2] = { work + An, work };
755  Index n = An, base = -n, odd = 0, c = 0, b = 1;
756  Index i = 0;
757  for( ;; )
758  {
760  if( ( c & 1 ) == 0 )
761  {
762  base += n, n += ( odd & 1 ), c |= 1, b ^= 1;
763  while( n > 3 )
764  odd <<= 1, odd |= ( n & 1 ), n >>= 1, c <<= 1, b ^= 1;
765  }
766  else
767  base -= n - ( odd & 1 ), n <<= 1, n -= ( odd & 1 ), odd >>= 1, c >>= 1;
768  if( c == 0 ) break;
769  p = buf[b] + base;
770  if( n == 2 )
771  {
772  Value v[2];
773  v[0] = *A, A += stride, v[1] = *A, A += stride;
774  if( v[1] < v[0] )
775  p[0].v = v[1], p[0].i = i + 1, p[1].v = v[0], p[1].i = i;
776  else
777  p[0].v = v[0], p[0].i = i, p[1].v = v[1], p[1].i = i + 1;
778  i += 2;
779  }
780  else if( n == 3 )
781  {
782  Value v[3];
783  v[0] = *A, A += stride, v[1] = *A, A += stride, v[2] = *A, A += stride;
784  if( v[1] < v[0] )
785  {
786  if( v[2] < v[1] )
787  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;
788  else
789  {
790  if( v[2] < v[0] )
791  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;
792  else
793  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;
794  }
795  }
796  else
797  {
798  if( v[2] < v[0] )
799  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;
800  else
801  {
802  if( v[2] < v[1] )
803  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;
804  else
805  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;
806  }
807  }
808  i += 3;
809  }
810  else
811  {
812  const Index na = n >> 1, nb = ( n + 1 ) >> 1;
813  const SortData< Value >*ap = buf[b ^ 1] + base, *ae = ap + na;
814  SortData< Value >*bp = p + na, *be = bp + nb;
815  for( ;; )
816  {
817  if( bp->v < ap->v )
818  {
819  *p++ = *bp++;
820  if( bp != be ) continue;
821  do
822  *p++ = *ap++;
823  while( ap != ae );
824  break;
825  }
826  else
827  {
828  *p++ = *ap++;
829  if( ap == ae ) break;
830  }
831  }
832  }
833  }
834  {
835  const SortData< Value >*p = buf[0], *pe = p + An;
836  do
837  *idx++ = ( p++ )->i;
838  while( p != pe );
839  }
840 }
841 
842 template < class Value >
843 void TupleList::index_sort( const Value* A, Index n, Index stride, Index* idx, SortData< Value >* work )
844 {
845  if( n < DIGIT_VALUES )
846  {
847  if( n == 0 ) return;
848  if( n == 1 )
849  *idx = 0;
850  else
851  merge_index_sort( A, n, stride, idx, work );
852  }
853  else
854  radix_index_sort( A, n, stride, idx, work );
855 }
856 
857 #undef DIGIT_BITS
858 #undef DIGIT_VALUES
859 #undef DIGIT_MASK
860 #undef CEILDIV
861 #undef DIGITS
862 #undef VALUE_BITS
863 #undef COUNT_SIZE
864 #undef sort_data_long
865 
866 } // namespace moab