MsmMap.h

Go to the documentation of this file.
00001 
00007 #ifndef MSMMAP_H
00008 #define MSMMAP_H
00009 
00010 // SSE and AVX vector intrinsics and memory alignment macros
00011 #if defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE)
00012 #include <emmintrin.h>  // SSE2
00013 #if defined(__INTEL_COMPILER)
00014 #define __align(X) __declspec(align(X) )
00015 #elif defined(__PGI)
00016 #define __align(X)  __attribute__((aligned(X) ))
00017 #define MISSING_mm_cvtsd_f64
00018 #elif defined(__GNUC__)
00019 #define __align(X)  __attribute__((aligned(X) ))
00020 #if (__GNUC__ < 4)
00021 #define MISSING_mm_cvtsd_f64
00022 #endif
00023 #else
00024 #define __align(X) __declspec(align(X) )
00025 #endif
00026 #endif
00027 
00028 // migration of MSM computes not currently enabled
00029 #define MSM_MIGRATION
00030 #undef MSM_MIGRATION
00031 
00032 // for debugging MSM migration
00033 #define DEBUG_MSM_MIGRATE
00034 #undef DEBUG_MSM_MIGRATE
00035 
00036 #define MSM_MAX_BLOCK_SIZE 8
00037 #define MSM_MAX_BLOCK_VOLUME \
00038   (MSM_MAX_BLOCK_SIZE * MSM_MAX_BLOCK_SIZE * MSM_MAX_BLOCK_SIZE)
00039 
00040 #define MSM_C1VECTOR_MAX_BLOCK_SIZE (MSM_MAX_BLOCK_SIZE / 2)
00041 #define MSM_C1VECTOR_MAX_BLOCK_VOLUME \
00042   (MSM_C1VECTOR_MAX_BLOCK_SIZE * \
00043    MSM_C1VECTOR_MAX_BLOCK_SIZE * \
00044    MSM_C1VECTOR_MAX_BLOCK_SIZE)
00045 
00046 #define DEBUG_MSM
00047 #undef DEBUG_MSM
00048 
00049 #define DEBUG_MSM_VERBOSE
00050 #undef DEBUG_MSM_VERBOSE
00051 
00052 #define DEBUG_MSM_GRID
00053 #undef DEBUG_MSM_GRID
00054 
00055 // assert macro
00056 #undef ASSERT
00057 #ifdef DEBUG_MSM
00058 #define ASSERT(expr) \
00059   do { \
00060     if ( !(expr) ) { \
00061       char msg[100]; \
00062       snprintf(msg, sizeof(msg), "ASSERT: \"%s\" " \
00063           "(%s, %d)\n", #expr, __FILE__, __LINE__); \
00064       NAMD_die(msg); \
00065     } \
00066   } while (0)
00067 #else
00068 #define ASSERT(expr)
00069 #endif 
00070 
00071 
00072 // employ mixed precision
00073 // (but allow easy change to all double precision for comparison)
00074 typedef float Float;
00075 typedef double Double;
00076 
00077 
00079   //
00080   // Vector and matrix elements for C1 Hermite interpolation.
00081   //
00083 
00084   enum { C1_VECTOR_SIZE = 8, C1_MATRIX_SIZE = 8*8 };
00085 
00086   struct C1Vector {
00087     Float velem[C1_VECTOR_SIZE];
00088     C1Vector(Float r=0) { set(r); }
00089     void set(Float r) {
00090       for (int n=0;  n < C1_VECTOR_SIZE;  n++)  velem[n] = r;
00091     }
00092     C1Vector& operator+=(const C1Vector& v) {
00093       for (int n=0;  n < C1_VECTOR_SIZE;  n++)  velem[n] += v.velem[n];
00094       return(*this);
00095     }
00096     friend Float operator*(const C1Vector& u, const C1Vector& v) {
00097       Float r=0;
00098       for (int n=0;  n < C1_VECTOR_SIZE;  n++)  r += u.velem[n] * v.velem[n];
00099       return r;
00100     }
00101     friend C1Vector operator+(const C1Vector& u, const C1Vector& v) {
00102       C1Vector w;
00103       for (int n=0;  n < C1_VECTOR_SIZE;  n++) {
00104         w.velem[n] = u.velem[n] + v.velem[n];
00105       }
00106       return w;
00107     }
00108   };
00109 
00110   struct C1Matrix {
00111     Float melem[C1_MATRIX_SIZE];
00112     C1Matrix(Float r=0) { set(r); }
00113     void set(Float r) {
00114       for (int n=0;  n < C1_MATRIX_SIZE;  n++)  melem[n] = 0;
00115     }
00116     friend C1Vector operator*(const C1Matrix& m, const C1Vector& u) {
00117       C1Vector v;
00118 
00119       // XXX not tested yet
00120 #if 1 && (defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE))
00121       // Hand-coded SSE2 vectorization
00122       // This loop requires that the single-precision input arrays be 
00123       // aligned on 16-byte boundaries, such that array[index % 4 == 0] 
00124       // can be safely accessed with aligned load/store operations
00125       for (int k=0, j=0;  j < C1_VECTOR_SIZE;  j++) {
00126         __m128 melem4 = _mm_load_ps(&m.melem[k]);
00127         __m128 uelem4 = _mm_load_ps(&u.velem[0]);
00128         __m128 tmp4 = _mm_mul_ps(melem4, uelem4); 
00129         melem4 = _mm_load_ps(&m.melem[k+4]);
00130         uelem4 = _mm_load_ps(&u.velem[4]);
00131         tmp4 = _mm_add_ps(tmp4, _mm_mul_ps(melem4, uelem4)); 
00132 
00133         // do a 4-element reduction and accumulate result
00134         __m128 sum4 = tmp4;
00135         sum4 = _mm_shuffle_ps(sum4, sum4, _MM_SHUFFLE(2, 3, 0, 1));
00136         sum4 = _mm_add_ps(sum4, tmp4);
00137         tmp4 = sum4;
00138         sum4 = _mm_shuffle_ps(sum4, sum4, _MM_SHUFFLE(1, 0, 3, 2));
00139         sum4 = _mm_add_ps(sum4, tmp4);
00140 
00141         // all 4 elements are now set to the sum
00142         float sum;
00143         _mm_store_ss(&sum, sum4); // store lowest element
00144         v.velem[j] += sum;
00145         k+=8;
00146       }
00147 #elif 0 && (defined(__AVX__) && ! defined(NAMD_DISABLE_SSE))
00148       // Hand-coded AVX vectorization
00149       // This loop requires that the single-precision input arrays be 
00150       // aligned on 32-byte boundaries, such that array[index % 8 == 0] 
00151       // can be safely accessed with aligned load/store operations
00152       for (int k=0, j=0;  j < C1_VECTOR_SIZE;  j++) {
00153         __m256 melem8 = _mm256_load_ps(&m.melem[k]);
00154         __m256 uelem8 = _mm256_load_ps(&u.velem[0]);
00155         __m256 tmp8 = _mm256_mul_ps(melem8, uelem8); 
00156 
00157         // XXX this still needs to be rewritten a bit for AVX
00158         // do an 8-element reduction and accumulate result
00159         __m256 sum8 = tmp8;
00160         sum8 = _mm256_hadd_ps(sum8, sum8);
00161         sum8 = _mm256_hadd_ps(sum8, sum8);
00162         tmp8 = sum8;
00163         tmp8 = _mm256_permute2f128_ps(tmp8, tmp8, 1);
00164         sum8 = _mm256_hadd_ps(tmp8, sum8);
00165 
00166         // all 8 elements are now set to the sum
00167         float sum;
00168         _mm_store_ss(&sum, sum8); // store lowest element
00169         v.velem[j] += sum;
00170         k+=8;
00171       }
00172 #else
00173 #if defined(__INTEL_COMPILER)
00174 #pragma vector always
00175 #endif
00176       for (int k=0, j=0;  j < C1_VECTOR_SIZE;  j++) {
00177         for (int i = 0;  i < C1_VECTOR_SIZE;  i++, k++) {
00178           v.velem[j] += m.melem[k] * u.velem[i];
00179         }
00180       }
00181 #endif
00182       return v;
00183     }
00184   };
00185 
00186   // index vector based on mixed partial derivatives in x,y,z
00187   enum { D000=0, D100, D010, D001, D110, D101, D011, D111 };
00188 
00189   // index matrix using 2 vector indexes, row-major column ordering,
00190   // defining partial derivatives of g(xj,yj,zj,xi,yi,zi)
00191 #define C1INDEX(drj,dri)  ((drj)*C1_VECTOR_SIZE + (dri))
00192 
00193 
00194 namespace msm {
00195 
00197   //
00198   // Resizable Array class
00199   //
00201 
00202   template <class T> class Array;
00203 
00204   template <class T>
00205   void swap(Array<T>& s, Array<T>& t);
00206 
00207   template <class T>
00208   class Array {
00209     public:
00210       Array() : abuffer(0), alen(0), amax(0) { }
00211       Array(int n) : abuffer(0), alen(0), amax(0) { resize(n); }
00212       Array(const Array& a) : abuffer(0), alen(0), amax(0) { copy(a); }
00213       ~Array() { setmax(0); }
00214       Array& operator=(const Array& a) {
00215         if (this != &a) copy(a);  // don't allow self-assignment
00216         return(*this);
00217       }
00218       int len() const { return alen; }
00219       int max() const { return amax; }
00220       const T& operator[](int i) const {
00221 #ifdef DEBUG_MSM
00222         return elem(i);
00223 #else
00224         return abuffer[i];
00225 #endif
00226       }
00227       const T& elem(int i) const {
00228         if (i < 0 || i >= alen) {
00229           char msg[100];
00230           snprintf(msg, sizeof(msg), "Array index:  alen=%d, i=%d\n", alen, i);
00231           NAMD_die(msg);
00232         }
00233         return abuffer[i];
00234       }
00235       T& operator[](int i) {
00236 #ifdef DEBUG_MSM
00237         return elem(i);
00238 #else
00239         return abuffer[i];
00240 #endif
00241       }
00242       T& elem(int i) {
00243         if (i < 0 || i >= alen) {
00244           char msg[100];
00245           snprintf(msg, sizeof(msg), "Array index:  alen=%d, i=%d\n", alen, i);
00246           NAMD_die(msg);
00247         }
00248         return abuffer[i];
00249       }
00250       void append(const T& t) {
00251         if (alen==amax) setmax(2*amax+1);
00252         abuffer[alen++] = t;
00253       }
00254       void resize(int n) {
00255         if (n > amax) setmax(n);
00256         alen = n;
00257       }
00258       void setmax(int m);
00259       const T *buffer() const { return abuffer; }
00260       T *buffer() { return abuffer; }
00261       const T *buffer(int& n) const { n = alen; return abuffer; }
00262       T *buffer(int& n) { n = alen; return abuffer; }
00263       void reset(const T& t) {
00264         for (int n = 0;  n < alen;  n++)  abuffer[n] = t;
00265       }
00266       friend void swap<T>(Array&, Array&);
00267 #ifdef DEBUG_MSM
00268       void print(const char *s=0) const {
00269         if (s) printf("PRINTING DATA FOR ARRAY \"%s\":\n", s);
00270         printf("abuffer=%p\n  alen=%d  amax=%d\n",
00271             (void *) abuffer, alen, amax);
00272       }
00273 #endif
00274     protected:
00275       T *abuffer;
00276       int alen, amax;
00277       void copy(const Array& a);
00278   };
00279 
00280   template <class T>
00281   void Array<T>::setmax(int m) {
00282     if (m == amax) return;
00283     else if (m > 0) {
00284       T *newbuffer = new T[m];
00285       if ( ! newbuffer) {
00286         char msg[100];
00287         snprintf(msg, sizeof(msg),
00288             "Can't allocate %lu KB for msm::Array\n",
00289             (unsigned long)(m * sizeof(T) / 1024));
00290         NAMD_die(msg);
00291       }
00292       if (alen > m) alen = m;  // new buffer is shorter than old buffer
00293       for (int i = 0;  i < alen;  i++) {
00294         newbuffer[i] = abuffer[i];
00295       }
00296       delete[] abuffer;
00297       abuffer = newbuffer;
00298       amax = m;
00299     }
00300     else {  // consider m == 0
00301       delete[] abuffer;
00302       abuffer = 0;
00303       alen = 0;
00304       amax = 0;
00305     }
00306   }
00307 
00308   template <class T>
00309   void Array<T>::copy(const Array<T>& a) {
00310     setmax(a.amax);
00311     alen = a.alen;
00312     for (int i = 0;  i < alen;  i++) {
00313       abuffer[i] = a.abuffer[i];
00314     }
00315   }
00316 
00317   // swap arrays without duplicating memory buffer
00318   template <class T>
00319   void swap(Array<T>& s, Array<T>& t) {
00320     T *tmpbuffer = s.abuffer;  s.abuffer = t.abuffer;  t.abuffer = tmpbuffer;
00321     tmpbuffer = 0;
00322     int tmpn = s.alen;  s.alen = t.alen;  t.alen = tmpn;
00323     tmpn = s.amax;  s.amax = t.amax;  t.amax = tmpn;
00324     tmpn = s.astate;  s.astate = t.astate;  t.astate = tmpn;
00325   }
00326 
00327 
00329   //
00330   // Priority queue for static load balancing of work
00331   //
00333 
00334   // smallest value has priority
00335   // T must have partial ordering operator<=() defined along with assignment
00336   template <class T>
00337   class PriorityQueue {
00338     public:
00339       PriorityQueue(int nelems=0) {
00340         if (nelems > 0)  init(nelems);
00341       }
00342       void init(int nelems) {
00343         a.resize(nelems);  // pre-allocate space
00344         a.resize(0);       // nothing stored yet (does not free memory)
00345       }
00346       void clear() {
00347         a.resize(0);
00348       }
00349       void insert(const T& t) {
00350         a.append(t);
00351         upheap();
00352       }
00353       void remove(T& t) {
00354         int last = a.len() - 1;
00355         if (last < 0) return;
00356         t = a[0];
00357         if (last > 0) a[0] = a[last];
00358         a.resize(last);  // remove last element from array
00359         downheap();
00360       }
00361     private:
00362       // bubble up last element to a correct position
00363       void upheap() {
00364         int n = a.len() - 1;
00365         while (n > 0) {
00366           int parent = (n-1) / 2;
00367           if (a[parent] <= a[n]) break;
00368           T tmp = a[parent];
00369           a[parent] = a[n];
00370           a[n] = tmp;
00371           n = parent;
00372         }
00373       }
00374       // trickle down first element to a correct position
00375       void downheap() {
00376         int n = 0;
00377         int len = a.len();
00378         int left = 2*n + 1;
00379         int right = left + 1;
00380         while (left < len) {
00381           if (right < len && a[right] <= a[left]) {
00382             if (a[n] <= a[right]) break;
00383             T tmp = a[right];
00384             a[right] = a[n];
00385             a[n] = tmp;
00386             n = right;
00387           }
00388           else {
00389             if (a[n] <= a[left]) break;
00390             T tmp = a[left];
00391             a[left] = a[n];
00392             a[n] = tmp;
00393             n = left;
00394           }
00395           left = 2*n + 1;
00396           right = left + 1;
00397         }
00398       }
00399       Array<T> a;
00400   };
00401 
00402 
00404   //
00405   // Grid is 3D lattice of grid points with user-definable index ranges.
00406   //
00408 
00409   // 3-integer vector, used for indexing from a 3D grid
00410   struct Ivec {
00411     int i, j, k;
00412     Ivec(int n=0) : i(n), j(n), k(n) { }
00413     Ivec(int ni, int nj, int nk) : i(ni), j(nj), k(nk) { }
00414     int operator==(const Ivec& n) { return(i==n.i && j==n.j && k==n.k); }
00415 #ifdef MSM_MIGRATION
00416     void pup(PUP::er& p) {
00417       p|i, p|j, p|k;
00418     }
00419 #endif
00420   };
00421 
00422   // index range for 3D lattice of grid points
00423   class IndexRange {
00424     public:
00425       IndexRange() : nlower(0), nextent(0) { }
00426       void set(int pia, int pni, int pja, int pnj, int pka, int pnk) {
00427         ASSERT(pni >= 0 && pnj >= 0 && pnk >= 0);
00428         nlower = Ivec(pia, pja, pka);
00429         nextent = Ivec(pni, pnj, pnk);
00430       }
00431       void setbounds(int pia, int pib, int pja, int pjb, int pka, int pkb) {
00432         set(pia, pib-pia+1, pja, pjb-pja+1, pka, pkb-pka+1);
00433       }
00434       int ia() const { return nlower.i; }
00435       int ib() const { return nlower.i + nextent.i - 1; }
00436       int ja() const { return nlower.j; }
00437       int jb() const { return nlower.j + nextent.j - 1; }
00438       int ka() const { return nlower.k; }
00439       int kb() const { return nlower.k + nextent.k - 1; }
00440       int ni() const { return nextent.i; }
00441       int nj() const { return nextent.j; }
00442       int nk() const { return nextent.k; }
00443       int nn() const { return nextent.i * nextent.j * nextent.k; }
00444       Ivec lower() const { return nlower; }
00445       Ivec extent() const { return nextent; }
00446       int operator<=(const IndexRange& n) {
00447         // true if this IndexRange fits inside n
00448         return ( ia() >= n.ia() && ib() <= n.ib() &&
00449                  ja() >= n.ja() && jb() <= n.jb() &&
00450                  ka() >= n.ka() && kb() <= n.kb() );
00451       }
00452 #ifdef MSM_MIGRATION
00453       void pup(PUP::er& p) {
00454         p|nlower, p|nextent;
00455       }
00456 #endif
00457     protected:
00458       Ivec nlower;   // index for lowest corner of rectangular lattice
00459       Ivec nextent;  // extent of lattice along each dimension
00460   };
00461 
00462   // storage and indexing for 3D lattice of grid points
00463   // with fixed buffer storage no larger than size of block
00464   template <class T> class Grid;
00465 
00466   template <class T, int N>
00467   class GridFixed : public IndexRange {
00468     friend class Grid<T>;
00469     public:
00470       GridFixed() { }
00471       void init(const IndexRange& n) {
00472         nlower = n.lower();
00473         nextent = n.extent();
00474         ASSERT(nextent.i * nextent.j * nextent.k <= N);
00475       }
00476       void set(int pia, int pni, int pja, int pnj, int pka, int pnk) {
00477         IndexRange::set(pia, pni, pja, pnj, pka, pnk);
00478         ASSERT(nextent.i * nextent.j * nextent.k <= N);
00479       }
00480       void setbounds(int pia, int pib, int pja, int pjb, int pka, int pkb) {
00481         IndexRange::setbounds(pia, pib, pja, pjb, pka, pkb);
00482         ASSERT(nextent.i * nextent.j * nextent.k <= N);
00483       }
00484       const T& operator()(int i, int j, int k) const {
00485 #ifdef DEBUG_MSM
00486         return elem(i,j,k);
00487 #else
00488         return gdata[flatindex(i,j,k)];
00489 #endif
00490       }
00491       const T& operator()(const Ivec& n) const {
00492         return this->operator()(n.i, n.j, n.k);
00493       }
00494       const T& elem(int i, int j, int k) const {
00495         if (i<ia() || i>ib() || j<ja() || j>jb() || k<ka() || k>kb()) {
00496           char msg[200];
00497           snprintf(msg, sizeof(msg), "Grid indexing:\n"
00498               "ia=%d, ib=%d, i=%d\n"
00499               "ja=%d, jb=%d, j=%d\n"
00500               "ka=%d, kb=%d, k=%d\n",
00501               ia(), ib(), i, ja(), jb(), j, ka(), kb(), k);
00502           NAMD_die(msg);
00503         }
00504         return gdata[flatindex(i,j,k)];
00505       }
00506       T& operator()(int i, int j, int k) {
00507 #ifdef DEBUG_MSM
00508         return elem(i,j,k);
00509 #else
00510         return gdata[flatindex(i,j,k)];
00511 #endif
00512       }
00513       T& operator()(const Ivec& n) {
00514         return this->operator()(n.i, n.j, n.k);
00515       }
00516       T& elem(int i, int j, int k) {
00517         if (i<ia() || i>ib() || j<ja() || j>jb() || k<ka() || k>kb()) {
00518           char msg[200];
00519           snprintf(msg, sizeof(msg), "Grid indexing:\n"
00520               "ia=%d, ib=%d, i=%d\n"
00521               "ja=%d, jb=%d, j=%d\n"
00522               "ka=%d, kb=%d, k=%d\n",
00523               ia(), ib(), i, ja(), jb(), j, ka(), kb(), k);
00524           NAMD_die(msg);
00525         }
00526         return gdata[flatindex(i,j,k)];
00527       }
00528       int flatindex(int i, int j, int k) const {
00529         return ((k-ka())*nj() + (j-ja()))*ni() + (i-ia());
00530       }
00531       const T *buffer() const { return gdata; }
00532       T *buffer() { return gdata; }
00533 
00534       // use to zero out grid
00535       void reset(const T& t) {
00536         int len = nn();
00537         for (int n = 0;  n < len;  n++) { gdata[n] = t; }
00538       }
00539 
00540       // use to modify the indexing by changing lower corner
00541       void updateLower(const Ivec& n) { nlower = n; }
00542 
00543       // accumulate another grid into this grid
00544       // the grid to be added must fit within this grid's index range
00545       GridFixed<T,N>& operator+=(const GridFixed<T,N>& g) {
00546         ASSERT(IndexRange(g) <= IndexRange(*this));
00547         int gni = g.nextent.i;
00548         int gnj = g.nextent.j;
00549         int gnk = g.nextent.k;
00550         int index = 0;
00551         int ni = nextent.i;
00552         int nij = nextent.i * nextent.j;
00553         int koff = (g.nlower.k - nlower.k) * nij
00554           + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
00555         const T *gbuf = g.gdata.buffer();
00556         T *buf = gdata.buffer();
00557         for (int k = 0;  k < gnk;  k++) {
00558           int jkoff = k * nij + koff;
00559           for (int j = 0;  j < gnj;  j++) {
00560             int ijkoff = j * ni + jkoff;
00561             for (int i = 0;  i < gni;  i++, index++) {
00562               buf[i + ijkoff] += gbuf[index];
00563             }
00564           }
00565         }
00566         return(*this);
00567       }
00568 
00569       // extract a subgrid from this grid
00570       // subgrid must fit within this grid's index range
00571       void extract(GridFixed<T,N>& g) {
00572         ASSERT(IndexRange(g) <= IndexRange(*this));
00573         int gni = g.nextent.i;
00574         int gnj = g.nextent.j;
00575         int gnk = g.nextent.k;
00576         int index = 0;
00577         int ni = nextent.i;
00578         int nij = nextent.i * nextent.j;
00579         int koff = (g.nlower.k - nlower.k) * nij
00580           + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
00581         T *gbuf = g.gdata.buffer();
00582         const T *buf = gdata.buffer();
00583         for (int k = 0;  k < gnk;  k++) {
00584           int jkoff = k * nij + koff;
00585           for (int j = 0;  j < gnj;  j++) {
00586             int ijkoff = j * ni + jkoff;
00587             for (int i = 0;  i < gni;  i++, index++) {
00588               gbuf[index] = buf[i + ijkoff];
00589             }
00590           }
00591         }
00592       }
00593 
00594     private:
00595       T gdata[N];
00596   };
00597 
00598   // storage and indexing for 3D lattice of grid points
00599   template <class T>
00600   class Grid : public IndexRange {
00601     public:
00602       Grid() { }
00603       void init(const IndexRange& n) {
00604         nlower = n.lower();
00605         nextent = n.extent();
00606         gdata.resize(nn());
00607       }
00608       void set(int pia, int pni, int pja, int pnj, int pka, int pnk) {
00609         IndexRange::set(pia, pni, pja, pnj, pka, pnk);
00610         gdata.resize(nn());
00611       }
00612       void setbounds(int pia, int pib, int pja, int pjb, int pka, int pkb) {
00613         IndexRange::setbounds(pia, pib, pja, pjb, pka, pkb);
00614         gdata.resize(nn());
00615       }
00616       void resize(int n) { // reserve space but don't set grid indexing
00617         gdata.resize(n);
00618       }
00619       const T& operator()(int i, int j, int k) const {
00620 #ifdef DEBUG_MSM
00621         return elem(i,j,k);
00622 #else
00623         return gdata[flatindex(i,j,k)];
00624 #endif
00625       }
00626       const T& operator()(const Ivec& n) const {
00627         return this->operator()(n.i, n.j, n.k);
00628       }
00629       const T& elem(int i, int j, int k) const {
00630         if (i<ia() || i>ib() || j<ja() || j>jb() || k<ka() || k>kb()) {
00631           char msg[200];
00632           snprintf(msg, sizeof(msg), "Grid indexing:\n"
00633               "ia=%d, ib=%d, i=%d\n"
00634               "ja=%d, jb=%d, j=%d\n"
00635               "ka=%d, kb=%d, k=%d\n",
00636               ia(), ib(), i, ja(), jb(), j, ka(), kb(), k);
00637           NAMD_die(msg);
00638         }
00639         return gdata[flatindex(i,j,k)];
00640       }
00641       T& operator()(int i, int j, int k) {
00642 #ifdef DEBUG_MSM
00643         return elem(i,j,k);
00644 #else
00645         return gdata[flatindex(i,j,k)];
00646 #endif
00647       }
00648       T& operator()(const Ivec& n) {
00649         return this->operator()(n.i, n.j, n.k);
00650       }
00651       T& elem(int i, int j, int k) {
00652         if (i<ia() || i>ib() || j<ja() || j>jb() || k<ka() || k>kb()) {
00653           char msg[200];
00654           snprintf(msg, sizeof(msg), "Grid indexing:\n"
00655               "ia=%d, ib=%d, i=%d\n"
00656               "ja=%d, jb=%d, j=%d\n"
00657               "ka=%d, kb=%d, k=%d\n",
00658               ia(), ib(), i, ja(), jb(), j, ka(), kb(), k);
00659           NAMD_die(msg);
00660         }
00661         return gdata[flatindex(i,j,k)];
00662       }
00663       int flatindex(int i, int j, int k) const {
00664         return ((k-ka())*nj() + (j-ja()))*ni() + (i-ia());
00665       }
00666       const Array<T>& data() const { return gdata; }
00667       Array<T>& data() { return gdata; }
00668 
00669       // use to zero out grid
00670       void reset(const T& t) {
00671         T *buf = gdata.buffer();
00672         int len = nn();
00673         for (int n = 0;  n < len;  n++) { buf[n] = t; }
00674       }
00675 
00676       // use to modify the indexing by changing lower corner
00677       void updateLower(const Ivec& n) { nlower = n; }
00678 
00679       // accumulate another grid into this grid
00680       // the grid to be added must fit within this grid's index range
00681       Grid<T>& operator+=(const Grid<T>& g) {
00682 #if 1
00683         ASSERT(IndexRange(g) <= IndexRange(*this));
00684 #else
00685         if ( ! (IndexRange(g) <= IndexRange(*this)) ) {
00686           Grid<T> tmp = *this;
00687           // expand myself to hold sum
00688           int ia = nlower.i;
00689           if (ia > g.nlower.i) ia = g.nlower.i;
00690           int ja = nlower.j;
00691           if (ja > g.nlower.j) ja = g.nlower.j;
00692           int ka = nlower.k;
00693           if (ka > g.nlower.k) ka = g.nlower.k;
00694           int ib1 = nlower.i + nextent.i;
00695           int gib1 = g.nlower.i + g.nextent.i;
00696           if (ib1 < gib1) ib1 = gib1;
00697           int jb1 = nlower.j + nextent.j;
00698           int gjb1 = g.nlower.j + g.nextent.j;
00699           if (jb1 < gjb1) jb1 = gjb1;
00700           int kb1 = nlower.k + nextent.k;
00701           int gkb1 = g.nlower.k + g.nextent.k;
00702           if (kb1 < gkb1) kb1 = gkb1;
00703           setbounds(ia, ib1-1, ja, jb1-1, ka, kb1-1);
00704           reset(0);  // make sure constructor for T accepts "0" as its zero
00705           // now copy "tmp" grid elements into my expanded self
00706           int index = 0;
00707           int ni = nextent.i;
00708           int nij = nextent.i * nextent.j;
00709           int koff = (tmp.nlower.k - nlower.k) * nij
00710             + (tmp.nlower.j - nlower.j) * ni + (tmp.nlower.i - nlower.i);
00711           const T *gbuf = tmp.gdata.buffer();
00712           T *buf = gdata.buffer();
00713           for (int k = 0;  k < tmp.nextent.k;  k++) {
00714             int jkoff = k * nij + koff;
00715             for (int j = 0;  j < tmp.nextent.j;  j++) {
00716               int ijkoff = j * ni + jkoff;
00717               for (int i = 0;  i < tmp.nextent.i;  i++, index++) {
00718                 buf[i + ijkoff] = gbuf[index];
00719               }
00720             }
00721           }
00722         }
00723 #endif
00724         int gni = g.nextent.i;
00725         int gnj = g.nextent.j;
00726         int gnk = g.nextent.k;
00727         int index = 0;
00728         int ni = nextent.i;
00729         int nij = nextent.i * nextent.j;
00730         int koff = (g.nlower.k - nlower.k) * nij
00731           + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
00732         const T *gbuf = g.gdata.buffer();
00733         T *buf = gdata.buffer();
00734         for (int k = 0;  k < gnk;  k++) {
00735           int jkoff = k * nij + koff;
00736           for (int j = 0;  j < gnj;  j++) {
00737             int ijkoff = j * ni + jkoff;
00738             for (int i = 0;  i < gni;  i++, index++) {
00739               buf[i + ijkoff] += gbuf[index];
00740             }
00741           }
00742         }
00743         return(*this);
00744       }
00745 
00746       // extract a subgrid from this grid
00747       // subgrid must fit within this grid's index range
00748       void extract(Grid<T>& g) {
00749         ASSERT(IndexRange(g) <= IndexRange(*this));
00750         int gni = g.nextent.i;
00751         int gnj = g.nextent.j;
00752         int gnk = g.nextent.k;
00753         int index = 0;
00754         int ni = nextent.i;
00755         int nij = nextent.i * nextent.j;
00756         int koff = (g.nlower.k - nlower.k) * nij
00757           + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
00758         T *gbuf = g.gdata.buffer();
00759         const T *buf = gdata.buffer();
00760         for (int k = 0;  k < gnk;  k++) {
00761           int jkoff = k * nij + koff;
00762           for (int j = 0;  j < gnj;  j++) {
00763             int ijkoff = j * ni + jkoff;
00764             for (int i = 0;  i < gni;  i++, index++) {
00765               gbuf[index] = buf[i + ijkoff];
00766             }
00767           }
00768         }
00769       }
00770 
00771       // accumulate a fixed size grid into this grid
00772       // the grid to be added must fit within this grid's index range
00773       template <int N>
00774       Grid<T>& operator+=(const GridFixed<T,N>& g) {
00775         ASSERT(IndexRange(g) <= IndexRange(*this));
00776         int gni = g.nextent.i;
00777         int gnj = g.nextent.j;
00778         int gnk = g.nextent.k;
00779         int index = 0;
00780         int ni = nextent.i;
00781         int nij = nextent.i * nextent.j;
00782         int koff = (g.nlower.k - nlower.k) * nij
00783           + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
00784         const T *gbuf = g.buffer();
00785         T *buf = gdata.buffer();
00786         for (int k = 0;  k < gnk;  k++) {
00787           int jkoff = k * nij + koff;
00788           for (int j = 0;  j < gnj;  j++) {
00789             int ijkoff = j * ni + jkoff;
00790             for (int i = 0;  i < gni;  i++, index++) {
00791               buf[i + ijkoff] += gbuf[index];
00792             }
00793           }
00794         }
00795         return(*this);
00796       }
00797 
00798       // extract a subgrid from this grid
00799       // subgrid must fit within this grid's index range
00800       template <int N>
00801       void extract(GridFixed<T,N>& g) {
00802         ASSERT(IndexRange(g) <= IndexRange(*this));
00803         int gni = g.nextent.i;
00804         int gnj = g.nextent.j;
00805         int gnk = g.nextent.k;
00806         int index = 0;
00807         int ni = nextent.i;
00808         int nij = nextent.i * nextent.j;
00809         int koff = (g.nlower.k - nlower.k) * nij
00810           + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
00811         T *gbuf = g.buffer();
00812         const T *buf = gdata.buffer();
00813         for (int k = 0;  k < gnk;  k++) {
00814           int jkoff = k * nij + koff;
00815           for (int j = 0;  j < gnj;  j++) {
00816             int ijkoff = j * ni + jkoff;
00817             for (int i = 0;  i < gni;  i++, index++) {
00818               gbuf[index] = buf[i + ijkoff];
00819             }
00820           }
00821         }
00822       }
00823 
00824     private:
00825       Array<T> gdata;
00826   }; // Grid
00827 
00828 
00830   //
00831   // Map object 
00832   //
00834 
00835   // index a block from the MSM grid hierarchy
00836   struct BlockIndex {
00837     int level;
00838     Ivec n;
00839     BlockIndex() : level(0), n(0) { }
00840     BlockIndex(int ll, const Ivec& nn) : level(ll), n(nn) { }
00841 #ifdef MSM_MIGRATION
00842     void pup(PUP::er& p) {
00843       p|level, p|n;
00844     }
00845 #endif
00846   };
00847 
00848   // for uppermost levels of hierarchy
00849   // fold out image charges along periodic boundaries
00850   // to fill up desired block size
00851   struct FoldFactor {
00852     int active;   // is some numrep dimension > 1?
00853     Ivec numrep;  // number of replications along each dimension
00854     FoldFactor() : active(0), numrep(1) { }
00855     FoldFactor(int i, int j, int k) { set(i,j,k); }
00856     void set(int i, int j, int k) {
00857       if (i <= 0) i = 1;
00858       if (j <= 0) j = 1;
00859       if (k <= 0) k = 1;
00860       if (i > 1 || j > 1 || k > 1) active = 1;
00861       numrep = Ivec(i, j, k);
00862     }
00863   };
00864 
00865   // sending part of an extended grid calculation to another block
00866   struct BlockSend {
00867     BlockIndex nblock;       // relative block index
00868     IndexRange nrange;       // relative grid index range
00869     BlockIndex nblock_wrap;  // true block index
00870     IndexRange nrange_wrap;  // true grid index range
00871     void reset() {
00872       nblock = BlockIndex();
00873       nrange = IndexRange();
00874       nblock_wrap = BlockIndex();
00875       nrange_wrap = IndexRange();
00876     } // reset
00877 #ifdef MSM_MIGRATION
00878     void pup(PUP::er& p) {
00879       p|nblock, p|nrange, p|nblock_wrap, p|nrange_wrap;
00880     }
00881 #endif
00882   };
00883 
00884   struct PatchSend {
00885     IndexRange nrange;         // true grid index range from my block
00886     IndexRange nrange_unwrap;  // relative grid index range for patch
00887     int patchID;               // patch ID
00888     void reset() {
00889       nrange = IndexRange();
00890       nrange_unwrap = IndexRange();
00891       patchID = -1;
00892     } // reset
00893   };
00894 
00895   // one PatchDiagram for each patch
00896   // maintain a Grid of PatchDiagram, indexed by patch ID
00897   struct PatchDiagram {
00898     IndexRange nrange;       // shows subset of MSM h-grid covering this patch
00899     Array<BlockSend> send;   // array of blocks to which this patch sends
00900     int numRecvs;            // number of blocks from which this patch receives
00901     void reset() {
00902       nrange = IndexRange();
00903       send.resize(0);
00904       numRecvs = 0;
00905     } // reset
00906   };
00907 
00908   // one BlockDiagram for each block of each level of each MSM grid
00909   // maintain a Grid of BlockDiagram for each level
00910   struct BlockDiagram {
00911     IndexRange nrange;            // subset of MSM grid for this block
00912     IndexRange nrangeCutoff;      // expanded subgrid for cutoff calculation
00913     IndexRange nrangeRestricted;  // (level+1) subgrid for restriction
00914     IndexRange nrangeProlongated; // (level-1) subgrid for prolongation
00915     Array<BlockSend> sendUp;      // send up charge to blocks on (level+1)
00916     Array<BlockSend> sendAcross;  // send across potential to blocks on (level)
00917     Array<int> indexGridCutoff;   // index of MsmGridCutoff chares to calculate
00918                                   // each charge -> potential block interaction
00919     Array<int> recvGridCutoff;    // index of MsmGridCutoff chares contributing
00920                                   // back into my potential block
00921     Array<BlockSend> sendDown;    // send down potential to blocks on (level-1)
00922     Array<PatchSend> sendPatch;   // send my (level=0) potential block to patch
00923     int numRecvsCharge;           // number of expected receives of charge
00924     int numRecvsPotential;        // number of expected receives of potential
00925 
00926     void reset() {
00927       nrange = IndexRange();
00928       nrangeCutoff = IndexRange();
00929       nrangeRestricted = IndexRange();
00930       nrangeProlongated = IndexRange();
00931       sendUp.resize(0);
00932       sendAcross.resize(0);
00933       sendDown.resize(0);
00934       sendPatch.resize(0);
00935       numRecvsCharge = 0;
00936       numRecvsPotential = 0;
00937     } // reset
00938   };
00939 
00940 
00941   struct Map {
00942     Array<IndexRange> gridrange;  // dimensions for each MSM grid level
00943 
00944     Array<Grid<Float> > gc;       // grid constant weights for each level
00945     Array<Grid<Float> > gvc;      // virial grid weights for each level
00946     Grid<Float> grespro;          // restriction / prolongation nonzero stencil
00947                                   // requires correct IndexOffset array
00948 
00949     Array<Grid<C1Matrix> > gc_c1hermite;    // grid constant weights C1 Hermite
00950     Array<Grid<C1Matrix> > gvc_c1hermite;   // virial grid weights C1 Hermite
00951     Array<Grid<C1Matrix> > gres_c1hermite;  // restriction stencil C1 Hermite
00952     Array<Grid<C1Matrix> > gpro_c1hermite;  // prolongation stencil C1 Hermite
00953                                             // requires index offsets
00954 
00955     Array<PatchDiagram> patchList;
00956     Array<Grid<BlockDiagram> > blockLevel;
00957 
00958     int ispx, ispy, ispz;         // is periodic in x, y, z?
00959 
00960     Array<int> bsx, bsy, bsz;     // block size in x, y, z for each level
00961 
00962     Array<FoldFactor> foldfactor; // for uppermost grid levels
00963       // replicate periodic dimensions in order to fill up block size
00964 
00965     // clip index to grid level, using periodicity flags
00966     Ivec clipIndexToLevel(const Ivec& n, int level) const {
00967       ASSERT(level >= 0 && level < gridrange.len());
00968       Ivec pn(n);
00969       if ( ! ispx) {
00970         int a = gridrange[level].ia();
00971         int b = gridrange[level].ib();
00972         if (pn.i < a) pn.i = a;
00973         if (pn.i > b) pn.i = b;
00974       }
00975       if ( ! ispy) {
00976         int a = gridrange[level].ja();
00977         int b = gridrange[level].jb();
00978         if (pn.j < a) pn.j = a;
00979         if (pn.j > b) pn.j = b;
00980       }
00981       if ( ! ispz) {
00982         int a = gridrange[level].ka();
00983         int b = gridrange[level].kb();
00984         if (pn.k < a) pn.k = a;
00985         if (pn.k > b) pn.k = b;
00986       }
00987       return pn;
00988     }
00989 
00990     // determine relative (unwrapped) block index for the given grid index
00991     BlockIndex blockOfGridIndex(const Ivec& n, int level) const {
00992       ASSERT(level >= 0 && level < gridrange.len());
00993       BlockIndex bn;
00994       // we want floor((i - ia) / bsx), etc.
00995       // modify case i < ia to avoid integer division of negative numbers
00996       int d = n.i - gridrange[level].ia();
00997       bn.n.i = (d >= 0 ? d / bsx[level] : -((-d+bsx[level]-1) / bsx[level]));
00998       d = n.j - gridrange[level].ja();
00999       bn.n.j = (d >= 0 ? d / bsy[level] : -((-d+bsy[level]-1) / bsy[level]));
01000       d = n.k - gridrange[level].ka();
01001       bn.n.k = (d >= 0 ? d / bsz[level] : -((-d+bsz[level]-1) / bsz[level]));
01002       bn.level = level;
01003       return bn;
01004     }
01005 
01006     // determine relative (unwrapped) block index for the given grid index
01007     // for unfolded replication of image charges
01008     BlockIndex blockOfGridIndexFold(const Ivec& n, int level) const {
01009       ASSERT(level >= 0 && level < gridrange.len());
01010       BlockIndex bn;
01011       int bsi = foldfactor[level].numrep.i * bsx[level];
01012       int bsj = foldfactor[level].numrep.j * bsy[level];
01013       int bsk = foldfactor[level].numrep.k * bsz[level];
01014       // we want floor((i - ia) / bsx), etc.
01015       // modify case i < ia to avoid integer division of negative numbers
01016       int d = n.i - gridrange[level].ia();
01017       bn.n.i = (d >= 0 ? d / bsi : -((-d+bsi-1) / bsi));
01018       d = n.j - gridrange[level].ja();
01019       bn.n.j = (d >= 0 ? d / bsj : -((-d+bsj-1) / bsj));
01020       d = n.k - gridrange[level].ka();
01021       bn.n.k = (d >= 0 ? d / bsk : -((-d+bsk-1) / bsk));
01022       bn.level = level;
01023       return bn;
01024     }
01025 
01026     // find the natural index range of the given relative block number
01027     IndexRange indexRangeOfBlock(const BlockIndex& nb) const {
01028       ASSERT(nb.level >= 0 && nb.level < gridrange.len());
01029       IndexRange nr;
01030       int ia = nb.n.i * bsx[nb.level] + gridrange[nb.level].ia();
01031       int ja = nb.n.j * bsy[nb.level] + gridrange[nb.level].ja();
01032       int ka = nb.n.k * bsz[nb.level] + gridrange[nb.level].ka();
01033       nr.set(ia, bsx[nb.level], ja, bsy[nb.level], ka, bsz[nb.level]);
01034       return nr;
01035     }
01036 
01037     // find the natural index range of the given relative block number
01038     // for unfolded replication of image charges
01039     IndexRange indexRangeOfBlockFold(const BlockIndex& nb) const {
01040       ASSERT(nb.level >= 0 && nb.level < gridrange.len());
01041       int bsi = foldfactor[nb.level].numrep.i * bsx[nb.level];
01042       int bsj = foldfactor[nb.level].numrep.j * bsy[nb.level];
01043       int bsk = foldfactor[nb.level].numrep.k * bsz[nb.level];
01044       IndexRange nr;
01045       int ia = nb.n.i * bsi + gridrange[nb.level].ia();
01046       int ja = nb.n.j * bsj + gridrange[nb.level].ja();
01047       int ka = nb.n.k * bsk + gridrange[nb.level].ka();
01048       nr.set(ia, bsi, ja, bsj, ka, bsk);
01049       return nr;
01050     }
01051 
01052     // clip the natural block index range to not exceed the given index range
01053     IndexRange clipBlockToIndexRange(const BlockIndex& nb,
01054         const IndexRange& nrange) const {
01055       IndexRange nr = indexRangeOfBlock(nb);
01056       int nia = nrange.ia();
01057       int nib = nrange.ib();
01058       int nja = nrange.ja();
01059       int njb = nrange.jb();
01060       int nka = nrange.ka();
01061       int nkb = nrange.kb();
01062       int ia = nr.ia();
01063       if (ia < nia) ia = nia;
01064       int ib = nr.ib();
01065       if (ib > nib) ib = nib;
01066       int ja = nr.ja();
01067       if (ja < nja) ja = nja;
01068       int jb = nr.jb();
01069       if (jb > njb) jb = njb;
01070       int ka = nr.ka();
01071       if (ka < nka) ka = nka;
01072       int kb = nr.kb();
01073       if (kb > nkb) kb = nkb;
01074       nr.setbounds(ia, ib, ja, jb, ka, kb);
01075       return nr;
01076     }
01077 
01078     // clip the natural block index range to not exceed the given index range
01079     // for unfolded replication of image charges
01080     IndexRange clipBlockToIndexRangeFold(const BlockIndex& nb,
01081         const IndexRange& nrange) const {
01082       IndexRange nr = indexRangeOfBlockFold(nb);
01083       int nia = nrange.ia();
01084       int nib = nrange.ib();
01085       int nja = nrange.ja();
01086       int njb = nrange.jb();
01087       int nka = nrange.ka();
01088       int nkb = nrange.kb();
01089       int ia = nr.ia();
01090       if (ia < nia) ia = nia;
01091       int ib = nr.ib();
01092       if (ib > nib) ib = nib;
01093       int ja = nr.ja();
01094       if (ja < nja) ja = nja;
01095       int jb = nr.jb();
01096       if (jb > njb) jb = njb;
01097       int ka = nr.ka();
01098       if (ka < nka) ka = nka;
01099       int kb = nr.kb();
01100       if (kb > nkb) kb = nkb;
01101       nr.setbounds(ia, ib, ja, jb, ka, kb);
01102       return nr;
01103     }
01104 
01105     // set the nblock_wrap and nrange_wrap fields based on periodicity
01106     void wrapBlockSend(BlockSend& bs) const {
01107       BlockIndex nb = bs.nblock;
01108       IndexRange nr = bs.nrange;
01109       int level = bs.nblock.level;
01110       ASSERT(level >= 0 && level < blockLevel.len());
01111       int ni = blockLevel[level].ni();
01112       int nj = blockLevel[level].nj();
01113       int nk = blockLevel[level].nk();
01114       int di=0, dj=0, dk=0;
01115       if (ispx) {
01116         while (nb.n.i < 0) {
01117           nb.n.i += ni;
01118           di += ni * bsx[level];
01119         }
01120         while (nb.n.i >= ni) {
01121           nb.n.i -= ni;
01122           di -= ni * bsx[level];
01123         }
01124       }
01125       if (ispy) {
01126         while (nb.n.j < 0) {
01127           nb.n.j += nj;
01128           dj += nj * bsy[level];
01129         }
01130         while (nb.n.j >= nj) {
01131           nb.n.j -= nj;
01132           dj -= nj * bsy[level];
01133         }
01134       }
01135       if (ispz) {
01136         while (nb.n.k < 0) {
01137           nb.n.k += nk;
01138           dk += nk * bsz[level];
01139         }
01140         while (nb.n.k >= nk) {
01141           nb.n.k -= nk;
01142           dk -= nk * bsz[level];
01143         }
01144       }
01145       int ia = nr.ia();
01146       int ib = nr.ib();
01147       int ja = nr.ja();
01148       int jb = nr.jb();
01149       int ka = nr.ka();
01150       int kb = nr.kb();
01151       nr.setbounds(ia + di, ib + di, ja + dj, jb + dj, ka + dk, kb + dk);
01152       bs.nblock_wrap = nb;
01153       bs.nrange_wrap = nr;
01154     }
01155 
01156     // set the nblock_wrap and nrange_wrap fields based on periodicity
01157     // for unfolded replication of image charges
01158     void wrapBlockSendFold(BlockSend& bs) const {
01159       BlockIndex nb = bs.nblock;
01160       IndexRange nr = bs.nrange;
01161       int level = bs.nblock.level;
01162       ASSERT(level >= 0 && level < blockLevel.len());
01163       int foldi = foldfactor[level].numrep.i;
01164       int foldj = foldfactor[level].numrep.j;
01165       int foldk = foldfactor[level].numrep.k;
01166       int ni = blockLevel[level].ni();
01167       int nj = blockLevel[level].nj();
01168       int nk = blockLevel[level].nk();
01169       int bsi = foldi * bsx[level];
01170       int bsj = foldj * bsy[level];
01171       int bsk = foldk * bsz[level];
01172       int di=0, dj=0, dk=0;
01173       if (ispx) {
01174         while (nb.n.i < 0) {
01175           nb.n.i += ni;
01176           di += ni * bsi;
01177         }
01178         while (nb.n.i >= ni) {
01179           nb.n.i -= ni;
01180           di -= ni * bsi;
01181         }
01182       }
01183       if (ispy) {
01184         while (nb.n.j < 0) {
01185           nb.n.j += nj;
01186           dj += nj * bsj;
01187         }
01188         while (nb.n.j >= nj) {
01189           nb.n.j -= nj;
01190           dj -= nj * bsj;
01191         }
01192       }
01193       if (ispz) {
01194         while (nb.n.k < 0) {
01195           nb.n.k += nk;
01196           dk += nk * bsk;
01197         }
01198         while (nb.n.k >= nk) {
01199           nb.n.k -= nk;
01200           dk -= nk * bsk;
01201         }
01202       }
01203       int ia = nr.ia();
01204       int ib = nr.ib();
01205       int ja = nr.ja();
01206       int jb = nr.jb();
01207       int ka = nr.ka();
01208       int kb = nr.kb();
01209       nr.setbounds(ia + di, ib + di, ja + dj, jb + dj, ka + dk, kb + dk);
01210       bs.nblock_wrap = nb;
01211       bs.nrange_wrap = nr;
01212     }
01213 
01214     void wrapBlockIndex(BlockIndex& bn) const {
01215       int level = bn.level;
01216       int ni = blockLevel[level].ni();
01217       int nj = blockLevel[level].nj();
01218       int nk = blockLevel[level].nk();
01219       if (ispx) {
01220         while (bn.n.i < 0) {
01221           bn.n.i += ni;
01222         }
01223         while (bn.n.i >= ni) {
01224           bn.n.i -= ni;
01225         }
01226       }
01227       if (ispy) {
01228         while (bn.n.j < 0) {
01229           bn.n.j += nj;
01230         }
01231         while (bn.n.j >= nj) {
01232           bn.n.j -= nj;
01233         }
01234       }
01235       if (ispz) {
01236         while (bn.n.k < 0) {
01237           bn.n.k += nk;
01238         }
01239         while (bn.n.k >= nk) {
01240           bn.n.k -= nk;
01241         }
01242       }
01243     }
01244 
01245   }; // Map
01246 
01247 
01248   struct AtomCoord {
01249     Position position;
01250     Real charge;
01251     int id;
01252   };
01253 
01254   typedef Array<AtomCoord> AtomCoordArray;
01255   typedef Array<Force> ForceArray;
01256 
01257   struct PatchData;
01258   typedef Array<PatchData *> PatchPtrArray;
01259 
01260 } // namespace msm
01261 
01262 #endif // MSMMAP_H

Generated on Thu Sep 21 01:17:13 2017 for NAMD by  doxygen 1.4.7