 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
23 TupleList::buffer::buffer( size_t sz )
24 {
25 ptr = NULL;
26 buffSize = 0;
27 this->buffer_init_( sz, __FILE__ );
28 }
29
30 TupleList::buffer::buffer()
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
58 void TupleList::buffer::reset()
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
71 TupleList::TupleList()
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 {
75 disableWriteAccess();
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
137 ErrorCode TupleList::resize( uint maxIn )
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 );
149 return moab::MB_MEMORY_ALLOCATION_FAILED;
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 );
160 return moab::MB_MEMORY_ALLOCATION_FAILED;
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 );
171 return moab::MB_MEMORY_ALLOCATION_FAILED;
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 );
182 return moab::MB_MEMORY_ALLOCATION_FAILED;
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
205 void TupleList::reset()
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
218 disableWriteAccess();
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
226 void TupleList::reserve()
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
405 void TupleList::enableWriteAccess()
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
415 void TupleList::disableWriteAccess()
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
432 uint TupleList::inc_n()
433 {
434 // Check for direct write access
435 if( !writeEnabled )
436 {
437 enableWriteAccess();
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 {
448 enableWriteAccess();
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
550 ErrorCode TupleList::sort( uint key, TupleList::buffer* buf )
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
641 void TupleList::radix_offsets( Index* c )
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 >
679 void TupleList::radix_index_pass_m( const SortData< Value >* src,
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 >
693 void TupleList::radix_index_pass_e( const SortData< Value >* src,
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 {
743 SortData< Value >* t;
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 {
759 SortData< Value >* p;
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
1.9.1.