NAMD
SortAtoms.C
Go to the documentation of this file.
1 
7 #include "SortAtoms.h"
8 #include "NamdTypes.h"
9 #include <algorithm>
10 
11 // #include "charm++.h"
12 
13 
14 struct sortop_base {
15  const FullAtom * const a;
16  sortop_base(const FullAtom* atoms) : a(atoms) { }
17 };
18 
19 struct sortop_x : public sortop_base {
20  sortop_x(const FullAtom* atoms) : sortop_base(atoms) { }
21  bool operator() (int i, int j) const {
22  return ( a[i].position.x < a[j].position.x );
23  }
24 };
25 
26 struct sortop_y : public sortop_base {
27  sortop_y(const FullAtom* atoms) : sortop_base(atoms) { }
28  bool operator() (int i, int j) const {
29  return ( a[i].position.y < a[j].position.y );
30  }
31 };
32 
33 struct sortop_z : public sortop_base {
34  sortop_z(const FullAtom* atoms) : sortop_base(atoms) { }
35  bool operator() (int i, int j) const {
36  return ( a[i].position.z < a[j].position.z );
37  }
38 };
39 
40 
41 static void partition(int *order, const FullAtom *atoms, int begin, int end) {
42 
43  // Applies orthogonal recursive bisection with splittings limited
44  // to multiples of 32 for warps and a final split on multiples of 16.
45 
46 #ifdef NAMD_AVXTILES
47 #define PARTWIDTH 16
48 #else
49 #define PARTWIDTH 32
50 #endif
51 
52  int split;
53  // must be a multiple of 32 or 16 between begin and end to split at
54  if ( begin/PARTWIDTH < (end-1)/PARTWIDTH ) {
55  // find a multiple of 32 near the median
56  split = ((begin + end + PARTWIDTH) / (PARTWIDTH*2)) * PARTWIDTH;
57  } else if ( begin/(PARTWIDTH/2) < (end-1)/(PARTWIDTH/2) ) {
58  // find a multiple of 16 near the median
59  split = ((begin + end + (PARTWIDTH/2)) / PARTWIDTH) * (PARTWIDTH/2);
60  } else {
61  return;
62  }
63 
64  BigReal xmin, ymin, zmin, xmax, ymax, zmax;
65  {
66  const Position &pos = atoms[order[begin]].position;
67  xmin = pos.x;
68  ymin = pos.y;
69  zmin = pos.z;
70  xmax = pos.x;
71  ymax = pos.y;
72  zmax = pos.z;
73  }
74  for ( int i=begin+1; i<end; ++i ) {
75  const Position &pos = atoms[order[i]].position;
76  if ( pos.x < xmin ) { xmin = pos.x; }
77  if ( pos.y < ymin ) { ymin = pos.y; }
78  if ( pos.z < zmin ) { zmin = pos.z; }
79  if ( pos.x > xmax ) { xmax = pos.x; }
80  if ( pos.y > ymax ) { ymax = pos.y; }
81  if ( pos.z > zmax ) { zmax = pos.z; }
82  }
83  xmax -= xmin;
84  ymax -= ymin;
85  zmax -= zmin;
86 
87 #define NTH_ELEMENT(BEGIN,SPLIT,END,OP) std::nth_element(BEGIN,SPLIT,END,OP)
88 #if defined(NAMD_CUDA) && defined(__GNUC_PATCHLEVEL__)
89 #if __GNUC__ == 4 && __GNUC_MINOR__ == 8 && __GNUC_PATCHLEVEL__ == 2
90 #define NTH_ELEMENT(BEGIN,SPLIT,END,OP) std::sort(BEGIN,END,OP)
91 #warning gcc 4.8.2 std::nth_element would segfault (see http://gcc.gnu.org/bugzilla/show_bug.cgi?id=58800)
92 #endif
93 #endif
94 
95  if ( xmax >= ymax && xmax >= zmax ) {
96  NTH_ELEMENT(order+begin, order+split, order+end, sortop_x(atoms));
97  } else if ( ymax >= xmax && ymax >= zmax ) {
98  NTH_ELEMENT(order+begin, order+split, order+end, sortop_y(atoms));
99  } else {
100  NTH_ELEMENT(order+begin, order+split, order+end, sortop_z(atoms));
101  }
102 
103  if ( split & PARTWIDTH/2 ) return;
104 
105 #undef PARTWIDTH
106 
107  // recursively partition before and after split
108  partition(order, atoms, begin, split);
109  partition(order, atoms, split, end);
110 
111 }
112 
113 #if defined(NAMD_CUDA) && defined(__GNUC_PATCHLEVEL__)
114 #if __GNUC__ == 4 && __GNUC_MINOR__ == 8 && __GNUC_PATCHLEVEL__ == 2
115 // #error gcc 4.8.2 std::nth_element would segfault (see http://gcc.gnu.org/bugzilla/show_bug.cgi?id=58800)
116 #endif
117 #endif
118 
119 void sortAtomsForCUDA(int *order, const FullAtom *atoms, int nfree, int n) {
120 
121  // partition free atoms
122  // CkPrintf("%d %d\n", 0, nfree);
123  partition(order, atoms, 0, nfree);
124 
125  // partition fixed atoms
126  // CkPrintf("%d %d\n", nfree, n);
127  partition(order, atoms, nfree, n);
128 
129 }
130 
131 void sortAtomsForPatches(int *order, int *breaks,
132  const FullAtom *atoms, int nmgrps, int natoms,
133  int ni, int nj, int nk) {
134 
135 //CkPrintf("sorting %d atoms in %d groups to %d x %d x %d\n",
136 // natoms, nmgrps, nk, nj, ni);
137  std::sort(order, order+nmgrps, sortop_z(atoms));
138  int pid = 0;
139  int ibegin = 0;
140  int nai = 0;
141  for ( int ip=0; ip < ni; ++ip ) {
142  int naj = nai;
143  int targi = nai + (natoms - nai - 1) / (ni - ip) + 1;
144  int iend;
145  for ( iend=ibegin; iend<nmgrps; ++iend ) {
146  int mgs = atoms[order[iend]].migrationGroupSize;
147  if (nai + mgs <= targi) nai += mgs;
148  else break;
149  }
150 //CkPrintf(" Z %d %d (%d) %d\n", ibegin, iend, iend-ibegin, nai);
151  std::sort(order+ibegin, order+iend, sortop_y(atoms));
152  int jbegin = ibegin;
153  for ( int jp=0; jp < nj; ++jp ) {
154  int nak = naj;
155  int targj = naj + (nai - naj - 1) / (nj - jp) + 1;
156  int jend;
157  for ( jend=jbegin; jend<iend; ++jend ) {
158  int mgs = atoms[order[jend]].migrationGroupSize;
159  if (naj + mgs <= targj) naj += mgs;
160  else break;
161  }
162 
163 //CkPrintf(" Y %d %d (%d) %d\n", jbegin, jend, jend-jbegin, naj);
164  std::sort(order+jbegin, order+jend, sortop_x(atoms));
165  int kbegin = jbegin;
166  for ( int kp=0; kp < nk; ++kp ) {
167  int targk = nak + (naj - nak - 1) / (nk - kp) + 1;
168  int kend;
169  for ( kend=kbegin; kend<jend; ++kend ) {
170  int mgs = atoms[order[kend]].migrationGroupSize;
171  if (nak + mgs <= targk) nak += mgs;
172  else break;
173 //CkPrintf(" atom %d %d %.2f\n", atoms[order[kend]].id, mgs,
174 // atoms[order[kend]].position.x);
175  }
176 //CkPrintf(" X %d %d (%d) %d\n", kbegin, kend, kend-kbegin, nak);
177  breaks[pid++] = kend;
178  kbegin = kend;
179  }
180  jbegin = jend;
181  }
182  ibegin = iend;
183  }
184 
185 }
186 
sortop_x(const FullAtom *atoms)
Definition: SortAtoms.C:20
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:41
Definition: Vector.h:64
static __thread atom * atoms
BigReal z
Definition: Vector.h:66
const FullAtom *const a
Definition: SortAtoms.C:15
bool operator()(int i, int j) const
Definition: SortAtoms.C:35
Position position
Definition: NamdTypes.h:53
sortop_y(const FullAtom *atoms)
Definition: SortAtoms.C:27
void sortAtomsForCUDA(int *order, const FullAtom *atoms, int nfree, int n)
Definition: SortAtoms.C:119
sortop_base(const FullAtom *atoms)
Definition: SortAtoms.C:16
#define NTH_ELEMENT(BEGIN, SPLIT, END, OP)
#define order
Definition: PmeRealSpace.C:235
#define PARTWIDTH
BigReal x
Definition: Vector.h:66
std::vector< std::string > split(const std::string &text, std::string delimiter)
Definition: MoleculeQM.C:73
int migrationGroupSize
Definition: NamdTypes.h:117
BlockRadixSort::TempStorage sort
BigReal y
Definition: Vector.h:66
bool operator()(int i, int j) const
Definition: SortAtoms.C:21
bool operator()(int i, int j) const
Definition: SortAtoms.C:28
void sortAtomsForPatches(int *order, int *breaks, const FullAtom *atoms, int nmgrps, int natoms, int ni, int nj, int nk)
Definition: SortAtoms.C:131
double BigReal
Definition: common.h:112
sortop_z(const FullAtom *atoms)
Definition: SortAtoms.C:34