Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ParallelMergeMesh.cpp
Go to the documentation of this file.
2 #include "moab/Core.hpp"
3 #include "moab/CartVect.hpp"
4 #include "moab/BoundBox.hpp"
5 #include "moab/Skinner.hpp"
6 #include "moab/MergeMesh.hpp"
7 #include "moab/CN.hpp"
8 #include <cfloat>
9 #include <algorithm>
10 
11 #ifdef MOAB_HAVE_MPI
12 #include "moab_mpi.h"
13 #endif
14 
15 namespace moab
16 {
17 
18 // Constructor
19 /*Get Merge Data and tolerance*/
20 ParallelMergeMesh::ParallelMergeMesh( ParallelComm* pc, const double epsilon ) : myPcomm( pc ), myEps( epsilon )
21 {
22  myMB = pc->get_moab();
23  mySkinEnts.resize( 4 );
24 }
25 
26 // Have a wrapper function on the actual merge to avoid memory leaks
27 // Merges elements within a proximity of epsilon
28 ErrorCode ParallelMergeMesh::merge( EntityHandle levelset, bool skip_local_merge, int dim )
29 {
30  ErrorCode rval = PerformMerge( levelset, skip_local_merge, dim );MB_CHK_ERR( rval );
31  CleanUp();
32  return rval;
33 }
34 
35 // Perform the merge
36 ErrorCode ParallelMergeMesh::PerformMerge( EntityHandle levelset, bool skip_local_merge, int dim )
37 {
38  // Get the mesh dimension
39  ErrorCode rval;
40  if( dim < 0 )
41  {
42  rval = myMB->get_dimension( dim );MB_CHK_ERR( rval );
43  }
44 
45  // Get the local skin elements
46  rval = PopulateMySkinEnts( levelset, dim, skip_local_merge );
47  // If there is only 1 proc, we can return now
48  if( rval != MB_SUCCESS || myPcomm->size() == 1 )
49  {
50  return rval;
51  }
52 
53  // Determine the global bounding box
54  double gbox[6];
55  rval = GetGlobalBox( gbox );MB_CHK_ERR( rval );
56 
57  /* Assemble The Destination Tuples */
58  // Get a list of tuples which contain (toProc, handle, x,y,z)
59  myTup.initialize( 1, 0, 1, 3, mySkinEnts[0].size() );
60  rval = PopulateMyTup( gbox );MB_CHK_ERR( rval );
61 
62  /* Gather-Scatter Tuple
63  -tup comes out as (remoteProc,handle,x,y,z) */
64  myCD.initialize( myPcomm->comm() );
65 
66  // 1 represents dynamic tuple, 0 represents index of the processor to send to
67  myCD.gs_transfer( 1, myTup, 0 );
68 
69  /* Sort By X,Y,Z
70  -Utilizes a custom quick sort incorporating epsilon*/
72 
73  // Initialize another tuple list for matches
74  myMatches.initialize( 2, 0, 2, 0, mySkinEnts[0].size() );
75 
76  // ID the matching tuples
77  rval = PopulateMyMatches();MB_CHK_ERR( rval );
78 
79  // We can free up the tuple myTup now
80  myTup.reset();
81 
82  /*Gather-Scatter Again*/
83  // 1 represents dynamic list, 0 represents proc index to send tuple to
84  myCD.gs_transfer( 1, myMatches, 0 );
85  // We can free up the crystal router now
86  myCD.reset();
87 
88  // Sort the matches tuple list
89  SortMyMatches();
90 
91  // Tag the shared elements
92  rval = TagSharedElements( dim );MB_CHK_ERR( rval );
93 
94  // Free up the matches tuples
95  myMatches.reset();
96  return rval;
97 }
98 
99 // Sets mySkinEnts with all of the skin entities on the processor
100 ErrorCode ParallelMergeMesh::PopulateMySkinEnts( const EntityHandle meshset, int dim, bool skip_local_merge )
101 {
102  /*Merge Mesh Locally*/
103  // Get all dim dimensional entities
104  Range ents;
105  ErrorCode rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval );
106 
107  if( ents.empty() && dim == 3 )
108  {
109  dim--;
110  rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval ); // maybe dimension 2
111  }
112 
113  // Merge Mesh Locally
114  if( !skip_local_merge )
115  {
116  MergeMesh merger( myMB, false );
117  merger.merge_entities( ents, myEps );
118  // We can return if there is only 1 proc
119  if( rval != MB_SUCCESS || myPcomm->size() == 1 )
120  {
121  return rval;
122  }
123 
124  // Rebuild the ents range
125  ents.clear();
126  rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval );
127  }
128 
129  /*Get Skin
130  -Get Range of all dimensional entities
131  -skinEnts[i] is the skin entities of dimension i*/
132  Skinner skinner( myMB );
133  for( int skin_dim = dim; skin_dim >= 0; skin_dim-- )
134  {
135  rval = skinner.find_skin( meshset, ents, skin_dim, mySkinEnts[skin_dim] );MB_CHK_ERR( rval );
136  }
137  return MB_SUCCESS;
138 }
139 
140 // Determine the global assembly box
142 {
143  ErrorCode rval;
144 
145  /*Get Bounding Box*/
146  BoundBox box;
147  if( mySkinEnts[0].size() != 0 )
148  {
149  rval = box.update( *myMB, mySkinEnts[0] );MB_CHK_ERR( rval );
150  }
151 
152  // Invert the max
153  box.bMax *= -1;
154 
155  /*Communicate to all processors*/
156  MPI_Allreduce( (void*)&box, gbox, 6, MPI_DOUBLE, MPI_MIN, myPcomm->comm() );
157 
158  /*Assemble Global Bounding Box*/
159  // Flip the max back
160  for( int i = 3; i < 6; i++ )
161  {
162  gbox[i] *= -1;
163  }
164  return MB_SUCCESS;
165 }
166 
167 // Assemble the tuples with their processor destination
169 {
170  /*Figure out how do partition the global box*/
171  double lengths[3];
172  int parts[3];
173  ErrorCode rval = PartitionGlobalBox( gbox, lengths, parts );MB_CHK_ERR( rval );
174 
175  /* Get Skin Coordinates, Vertices */
176  double* x = new double[mySkinEnts[0].size()];
177  double* y = new double[mySkinEnts[0].size()];
178  double* z = new double[mySkinEnts[0].size()];
179  rval = myMB->get_coords( mySkinEnts[0], x, y, z );
180  if( rval != MB_SUCCESS )
181  {
182  // Prevent Memory Leak
183  delete[] x;
184  delete[] y;
185  delete[] z;
186  return rval;
187  }
188 
189  // Initialize variable to be used in the loops
190  std::vector< int > toProcs;
191  int xPart, yPart, zPart, xEps, yEps, zEps, baseProc;
192  unsigned long long tup_i = 0, tup_ul = 0, tup_r = 0, count = 0;
193  // These are boolean to determine if the vertex is on close enough to a given border
194  bool xDup, yDup, zDup;
195  bool canWrite = myTup.get_writeEnabled();
196  if( !canWrite ) myTup.enableWriteAccess();
197  // Go through each vertex
198  for( Range::iterator it = mySkinEnts[0].begin(); it != mySkinEnts[0].end(); ++it )
199  {
200  xDup = false;
201  yDup = false;
202  zDup = false;
203  // Figure out which x,y,z partition the element is in.
204  xPart = static_cast< int >( floor( ( x[count] - gbox[0] ) / lengths[0] ) );
205  xPart = ( xPart < parts[0] ? xPart : parts[0] - 1 ); // Make sure it stays within the bounds
206 
207  yPart = static_cast< int >( floor( ( y[count] - gbox[1] ) / lengths[1] ) );
208  yPart = ( yPart < parts[1] ? yPart : parts[1] - 1 ); // Make sure it stays within the bounds
209 
210  zPart = static_cast< int >( floor( ( z[count] - gbox[2] ) / lengths[2] ) );
211  zPart = ( zPart < parts[2] ? zPart : parts[2] - 1 ); // Make sure it stays within the bounds
212 
213  // Figure out the partition with the addition of Epsilon
214  xEps = static_cast< int >( floor( ( x[count] - gbox[0] + myEps ) / lengths[0] ) );
215  yEps = static_cast< int >( floor( ( y[count] - gbox[1] + myEps ) / lengths[1] ) );
216  zEps = static_cast< int >( floor( ( z[count] - gbox[2] + myEps ) / lengths[2] ) );
217 
218  // Figure out if the vertex needs to be sent to multiple procs
219  xDup = ( xPart != xEps && xEps < parts[0] );
220  yDup = ( yPart != yEps && yEps < parts[1] );
221  zDup = ( zPart != zEps && zEps < parts[2] );
222 
223  // Add appropriate processors to the vector
224  baseProc = xPart + yPart * parts[0] + zPart * parts[0] * parts[1];
225  toProcs.push_back( baseProc );
226  if( xDup )
227  {
228  toProcs.push_back( baseProc + 1 ); // Get partition to the right
229  }
230  if( yDup )
231  {
232  // Partition up 1
233  toProcs.push_back( baseProc + parts[0] );
234  }
235  if( zDup )
236  {
237  // Partition above 1
238  toProcs.push_back( baseProc + parts[0] * parts[1] );
239  }
240  if( xDup && yDup )
241  {
242  // Partition up 1 and right 1
243  toProcs.push_back( baseProc + parts[0] + 1 );
244  }
245  if( xDup && zDup )
246  {
247  // Partition right 1 and above 1
248  toProcs.push_back( baseProc + parts[0] * parts[1] + 1 );
249  }
250  if( yDup && zDup )
251  {
252  // Partition up 1 and above 1
253  toProcs.push_back( baseProc + parts[0] * parts[1] + parts[0] );
254  }
255  if( xDup && yDup && zDup )
256  {
257  // Partition right 1, up 1, and above 1
258  toProcs.push_back( baseProc + parts[0] * parts[1] + parts[0] + 1 );
259  }
260  // Grow the tuple list if necessary
261  while( myTup.get_n() + toProcs.size() >= myTup.get_max() )
262  {
263  myTup.resize( myTup.get_max() ? myTup.get_max() + myTup.get_max() / 2 + 1 : 2 );
264  }
265 
266  // Add each proc as a tuple
267  for( std::vector< int >::iterator proc = toProcs.begin(); proc != toProcs.end(); ++proc )
268  {
269  myTup.vi_wr[tup_i++] = *proc;
270  myTup.vul_wr[tup_ul++] = *it;
271  myTup.vr_wr[tup_r++] = x[count];
272  myTup.vr_wr[tup_r++] = y[count];
273  myTup.vr_wr[tup_r++] = z[count];
274  myTup.inc_n();
275  }
276  count++;
277  toProcs.clear();
278  }
279  delete[] x;
280  delete[] y;
281  delete[] z;
282  if( !canWrite ) myTup.disableWriteAccess();
283  return MB_SUCCESS;
284 }
285 
286 // Partition the global box by the number of procs
287 ErrorCode ParallelMergeMesh::PartitionGlobalBox( double* gbox, double* lengths, int* parts )
288 {
289  // Determine the length of each side
290  double xLen = gbox[3] - gbox[0];
291  double yLen = gbox[4] - gbox[1];
292  double zLen = gbox[5] - gbox[2];
293  // avoid division by zero for deciding the way to partition the global box
294  // make all sides of the box at least myEps
295  // it can be zero in some cases, so not possible to compare the ratios :) for best division
296  if( xLen < myEps ) xLen = myEps;
297  if( yLen < myEps ) yLen = myEps;
298  if( zLen < myEps ) zLen = myEps;
299  unsigned numProcs = myPcomm->size();
300 
301  // Partition sides from the longest to shortest lengths
302  // If x is the longest side
303  if( xLen >= yLen && xLen >= zLen )
304  {
305  parts[0] = PartitionSide( xLen, yLen * zLen, numProcs, true );
306  numProcs /= parts[0];
307  // If y is second longest
308  if( yLen >= zLen )
309  {
310  parts[1] = PartitionSide( yLen, zLen, numProcs, false );
311  parts[2] = numProcs / parts[1];
312  }
313  // If z is the longer
314  else
315  {
316  parts[2] = PartitionSide( zLen, yLen, numProcs, false );
317  parts[1] = numProcs / parts[2];
318  }
319  }
320  // If y is the longest side
321  else if( yLen >= zLen )
322  {
323  parts[1] = PartitionSide( yLen, xLen * zLen, numProcs, true );
324  numProcs /= parts[1];
325  // If x is the second longest
326  if( xLen >= zLen )
327  {
328  parts[0] = PartitionSide( xLen, zLen, numProcs, false );
329  parts[2] = numProcs / parts[0];
330  }
331  // If z is the second longest
332  else
333  {
334  parts[2] = PartitionSide( zLen, xLen, numProcs, false );
335  parts[0] = numProcs / parts[2];
336  }
337  }
338  // If z is the longest side
339  else
340  {
341  parts[2] = PartitionSide( zLen, xLen * yLen, numProcs, true );
342  numProcs /= parts[2];
343  // If x is the second longest
344  if( xLen >= yLen )
345  {
346  parts[0] = PartitionSide( xLen, yLen, numProcs, false );
347  parts[1] = numProcs / parts[0];
348  }
349  // If y is the second longest
350  else
351  {
352  parts[1] = PartitionSide( yLen, xLen, numProcs, false );
353  parts[0] = numProcs / parts[1];
354  }
355  }
356 
357  // Divide up each side to give the lengths
358  lengths[0] = xLen / (double)parts[0];
359  lengths[1] = yLen / (double)parts[1];
360  lengths[2] = zLen / (double)parts[2];
361  return MB_SUCCESS;
362 }
363 
364 // Partition a side based on the length ratios
365 int ParallelMergeMesh::PartitionSide( double sideLen, double restLen, unsigned numProcs, bool altRatio )
366 {
367  // If theres only 1 processor, then just return 1
368  if( numProcs == 1 )
369  {
370  return 1;
371  }
372  // Initialize with the ratio of 1 proc
373  double ratio = -DBL_MAX;
374  unsigned factor = 1;
375  // We need to be able to save the last ratio and factor (for comparison)
376  double oldRatio = ratio;
377  double oldFactor = 1;
378 
379  // This is the ratio were shooting for
380  double goalRatio = sideLen / restLen;
381 
382  // Calculate the divisor and numerator power
383  // This avoid if statements in the loop and is useful since both calculations are similar
384  double divisor, p;
385  if( altRatio )
386  {
387  divisor = (double)numProcs * sideLen;
388  p = 3;
389  }
390  else
391  {
392  divisor = (double)numProcs;
393  p = 2;
394  }
395 
396  // Find each possible factor
397  for( unsigned i = 2; i <= numProcs / 2; i++ )
398  {
399  // If it is a factor...
400  if( numProcs % i == 0 )
401  {
402  // We need to save the past factor
403  oldRatio = ratio;
404  oldFactor = factor;
405  // There are 2 different ways to calculate the ratio:
406  // Comparing 1 side to 2 sides: (i*i*i)/(numProcs*x)
407  // Justification: We have a ratio x:y:z (side Lengths) == a:b:c (procs). So a=kx,
408  // b=ky, c=kz. Also, abc=n (numProcs) => bc = n/a. Also, a=kx => k=a/x => 1/k=x/a And so
409  // x/(yz) == (kx)/(kyz) == (kx)/(kykz(1/k)) == a/(bc(x/a)) == a/((n/a)(x/a)) == a^3/(nx).
410  // Comparing 1 side to 1 side: (i*i)/numprocs
411  // Justification: i/(n/i) == i^2/n
412  ratio = pow( (double)i, p ) / divisor;
413  factor = i;
414  // Once we have passed the goal ratio, we can break since we'll only move away from the
415  // goal ratio
416  if( ratio >= goalRatio )
417  {
418  break;
419  }
420  }
421  }
422  // If we haven't reached the goal ratio yet, check out factor = numProcs
423  if( ratio < goalRatio )
424  {
425  oldRatio = ratio;
426  oldFactor = factor;
427  factor = numProcs;
428  ratio = pow( (double)numProcs, p ) / divisor;
429  }
430 
431  // Figure out if our oldRatio is better than ratio
432  if( fabs( ratio - goalRatio ) > fabs( oldRatio - goalRatio ) )
433  {
434  factor = oldFactor;
435  }
436  // Return our findings
437  return factor;
438 }
439 
440 // Id the tuples that are matching
442 {
443  // Counters for accessing tuples more efficiently
444  unsigned long i = 0, mat_i = 0, mat_ul = 0, j = 0, tup_r = 0;
445  double eps2 = myEps * myEps;
446 
447  uint tup_mi, tup_ml, tup_mul, tup_mr;
448  myTup.getTupleSize( tup_mi, tup_ml, tup_mul, tup_mr );
449 
450  bool canWrite = myMatches.get_writeEnabled();
451  if( !canWrite ) myMatches.enableWriteAccess();
452 
453  while( ( i + 1 ) < myTup.get_n() )
454  {
455  // Proximity Comparison
456  double xi = myTup.vr_rd[tup_r], yi = myTup.vr_rd[tup_r + 1], zi = myTup.vr_rd[tup_r + 2];
457 
458  bool done = false;
459  while( !done )
460  {
461  j++;
462  tup_r += tup_mr;
463  if( j >= myTup.get_n() )
464  {
465  break;
466  }
467  CartVect cv( myTup.vr_rd[tup_r] - xi, myTup.vr_rd[tup_r + 1] - yi, myTup.vr_rd[tup_r + 2] - zi );
468  if( cv.length_squared() > eps2 )
469  {
470  done = true;
471  }
472  }
473  // Allocate the tuple list before adding matches
474  while( myMatches.get_n() + ( j - i ) * ( j - i - 1 ) >= myMatches.get_max() )
475  {
477  }
478 
479  // We now know that tuples [i to j) exclusive match.
480  // If n tuples match, n*(n-1) match tuples will be made
481  // tuples are of the form (proc1,proc2,handle1,handle2)
482  if( i + 1 < j )
483  {
484  int kproc = i * tup_mi;
485  unsigned long khand = i * tup_mul;
486  for( unsigned long k = i; k < j; k++ )
487  {
488  int lproc = kproc + tup_mi;
489  unsigned long lhand = khand + tup_mul;
490  for( unsigned long l = k + 1; l < j; l++ )
491  {
492  myMatches.vi_wr[mat_i++] = myTup.vi_rd[kproc]; // proc1
493  myMatches.vi_wr[mat_i++] = myTup.vi_rd[lproc]; // proc2
494  myMatches.vul_wr[mat_ul++] = myTup.vul_rd[khand]; // handle1
495  myMatches.vul_wr[mat_ul++] = myTup.vul_rd[lhand]; // handle2
496  myMatches.inc_n();
497 
498  myMatches.vi_wr[mat_i++] = myTup.vi_rd[lproc]; // proc1
499  myMatches.vi_wr[mat_i++] = myTup.vi_rd[kproc]; // proc2
500  myMatches.vul_wr[mat_ul++] = myTup.vul_rd[lhand]; // handle1
501  myMatches.vul_wr[mat_ul++] = myTup.vul_rd[khand]; // handle2
502  myMatches.inc_n();
503  lproc += tup_mi;
504  lhand += tup_mul;
505  }
506  kproc += tup_mi;
507  khand += tup_mul;
508  } // End for(int k...
509  }
510  i = j;
511  } // End while(i+1<tup.n)
512 
513  if( !canWrite ) myMatches.disableWriteAccess();
514  return MB_SUCCESS;
515 }
516 
517 // Sort the matching tuples so that vertices can be tagged accurately
519 {
520  TupleList::buffer buf( mySkinEnts[0].size() );
521  // Sorts are necessary to check for doubles
522  // Sort by remote handle
523  myMatches.sort( 3, &buf );
524  // Sort by matching proc
525  myMatches.sort( 1, &buf );
526  // Sort by local handle
527  myMatches.sort( 2, &buf );
528  buf.reset();
529  return MB_SUCCESS;
530 }
531 
532 // Tag the shared elements using existing PComm functionality
534 {
535  // Manipulate the matches list to tag vertices and entities
536  // Set up proc ents
537  Range proc_ents;
538  ErrorCode rval;
539 
540  // get the entities in the partition sets
541  for( Range::iterator rit = myPcomm->partitionSets.begin(); rit != myPcomm->partitionSets.end(); ++rit )
542  {
543  Range tmp_ents;
544  rval = myMB->get_entities_by_handle( *rit, tmp_ents, true );
545  if( MB_SUCCESS != rval )
546  {
547  return rval;
548  }
549  proc_ents.merge( tmp_ents );
550  }
551  if( myMB->dimension_from_handle( *proc_ents.rbegin() ) != myMB->dimension_from_handle( *proc_ents.begin() ) )
552  {
553  Range::iterator lower = proc_ents.lower_bound( CN::TypeDimensionMap[0].first ),
554  upper = proc_ents.upper_bound( CN::TypeDimensionMap[dim - 1].second );
555  proc_ents.erase( lower, upper );
556  }
557 
558  // This vector doesn't appear to be used but its in resolve_shared_ents
559  int maxp = -1;
560  std::vector< int > sharing_procs( MAX_SHARING_PROCS );
561  std::fill( sharing_procs.begin(), sharing_procs.end(), maxp );
562 
563  // get ents shared by 1 or n procs
564  std::map< std::vector< int >, std::vector< EntityHandle > > proc_nranges;
565  Range proc_verts;
566  rval = myMB->get_adjacencies( proc_ents, 0, false, proc_verts, Interface::UNION );
567  if( rval != MB_SUCCESS )
568  {
569  return rval;
570  }
571 
572  rval = myPcomm->tag_shared_verts( myMatches, proc_nranges, proc_verts );
573  if( rval != MB_SUCCESS )
574  {
575  return rval;
576  }
577 
578  // get entities shared by 1 or n procs
579  rval = myPcomm->get_proc_nvecs( dim, dim - 1, &mySkinEnts[0], proc_nranges );
580  if( rval != MB_SUCCESS )
581  {
582  return rval;
583  }
584 
585  // create the sets for each interface; store them as tags on
586  // the interface instance
587  Range iface_sets;
588  rval = myPcomm->create_interface_sets( proc_nranges );
589  if( rval != MB_SUCCESS )
590  {
591  return rval;
592  }
593  // establish comm procs and buffers for them
594  std::set< unsigned int > procs;
595  rval = myPcomm->get_interface_procs( procs, true );
596  if( rval != MB_SUCCESS )
597  {
598  return rval;
599  }
600 
601  // resolve shared entity remote handles; implemented in ghost cell exchange
602  // code because it's so similar
603  rval = myPcomm->exchange_ghost_cells( -1, -1, 0, true, true );
604  if( rval != MB_SUCCESS )
605  {
606  return rval;
607  }
608  // now build parent/child links for interface sets
609  rval = myPcomm->create_iface_pc_links();
610  return rval;
611 }
612 
613 // Make sure to free up any allocated data
614 // Need to avoid a double free
616 {
617  // The reset operation is now safe and avoids a double free()
618  myMatches.reset();
619  myTup.reset();
620  myCD.reset();
621 }
622 
623 // Simple quick sort to real
625 {
626  bool canWrite = tup.get_writeEnabled();
627  if( !canWrite ) tup.enableWriteAccess();
628 
629  uint mi, ml, mul, mr;
630  tup.getTupleSize( mi, ml, mul, mr );
631  PerformRealSort( tup, 0, tup.get_n(), eps, mr );
632 
633  if( !canWrite ) tup.disableWriteAccess();
634 }
635 
636 // Swap around tuples
637 void ParallelMergeMesh::SwapTuples( TupleList& tup, unsigned long a, unsigned long b )
638 {
639  if( a == b ) return;
640 
641  uint mi, ml, mul, mr;
642  tup.getTupleSize( mi, ml, mul, mr );
643 
644  // Swap mi
645  unsigned long a_val = a * mi, b_val = b * mi;
646  for( unsigned long i = 0; i < mi; i++ )
647  {
648  sint t = tup.vi_rd[a_val];
649  tup.vi_wr[a_val] = tup.vi_rd[b_val];
650  tup.vi_wr[b_val] = t;
651  a_val++;
652  b_val++;
653  }
654  // Swap ml
655  a_val = a * ml;
656  b_val = b * ml;
657  for( unsigned long i = 0; i < ml; i++ )
658  {
659  slong t = tup.vl_rd[a_val];
660  tup.vl_wr[a_val] = tup.vl_rd[b_val];
661  tup.vl_wr[b_val] = t;
662  a_val++;
663  b_val++;
664  }
665  // Swap mul
666  a_val = a * mul;
667  b_val = b * mul;
668  for( unsigned long i = 0; i < mul; i++ )
669  {
670  Ulong t = tup.vul_rd[a_val];
671  tup.vul_wr[a_val] = tup.vul_rd[b_val];
672  tup.vul_wr[b_val] = t;
673  a_val++;
674  b_val++;
675  }
676  // Swap mr
677  a_val = a * mr;
678  b_val = b * mr;
679  for( unsigned long i = 0; i < mr; i++ )
680  {
681  realType t = tup.vr_rd[a_val];
682  tup.vr_wr[a_val] = tup.vr_rd[b_val];
683  tup.vr_wr[b_val] = t;
684  a_val++;
685  b_val++;
686  }
687 }
688 
689 // Perform the sorting of a tuple by real
690 // To sort an entire tuple_list, call (tup,0,tup.n,epsilon)
692  unsigned long left,
693  unsigned long right,
694  double eps,
695  uint tup_mr )
696 {
697  // If list size is only 1 or 0 return
698  if( left + 1 >= right )
699  {
700  return;
701  }
702  unsigned long swap = left, tup_l = left * tup_mr, tup_t = tup_l + tup_mr;
703 
704  // Swap the median with the left position for a (hopefully) better split
705  SwapTuples( tup, left, ( left + right ) / 2 );
706 
707  // Partition the data
708  for( unsigned long t = left + 1; t < right; t++ )
709  {
710  // If the left value(pivot) is greater than t_val, swap it into swap
711  if( TupleGreaterThan( tup, tup_l, tup_t, eps, tup_mr ) )
712  {
713  swap++;
714  SwapTuples( tup, swap, t );
715  }
716  tup_t += tup_mr;
717  }
718 
719  // Swap so that position swap is in the correct position
720  SwapTuples( tup, left, swap );
721 
722  // Sort left and right of swap
723  PerformRealSort( tup, left, swap, eps, tup_mr );
724  PerformRealSort( tup, swap + 1, right, eps, tup_mr );
725 }
726 
727 // Note, this takes the actual tup.vr[] index (aka i*tup.mr)
729  unsigned long vrI,
730  unsigned long vrJ,
731  double eps,
732  uint tup_mr )
733 {
734  unsigned check = 0;
735  while( check < tup_mr )
736  {
737  // If the values are the same
738  if( fabs( tup.vr_rd[vrI + check] - tup.vr_rd[vrJ + check] ) <= eps )
739  {
740  check++;
741  continue;
742  }
743  // If I greater than J
744  else if( tup.vr_rd[vrI + check] > tup.vr_rd[vrJ + check] )
745  {
746  return true;
747  }
748  // If J greater than I
749  else
750  {
751  return false;
752  }
753  }
754  // All Values are the same
755  return false;
756 }
757 
758 } // End namespace moab