SortAtoms.C

Go to the documentation of this file.
00001 
00007 #include "SortAtoms.h"
00008 #include "NamdTypes.h"
00009 #include <algorithm>
00010 
00011 // #include "charm++.h"
00012 
00013 
00014 struct sortop_base {
00015   const FullAtom * const a;
00016   sortop_base(const FullAtom* atoms) : a(atoms) { }
00017 };
00018 
00019 struct sortop_x : public sortop_base {
00020   sortop_x(const FullAtom* atoms) : sortop_base(atoms) { }
00021   bool operator() (int i, int j) const {
00022     return ( a[i].position.x < a[j].position.x );
00023   }
00024 };
00025 
00026 struct sortop_y : public sortop_base {
00027   sortop_y(const FullAtom* atoms) : sortop_base(atoms) { }
00028   bool operator() (int i, int j) const {
00029     return ( a[i].position.y < a[j].position.y );
00030   }
00031 };
00032 
00033 struct sortop_z : public sortop_base {
00034   sortop_z(const FullAtom* atoms) : sortop_base(atoms) { }
00035   bool operator() (int i, int j) const {
00036     return ( a[i].position.z < a[j].position.z );
00037   }
00038 };
00039 
00040 
00041 static void partition(int *order, const FullAtom *atoms, int begin, int end) {
00042 
00043   //  Applies orthogonal recursive bisection with splittings limited
00044   //  to multiples of 32 for warps and a final split on multiples of 16.
00045 
00046   int split;
00047   // must be a multiple of 32 or 16 between begin and end to split at
00048   if ( begin/32 < (end-1)/32 ) {
00049     // find a multiple of 32 near the median
00050     split = ((begin + end + 32) / 64) * 32;
00051   } else if ( begin/16 < (end-1)/16 ) {
00052     // find a multiple of 16 near the median
00053     split = ((begin + end + 16) / 32) * 16;
00054   } else {
00055     return;
00056   }
00057 
00058   BigReal xmin, ymin, zmin, xmax, ymax, zmax;
00059   {
00060     const Position &pos = atoms[order[begin]].position;
00061     xmin = pos.x;
00062     ymin = pos.y;
00063     zmin = pos.z;
00064     xmax = pos.x;
00065     ymax = pos.y;
00066     zmax = pos.z;
00067   }
00068   for ( int i=begin+1; i<end; ++i ) {
00069     const Position &pos = atoms[order[i]].position;
00070     if ( pos.x < xmin ) { xmin = pos.x; }
00071     if ( pos.y < ymin ) { ymin = pos.y; }
00072     if ( pos.z < zmin ) { zmin = pos.z; }
00073     if ( pos.x > xmax ) { xmax = pos.x; }
00074     if ( pos.y > ymax ) { ymax = pos.y; }
00075     if ( pos.z > zmax ) { zmax = pos.z; }
00076   }
00077   xmax -= xmin;
00078   ymax -= ymin;
00079   zmax -= zmin;
00080 
00081 #define NTH_ELEMENT(BEGIN,SPLIT,END,OP) std::nth_element(BEGIN,SPLIT,END,OP)
00082 #if defined(NAMD_CUDA) && defined(__GNUC_PATCHLEVEL__)
00083 #if __GNUC__ == 4 && __GNUC_MINOR__ == 8 && __GNUC_PATCHLEVEL__ == 2
00084 #define NTH_ELEMENT(BEGIN,SPLIT,END,OP) std::sort(BEGIN,END,OP)
00085 #warning gcc 4.8.2 std::nth_element would segfault (see http://gcc.gnu.org/bugzilla/show_bug.cgi?id=58800)
00086 #endif
00087 #endif
00088 
00089   if ( xmax >= ymax && xmax >= zmax ) {
00090     NTH_ELEMENT(order+begin, order+split, order+end, sortop_x(atoms));
00091   } else if ( ymax >= xmax && ymax >= zmax ) {
00092     NTH_ELEMENT(order+begin, order+split, order+end, sortop_y(atoms));
00093   } else {
00094     NTH_ELEMENT(order+begin, order+split, order+end, sortop_z(atoms));
00095   }
00096 
00097   if ( split & 16 ) return;
00098 
00099   // recursively partition before and after split
00100   partition(order, atoms, begin, split);
00101   partition(order, atoms, split, end);
00102 
00103 }
00104 
00105 #if defined(NAMD_CUDA) && defined(__GNUC_PATCHLEVEL__)
00106 #if __GNUC__ == 4 && __GNUC_MINOR__ == 8 && __GNUC_PATCHLEVEL__ == 2
00107 // #error gcc 4.8.2 std::nth_element would segfault (see http://gcc.gnu.org/bugzilla/show_bug.cgi?id=58800)
00108 #endif
00109 #endif
00110 
00111 void sortAtomsForCUDA(int *order, const FullAtom *atoms, int nfree, int n) {
00112 
00113   // partition free atoms
00114   // CkPrintf("%d %d\n", 0, nfree);
00115   partition(order, atoms, 0, nfree);
00116 
00117   // partition fixed atoms
00118   // CkPrintf("%d %d\n", nfree, n);
00119   partition(order, atoms, nfree, n);
00120 
00121 }
00122 
00123 void sortAtomsForPatches(int *order, int *breaks,
00124                          const FullAtom *atoms, int nmgrps, int natoms,
00125                          int ni, int nj, int nk) {
00126 
00127 //CkPrintf("sorting %d atoms in %d groups to %d x %d x %d\n",
00128 //    natoms, nmgrps, nk, nj, ni);
00129   std::sort(order, order+nmgrps, sortop_z(atoms));
00130   int pid = 0;
00131   int ibegin = 0;
00132   int nai = 0;
00133   for ( int ip=0; ip < ni; ++ip ) {
00134     int naj = nai;
00135     int targi = nai + (natoms - nai - 1) / (ni - ip) + 1;
00136     int iend;
00137     for ( iend=ibegin; iend<nmgrps; ++iend ) { 
00138       int mgs = atoms[order[iend]].migrationGroupSize;
00139       if (nai + mgs <= targi) nai += mgs;
00140       else break;
00141     }
00142 //CkPrintf("  Z %d %d (%d) %d\n", ibegin, iend, iend-ibegin, nai);
00143     std::sort(order+ibegin, order+iend, sortop_y(atoms));
00144     int jbegin = ibegin;
00145     for ( int jp=0; jp < nj; ++jp ) {
00146       int nak = naj;
00147       int targj = naj + (nai - naj - 1) / (nj - jp) + 1;
00148       int jend;
00149       for ( jend=jbegin; jend<iend; ++jend ) { 
00150         int mgs = atoms[order[jend]].migrationGroupSize;
00151         if (naj + mgs <= targj) naj += mgs;
00152         else break;
00153       }
00154 
00155 //CkPrintf("    Y %d %d (%d) %d\n", jbegin, jend, jend-jbegin, naj);
00156       std::sort(order+jbegin, order+jend, sortop_x(atoms));
00157       int kbegin = jbegin;
00158       for ( int kp=0; kp < nk; ++kp ) {
00159         int targk = nak + (naj - nak - 1) / (nk - kp) + 1;
00160         int kend;  
00161         for ( kend=kbegin; kend<jend; ++kend ) {
00162           int mgs = atoms[order[kend]].migrationGroupSize;
00163           if (nak + mgs <= targk) nak += mgs;
00164           else break;
00165 //CkPrintf("        atom %d %d %.2f\n", atoms[order[kend]].id, mgs,
00166 //                  atoms[order[kend]].position.x);
00167         }
00168 //CkPrintf("      X %d %d (%d) %d\n", kbegin, kend, kend-kbegin, nak);
00169         breaks[pid++] = kend;
00170         kbegin = kend;
00171       }
00172       jbegin = jend;
00173     }
00174     ibegin = iend;
00175   }
00176 
00177 }
00178 

Generated on Sun Nov 19 01:17:15 2017 for NAMD by  doxygen 1.4.7