NAMD
MsmMap.h
Go to the documentation of this file.
1 
7 #ifndef MSMMAP_H
8 #define MSMMAP_H
9 
10 // SSE and AVX vector intrinsics and memory alignment macros
11 #if defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE)
12 #include <emmintrin.h> // SSE2
13 #if defined(__INTEL_COMPILER)
14 #define __align(X) __declspec(align(X) )
15 #elif defined(__PGI)
16 #define __align(X) __attribute__((aligned(X) ))
17 #define MISSING_mm_cvtsd_f64
18 #elif defined(__GNUC__)
19 #define __align(X) __attribute__((aligned(X) ))
20 #if (__GNUC__ < 4)
21 #define MISSING_mm_cvtsd_f64
22 #endif
23 #else
24 #define __align(X) __declspec(align(X) )
25 #endif
26 #endif
27 
28 // migration of MSM computes not currently enabled
29 #define MSM_MIGRATION
30 #undef MSM_MIGRATION
31 
32 // for debugging MSM migration
33 #define DEBUG_MSM_MIGRATE
34 #undef DEBUG_MSM_MIGRATE
35 
36 #define MSM_MAX_BLOCK_SIZE 8
37 #define MSM_MAX_BLOCK_VOLUME \
38  (MSM_MAX_BLOCK_SIZE * MSM_MAX_BLOCK_SIZE * MSM_MAX_BLOCK_SIZE)
39 
40 #define MSM_C1VECTOR_MAX_BLOCK_SIZE (MSM_MAX_BLOCK_SIZE / 2)
41 #define MSM_C1VECTOR_MAX_BLOCK_VOLUME \
42  (MSM_C1VECTOR_MAX_BLOCK_SIZE * \
43  MSM_C1VECTOR_MAX_BLOCK_SIZE * \
44  MSM_C1VECTOR_MAX_BLOCK_SIZE)
45 
46 #define DEBUG_MSM
47 #undef DEBUG_MSM
48 
49 #define DEBUG_MSM_VERBOSE
50 #undef DEBUG_MSM_VERBOSE
51 
52 #define DEBUG_MSM_GRID
53 #undef DEBUG_MSM_GRID
54 
55 // assert macro
56 #undef ASSERT
57 #ifdef DEBUG_MSM
58 #define ASSERT(expr) \
59  do { \
60  if ( !(expr) ) { \
61  char msg[100]; \
62  snprintf(msg, sizeof(msg), "ASSERT: \"%s\" " \
63  "(%s, %d)\n", #expr, __FILE__, __LINE__); \
64  NAMD_die(msg); \
65  } \
66  } while (0)
67 #else
68 #define ASSERT(expr)
69 #endif
70 
71 
72 // employ mixed precision
73 // (but allow easy change to all double precision for comparison)
74 typedef float Float;
75 typedef double Double;
76 
77 
79  //
80  // Vector and matrix elements for C1 Hermite interpolation.
81  //
83 
84  enum { C1_VECTOR_SIZE = 8, C1_MATRIX_SIZE = 8*8 };
85 
86  struct C1Vector {
88  C1Vector(Float r=0) { set(r); }
89  void set(Float r) {
90  for (int n=0; n < C1_VECTOR_SIZE; n++) velem[n] = r;
91  }
93  for (int n=0; n < C1_VECTOR_SIZE; n++) velem[n] += v.velem[n];
94  return(*this);
95  }
96  friend Float operator*(const C1Vector& u, const C1Vector& v) {
97  Float r=0;
98  for (int n=0; n < C1_VECTOR_SIZE; n++) r += u.velem[n] * v.velem[n];
99  return r;
100  }
101  friend C1Vector operator+(const C1Vector& u, const C1Vector& v) {
102  C1Vector w;
103  for (int n=0; n < C1_VECTOR_SIZE; n++) {
104  w.velem[n] = u.velem[n] + v.velem[n];
105  }
106  return w;
107  }
108  };
109 
110  struct C1Matrix {
112  C1Matrix(Float r=0) { set(r); }
113  void set(Float r) {
114  for (int n=0; n < C1_MATRIX_SIZE; n++) melem[n] = 0;
115  }
116  friend C1Vector operator*(const C1Matrix& m, const C1Vector& u) {
117  C1Vector v;
118 
119  // XXX not tested yet
120 #if 1 && (defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE))
121  // Hand-coded SSE2 vectorization
122  // This loop requires that the single-precision input arrays be
123  // aligned on 16-byte boundaries, such that array[index % 4 == 0]
124  // can be safely accessed with aligned load/store operations
125  for (int k=0, j=0; j < C1_VECTOR_SIZE; j++) {
126  __m128 melem4 = _mm_load_ps(&m.melem[k]);
127  __m128 uelem4 = _mm_load_ps(&u.velem[0]);
128  __m128 tmp4 = _mm_mul_ps(melem4, uelem4);
129  melem4 = _mm_load_ps(&m.melem[k+4]);
130  uelem4 = _mm_load_ps(&u.velem[4]);
131  tmp4 = _mm_add_ps(tmp4, _mm_mul_ps(melem4, uelem4));
132 
133  // do a 4-element reduction and accumulate result
134  __m128 sum4 = tmp4;
135  sum4 = _mm_shuffle_ps(sum4, sum4, _MM_SHUFFLE(2, 3, 0, 1));
136  sum4 = _mm_add_ps(sum4, tmp4);
137  tmp4 = sum4;
138  sum4 = _mm_shuffle_ps(sum4, sum4, _MM_SHUFFLE(1, 0, 3, 2));
139  sum4 = _mm_add_ps(sum4, tmp4);
140 
141  // all 4 elements are now set to the sum
142  float sum;
143  _mm_store_ss(&sum, sum4); // store lowest element
144  v.velem[j] += sum;
145  k+=8;
146  }
147 #elif 0 && (defined(__AVX__) && ! defined(NAMD_DISABLE_SSE))
148  // Hand-coded AVX vectorization
149  // This loop requires that the single-precision input arrays be
150  // aligned on 32-byte boundaries, such that array[index % 8 == 0]
151  // can be safely accessed with aligned load/store operations
152  for (int k=0, j=0; j < C1_VECTOR_SIZE; j++) {
153  __m256 melem8 = _mm256_load_ps(&m.melem[k]);
154  __m256 uelem8 = _mm256_load_ps(&u.velem[0]);
155  __m256 tmp8 = _mm256_mul_ps(melem8, uelem8);
156 
157  // XXX this still needs to be rewritten a bit for AVX
158  // do an 8-element reduction and accumulate result
159  __m256 sum8 = tmp8;
160  sum8 = _mm256_hadd_ps(sum8, sum8);
161  sum8 = _mm256_hadd_ps(sum8, sum8);
162  tmp8 = sum8;
163  tmp8 = _mm256_permute2f128_ps(tmp8, tmp8, 1);
164  sum8 = _mm256_hadd_ps(tmp8, sum8);
165 
166  // all 8 elements are now set to the sum
167  float sum;
168  _mm_store_ss(&sum, sum8); // store lowest element
169  v.velem[j] += sum;
170  k+=8;
171  }
172 #else
173 #if defined(__INTEL_COMPILER)
174 #pragma vector always
175 #endif
176  for (int k=0, j=0; j < C1_VECTOR_SIZE; j++) {
177  for (int i = 0; i < C1_VECTOR_SIZE; i++, k++) {
178  v.velem[j] += m.melem[k] * u.velem[i];
179  }
180  }
181 #endif
182  return v;
183  }
184  };
185 
186  // index vector based on mixed partial derivatives in x,y,z
187  enum { D000=0, D100, D010, D001, D110, D101, D011, D111 };
188 
189  // index matrix using 2 vector indexes, row-major column ordering,
190  // defining partial derivatives of g(xj,yj,zj,xi,yi,zi)
191 #define C1INDEX(drj,dri) ((drj)*C1_VECTOR_SIZE + (dri))
192 
193 
194 namespace msm {
195 
197  //
198  // Resizable Array class
199  //
201 
202  template <class T> class Array;
203 
204  template <class T>
205  void swap(Array<T>& s, Array<T>& t);
206 
207  template <class T>
208  class Array {
209  public:
210  Array() : abuffer(0), alen(0), amax(0) { }
211  Array(int n) : abuffer(0), alen(0), amax(0) { resize(n); }
212  Array(const Array& a) : abuffer(0), alen(0), amax(0) { copy(a); }
213  ~Array() { setmax(0); }
214  Array& operator=(const Array& a) {
215  if (this != &a) copy(a); // don't allow self-assignment
216  return(*this);
217  }
218  int len() const { return alen; }
219  int max() const { return amax; }
220  const T& operator[](int i) const {
221 #ifdef DEBUG_MSM
222  return elem(i);
223 #else
224  return abuffer[i];
225 #endif
226  }
227  const T& elem(int i) const {
228  if (i < 0 || i >= alen) {
229  char msg[100];
230  snprintf(msg, sizeof(msg), "Array index: alen=%d, i=%d\n", alen, i);
231  NAMD_die(msg);
232  }
233  return abuffer[i];
234  }
235  T& operator[](int i) {
236 #ifdef DEBUG_MSM
237  return elem(i);
238 #else
239  return abuffer[i];
240 #endif
241  }
242  T& elem(int i) {
243  if (i < 0 || i >= alen) {
244  char msg[100];
245  snprintf(msg, sizeof(msg), "Array index: alen=%d, i=%d\n", alen, i);
246  NAMD_die(msg);
247  }
248  return abuffer[i];
249  }
250  void append(const T& t) {
251  if (alen==amax) setmax(2*amax+1);
252  abuffer[alen++] = t;
253  }
254  void resize(int n) {
255  if (n > amax) setmax(n);
256  alen = n;
257  }
258  void setmax(int m);
259  const T *buffer() const { return abuffer; }
260  T *buffer() { return abuffer; }
261  const T *buffer(int& n) const { n = alen; return abuffer; }
262  T *buffer(int& n) { n = alen; return abuffer; }
263  void reset(const T& t) {
264  for (int n = 0; n < alen; n++) abuffer[n] = t;
265  }
266  friend void swap<T>(Array&, Array&);
267 #ifdef DEBUG_MSM
268  void print(const char *s=0) const {
269  if (s) printf("PRINTING DATA FOR ARRAY \"%s\":\n", s);
270  printf("abuffer=%p\n alen=%d amax=%d\n",
271  (void *) abuffer, alen, amax);
272  }
273 #endif
274  protected:
276  int alen, amax;
277  void copy(const Array& a);
278  };
279 
280  template <class T>
281  void Array<T>::setmax(int m) {
282  if (m == amax) return;
283  else if (m > 0) {
284  T *newbuffer = new T[m];
285  if ( ! newbuffer) {
286  char msg[100];
287  snprintf(msg, sizeof(msg),
288  "Can't allocate %lu KB for msm::Array\n",
289  (unsigned long)(m * sizeof(T) / 1024));
290  NAMD_die(msg);
291  }
292  if (alen > m) alen = m; // new buffer is shorter than old buffer
293  for (int i = 0; i < alen; i++) {
294  newbuffer[i] = abuffer[i];
295  }
296  delete[] abuffer;
297  abuffer = newbuffer;
298  amax = m;
299  }
300  else { // consider m == 0
301  delete[] abuffer;
302  abuffer = 0;
303  alen = 0;
304  amax = 0;
305  }
306  }
307 
308  template <class T>
309  void Array<T>::copy(const Array<T>& a) {
310  setmax(a.amax);
311  alen = a.alen;
312  for (int i = 0; i < alen; i++) {
313  abuffer[i] = a.abuffer[i];
314  }
315  }
316 
317  // swap arrays without duplicating memory buffer
318  template <class T>
319  void swap(Array<T>& s, Array<T>& t) {
320  T *tmpbuffer = s.abuffer; s.abuffer = t.abuffer; t.abuffer = tmpbuffer;
321  tmpbuffer = 0;
322  int tmpn = s.alen; s.alen = t.alen; t.alen = tmpn;
323  tmpn = s.amax; s.amax = t.amax; t.amax = tmpn;
324  tmpn = s.astate; s.astate = t.astate; t.astate = tmpn;
325  }
326 
327 
329  //
330  // Priority queue for static load balancing of work
331  //
333 
334  // smallest value has priority
335  // T must have partial ordering operator<=() defined along with assignment
336  template <class T>
338  public:
339  PriorityQueue(int nelems=0) {
340  if (nelems > 0) init(nelems);
341  }
342  void init(int nelems) {
343  a.resize(nelems); // pre-allocate space
344  a.resize(0); // nothing stored yet (does not free memory)
345  }
346  void clear() {
347  a.resize(0);
348  }
349  void insert(const T& t) {
350  a.append(t);
351  upheap();
352  }
353  void remove(T& t) {
354  int last = a.len() - 1;
355  if (last < 0) return;
356  t = a[0];
357  if (last > 0) a[0] = a[last];
358  a.resize(last); // remove last element from array
359  downheap();
360  }
361  private:
362  // bubble up last element to a correct position
363  void upheap() {
364  int n = a.len() - 1;
365  while (n > 0) {
366  int parent = (n-1) / 2;
367  if (a[parent] <= a[n]) break;
368  T tmp = a[parent];
369  a[parent] = a[n];
370  a[n] = tmp;
371  n = parent;
372  }
373  }
374  // trickle down first element to a correct position
375  void downheap() {
376  int n = 0;
377  int len = a.len();
378  int left = 2*n + 1;
379  int right = left + 1;
380  while (left < len) {
381  if (right < len && a[right] <= a[left]) {
382  if (a[n] <= a[right]) break;
383  T tmp = a[right];
384  a[right] = a[n];
385  a[n] = tmp;
386  n = right;
387  }
388  else {
389  if (a[n] <= a[left]) break;
390  T tmp = a[left];
391  a[left] = a[n];
392  a[n] = tmp;
393  n = left;
394  }
395  left = 2*n + 1;
396  right = left + 1;
397  }
398  }
399  Array<T> a;
400  };
401 
402 
404  //
405  // Grid is 3D lattice of grid points with user-definable index ranges.
406  //
408 
409  // 3-integer vector, used for indexing from a 3D grid
410  struct Ivec {
411  int i, j, k;
412  Ivec(int n=0) : i(n), j(n), k(n) { }
413  Ivec(int ni, int nj, int nk) : i(ni), j(nj), k(nk) { }
414  int operator==(const Ivec& n) { return(i==n.i && j==n.j && k==n.k); }
415 #ifdef MSM_MIGRATION
416  void pup(PUP::er& p) {
417  p|i, p|j, p|k;
418  }
419 #endif
420  };
421 
422  // index range for 3D lattice of grid points
423  class IndexRange {
424  public:
425  IndexRange() : nlower(0), nextent(0) { }
426  void set(int pia, int pni, int pja, int pnj, int pka, int pnk) {
427  ASSERT(pni >= 0 && pnj >= 0 && pnk >= 0);
428  nlower = Ivec(pia, pja, pka);
429  nextent = Ivec(pni, pnj, pnk);
430  }
431  void setbounds(int pia, int pib, int pja, int pjb, int pka, int pkb) {
432  set(pia, pib-pia+1, pja, pjb-pja+1, pka, pkb-pka+1);
433  }
434  int ia() const { return nlower.i; }
435  int ib() const { return nlower.i + nextent.i - 1; }
436  int ja() const { return nlower.j; }
437  int jb() const { return nlower.j + nextent.j - 1; }
438  int ka() const { return nlower.k; }
439  int kb() const { return nlower.k + nextent.k - 1; }
440  int ni() const { return nextent.i; }
441  int nj() const { return nextent.j; }
442  int nk() const { return nextent.k; }
443  int nn() const { return nextent.i * nextent.j * nextent.k; }
444  Ivec lower() const { return nlower; }
445  Ivec extent() const { return nextent; }
446  int operator<=(const IndexRange& n) {
447  // true if this IndexRange fits inside n
448  return ( ia() >= n.ia() && ib() <= n.ib() &&
449  ja() >= n.ja() && jb() <= n.jb() &&
450  ka() >= n.ka() && kb() <= n.kb() );
451  }
452 #ifdef MSM_MIGRATION
453  void pup(PUP::er& p) {
454  p|nlower, p|nextent;
455  }
456 #endif
457  protected:
458  Ivec nlower; // index for lowest corner of rectangular lattice
459  Ivec nextent; // extent of lattice along each dimension
460  };
461 
462  // storage and indexing for 3D lattice of grid points
463  // with fixed buffer storage no larger than size of block
464  template <class T> class Grid;
465 
466  template <class T, int N>
467  class GridFixed : public IndexRange {
468  friend class Grid<T>;
469  public:
470  GridFixed() { }
471  void init(const IndexRange& n) {
472  nlower = n.lower();
473  nextent = n.extent();
474  ASSERT(nextent.i * nextent.j * nextent.k <= N);
475  }
476  void set(int pia, int pni, int pja, int pnj, int pka, int pnk) {
477  IndexRange::set(pia, pni, pja, pnj, pka, pnk);
478  ASSERT(nextent.i * nextent.j * nextent.k <= N);
479  }
480  void setbounds(int pia, int pib, int pja, int pjb, int pka, int pkb) {
481  IndexRange::setbounds(pia, pib, pja, pjb, pka, pkb);
482  ASSERT(nextent.i * nextent.j * nextent.k <= N);
483  }
484  const T& operator()(int i, int j, int k) const {
485 #ifdef DEBUG_MSM
486  return elem(i,j,k);
487 #else
488  return gdata[flatindex(i,j,k)];
489 #endif
490  }
491  const T& operator()(const Ivec& n) const {
492  return this->operator()(n.i, n.j, n.k);
493  }
494  const T& elem(int i, int j, int k) const {
495  if (i<ia() || i>ib() || j<ja() || j>jb() || k<ka() || k>kb()) {
496  char msg[200];
497  snprintf(msg, sizeof(msg), "Grid indexing:\n"
498  "ia=%d, ib=%d, i=%d\n"
499  "ja=%d, jb=%d, j=%d\n"
500  "ka=%d, kb=%d, k=%d\n",
501  ia(), ib(), i, ja(), jb(), j, ka(), kb(), k);
502  NAMD_die(msg);
503  }
504  return gdata[flatindex(i,j,k)];
505  }
506  T& operator()(int i, int j, int k) {
507 #ifdef DEBUG_MSM
508  return elem(i,j,k);
509 #else
510  return gdata[flatindex(i,j,k)];
511 #endif
512  }
513  T& operator()(const Ivec& n) {
514  return this->operator()(n.i, n.j, n.k);
515  }
516  T& elem(int i, int j, int k) {
517  if (i<ia() || i>ib() || j<ja() || j>jb() || k<ka() || k>kb()) {
518  char msg[200];
519  snprintf(msg, sizeof(msg), "Grid indexing:\n"
520  "ia=%d, ib=%d, i=%d\n"
521  "ja=%d, jb=%d, j=%d\n"
522  "ka=%d, kb=%d, k=%d\n",
523  ia(), ib(), i, ja(), jb(), j, ka(), kb(), k);
524  NAMD_die(msg);
525  }
526  return gdata[flatindex(i,j,k)];
527  }
528  int flatindex(int i, int j, int k) const {
529  return ((k-ka())*nj() + (j-ja()))*ni() + (i-ia());
530  }
531  const T *buffer() const { return gdata; }
532  T *buffer() { return gdata; }
533 
534  // use to zero out grid
535  void reset(const T& t) {
536  int len = nn();
537  for (int n = 0; n < len; n++) { gdata[n] = t; }
538  }
539 
540  // use to modify the indexing by changing lower corner
541  void updateLower(const Ivec& n) { nlower = n; }
542 
543  // accumulate another grid into this grid
544  // the grid to be added must fit within this grid's index range
546  ASSERT(IndexRange(g) <= IndexRange(*this));
547  int gni = g.nextent.i;
548  int gnj = g.nextent.j;
549  int gnk = g.nextent.k;
550  int index = 0;
551  int ni = nextent.i;
552  int nij = nextent.i * nextent.j;
553  int koff = (g.nlower.k - nlower.k) * nij
554  + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
555  const T *gbuf = g.gdata.buffer();
556  T *buf = gdata.buffer();
557  for (int k = 0; k < gnk; k++) {
558  int jkoff = k * nij + koff;
559  for (int j = 0; j < gnj; j++) {
560  int ijkoff = j * ni + jkoff;
561  for (int i = 0; i < gni; i++, index++) {
562  buf[i + ijkoff] += gbuf[index];
563  }
564  }
565  }
566  return(*this);
567  }
568 
569  // extract a subgrid from this grid
570  // subgrid must fit within this grid's index range
572  ASSERT(IndexRange(g) <= IndexRange(*this));
573  int gni = g.nextent.i;
574  int gnj = g.nextent.j;
575  int gnk = g.nextent.k;
576  int index = 0;
577  int ni = nextent.i;
578  int nij = nextent.i * nextent.j;
579  int koff = (g.nlower.k - nlower.k) * nij
580  + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
581  T *gbuf = g.gdata.buffer();
582  const T *buf = gdata.buffer();
583  for (int k = 0; k < gnk; k++) {
584  int jkoff = k * nij + koff;
585  for (int j = 0; j < gnj; j++) {
586  int ijkoff = j * ni + jkoff;
587  for (int i = 0; i < gni; i++, index++) {
588  gbuf[index] = buf[i + ijkoff];
589  }
590  }
591  }
592  }
593 
594  private:
595  T gdata[N];
596  };
597 
598  // storage and indexing for 3D lattice of grid points
599  template <class T>
600  class Grid : public IndexRange {
601  public:
602  Grid() { }
603  void init(const IndexRange& n) {
604  nlower = n.lower();
605  nextent = n.extent();
606  gdata.resize(nn());
607  }
608  void set(int pia, int pni, int pja, int pnj, int pka, int pnk) {
609  IndexRange::set(pia, pni, pja, pnj, pka, pnk);
610  gdata.resize(nn());
611  }
612  void setbounds(int pia, int pib, int pja, int pjb, int pka, int pkb) {
613  IndexRange::setbounds(pia, pib, pja, pjb, pka, pkb);
614  gdata.resize(nn());
615  }
616  void resize(int n) { // reserve space but don't set grid indexing
617  gdata.resize(n);
618  }
619  const T& operator()(int i, int j, int k) const {
620 #ifdef DEBUG_MSM
621  return elem(i,j,k);
622 #else
623  return gdata[flatindex(i,j,k)];
624 #endif
625  }
626  const T& operator()(const Ivec& n) const {
627  return this->operator()(n.i, n.j, n.k);
628  }
629  const T& elem(int i, int j, int k) const {
630  if (i<ia() || i>ib() || j<ja() || j>jb() || k<ka() || k>kb()) {
631  char msg[200];
632  snprintf(msg, sizeof(msg), "Grid indexing:\n"
633  "ia=%d, ib=%d, i=%d\n"
634  "ja=%d, jb=%d, j=%d\n"
635  "ka=%d, kb=%d, k=%d\n",
636  ia(), ib(), i, ja(), jb(), j, ka(), kb(), k);
637  NAMD_die(msg);
638  }
639  return gdata[flatindex(i,j,k)];
640  }
641  T& operator()(int i, int j, int k) {
642 #ifdef DEBUG_MSM
643  return elem(i,j,k);
644 #else
645  return gdata[flatindex(i,j,k)];
646 #endif
647  }
648  T& operator()(const Ivec& n) {
649  return this->operator()(n.i, n.j, n.k);
650  }
651  T& elem(int i, int j, int k) {
652  if (i<ia() || i>ib() || j<ja() || j>jb() || k<ka() || k>kb()) {
653  char msg[200];
654  snprintf(msg, sizeof(msg), "Grid indexing:\n"
655  "ia=%d, ib=%d, i=%d\n"
656  "ja=%d, jb=%d, j=%d\n"
657  "ka=%d, kb=%d, k=%d\n",
658  ia(), ib(), i, ja(), jb(), j, ka(), kb(), k);
659  NAMD_die(msg);
660  }
661  return gdata[flatindex(i,j,k)];
662  }
663  int flatindex(int i, int j, int k) const {
664  return ((k-ka())*nj() + (j-ja()))*ni() + (i-ia());
665  }
666  const Array<T>& data() const { return gdata; }
667  Array<T>& data() { return gdata; }
668 
669  // use to zero out grid
670  void reset(const T& t) {
671  T *buf = gdata.buffer();
672  int len = nn();
673  for (int n = 0; n < len; n++) { buf[n] = t; }
674  }
675 
676  // use to modify the indexing by changing lower corner
677  void updateLower(const Ivec& n) { nlower = n; }
678 
679  // accumulate another grid into this grid
680  // the grid to be added must fit within this grid's index range
682 #if 1
683  ASSERT(IndexRange(g) <= IndexRange(*this));
684 #else
685  if ( ! (IndexRange(g) <= IndexRange(*this)) ) {
686  Grid<T> tmp = *this;
687  // expand myself to hold sum
688  int ia = nlower.i;
689  if (ia > g.nlower.i) ia = g.nlower.i;
690  int ja = nlower.j;
691  if (ja > g.nlower.j) ja = g.nlower.j;
692  int ka = nlower.k;
693  if (ka > g.nlower.k) ka = g.nlower.k;
694  int ib1 = nlower.i + nextent.i;
695  int gib1 = g.nlower.i + g.nextent.i;
696  if (ib1 < gib1) ib1 = gib1;
697  int jb1 = nlower.j + nextent.j;
698  int gjb1 = g.nlower.j + g.nextent.j;
699  if (jb1 < gjb1) jb1 = gjb1;
700  int kb1 = nlower.k + nextent.k;
701  int gkb1 = g.nlower.k + g.nextent.k;
702  if (kb1 < gkb1) kb1 = gkb1;
703  setbounds(ia, ib1-1, ja, jb1-1, ka, kb1-1);
704  reset(0); // make sure constructor for T accepts "0" as its zero
705  // now copy "tmp" grid elements into my expanded self
706  int index = 0;
707  int ni = nextent.i;
708  int nij = nextent.i * nextent.j;
709  int koff = (tmp.nlower.k - nlower.k) * nij
710  + (tmp.nlower.j - nlower.j) * ni + (tmp.nlower.i - nlower.i);
711  const T *gbuf = tmp.gdata.buffer();
712  T *buf = gdata.buffer();
713  for (int k = 0; k < tmp.nextent.k; k++) {
714  int jkoff = k * nij + koff;
715  for (int j = 0; j < tmp.nextent.j; j++) {
716  int ijkoff = j * ni + jkoff;
717  for (int i = 0; i < tmp.nextent.i; i++, index++) {
718  buf[i + ijkoff] = gbuf[index];
719  }
720  }
721  }
722  }
723 #endif
724  int gni = g.nextent.i;
725  int gnj = g.nextent.j;
726  int gnk = g.nextent.k;
727  int index = 0;
728  int ni = nextent.i;
729  int nij = nextent.i * nextent.j;
730  int koff = (g.nlower.k - nlower.k) * nij
731  + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
732  const T *gbuf = g.gdata.buffer();
733  T *buf = gdata.buffer();
734  for (int k = 0; k < gnk; k++) {
735  int jkoff = k * nij + koff;
736  for (int j = 0; j < gnj; j++) {
737  int ijkoff = j * ni + jkoff;
738  for (int i = 0; i < gni; i++, index++) {
739  buf[i + ijkoff] += gbuf[index];
740  }
741  }
742  }
743  return(*this);
744  }
745 
746  // extract a subgrid from this grid
747  // subgrid must fit within this grid's index range
748  void extract(Grid<T>& g) {
749  ASSERT(IndexRange(g) <= IndexRange(*this));
750  int gni = g.nextent.i;
751  int gnj = g.nextent.j;
752  int gnk = g.nextent.k;
753  int index = 0;
754  int ni = nextent.i;
755  int nij = nextent.i * nextent.j;
756  int koff = (g.nlower.k - nlower.k) * nij
757  + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
758  T *gbuf = g.gdata.buffer();
759  const T *buf = gdata.buffer();
760  for (int k = 0; k < gnk; k++) {
761  int jkoff = k * nij + koff;
762  for (int j = 0; j < gnj; j++) {
763  int ijkoff = j * ni + jkoff;
764  for (int i = 0; i < gni; i++, index++) {
765  gbuf[index] = buf[i + ijkoff];
766  }
767  }
768  }
769  }
770 
771  // accumulate a fixed size grid into this grid
772  // the grid to be added must fit within this grid's index range
773  template <int N>
775  ASSERT(IndexRange(g) <= IndexRange(*this));
776  int gni = g.nextent.i;
777  int gnj = g.nextent.j;
778  int gnk = g.nextent.k;
779  int index = 0;
780  int ni = nextent.i;
781  int nij = nextent.i * nextent.j;
782  int koff = (g.nlower.k - nlower.k) * nij
783  + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
784  const T *gbuf = g.buffer();
785  T *buf = gdata.buffer();
786  for (int k = 0; k < gnk; k++) {
787  int jkoff = k * nij + koff;
788  for (int j = 0; j < gnj; j++) {
789  int ijkoff = j * ni + jkoff;
790  for (int i = 0; i < gni; i++, index++) {
791  buf[i + ijkoff] += gbuf[index];
792  }
793  }
794  }
795  return(*this);
796  }
797 
798  // extract a subgrid from this grid
799  // subgrid must fit within this grid's index range
800  template <int N>
802  ASSERT(IndexRange(g) <= IndexRange(*this));
803  int gni = g.nextent.i;
804  int gnj = g.nextent.j;
805  int gnk = g.nextent.k;
806  int index = 0;
807  int ni = nextent.i;
808  int nij = nextent.i * nextent.j;
809  int koff = (g.nlower.k - nlower.k) * nij
810  + (g.nlower.j - nlower.j) * ni + (g.nlower.i - nlower.i);
811  T *gbuf = g.buffer();
812  const T *buf = gdata.buffer();
813  for (int k = 0; k < gnk; k++) {
814  int jkoff = k * nij + koff;
815  for (int j = 0; j < gnj; j++) {
816  int ijkoff = j * ni + jkoff;
817  for (int i = 0; i < gni; i++, index++) {
818  gbuf[index] = buf[i + ijkoff];
819  }
820  }
821  }
822  }
823 
824  private:
825  Array<T> gdata;
826  }; // Grid
827 
828 
830  //
831  // Map object
832  //
834 
835  // index a block from the MSM grid hierarchy
836  struct BlockIndex {
837  int level;
839  BlockIndex() : level(0), n(0) { }
840  BlockIndex(int ll, const Ivec& nn) : level(ll), n(nn) { }
841 #ifdef MSM_MIGRATION
842  void pup(PUP::er& p) {
843  p|level, p|n;
844  }
845 #endif
846  };
847 
848  // for uppermost levels of hierarchy
849  // fold out image charges along periodic boundaries
850  // to fill up desired block size
851  struct FoldFactor {
852  int active; // is some numrep dimension > 1?
853  Ivec numrep; // number of replications along each dimension
854  FoldFactor() : active(0), numrep(1) { }
855  FoldFactor(int i, int j, int k) { set(i,j,k); }
856  void set(int i, int j, int k) {
857  if (i <= 0) i = 1;
858  if (j <= 0) j = 1;
859  if (k <= 0) k = 1;
860  if (i > 1 || j > 1 || k > 1) active = 1;
861  numrep = Ivec(i, j, k);
862  }
863  };
864 
865  // sending part of an extended grid calculation to another block
866  struct BlockSend {
867  BlockIndex nblock; // relative block index
868  IndexRange nrange; // relative grid index range
869  BlockIndex nblock_wrap; // true block index
870  IndexRange nrange_wrap; // true grid index range
871  void reset() {
872  nblock = BlockIndex();
873  nrange = IndexRange();
876  } // reset
877 #ifdef MSM_MIGRATION
878  void pup(PUP::er& p) {
880  }
881 #endif
882  };
883 
884  struct PatchSend {
885  IndexRange nrange; // true grid index range from my block
886  IndexRange nrange_unwrap; // relative grid index range for patch
887  int patchID; // patch ID
888  void reset() {
889  nrange = IndexRange();
891  patchID = -1;
892  } // reset
893  };
894 
895  // one PatchDiagram for each patch
896  // maintain a Grid of PatchDiagram, indexed by patch ID
897  struct PatchDiagram {
898  IndexRange nrange; // shows subset of MSM h-grid covering this patch
899  Array<BlockSend> send; // array of blocks to which this patch sends
900  int numRecvs; // number of blocks from which this patch receives
901  void reset() {
902  nrange = IndexRange();
903  send.resize(0);
904  numRecvs = 0;
905  } // reset
906  };
907 
908  // one BlockDiagram for each block of each level of each MSM grid
909  // maintain a Grid of BlockDiagram for each level
910  struct BlockDiagram {
911  IndexRange nrange; // subset of MSM grid for this block
912  IndexRange nrangeCutoff; // expanded subgrid for cutoff calculation
913  IndexRange nrangeRestricted; // (level+1) subgrid for restriction
914  IndexRange nrangeProlongated; // (level-1) subgrid for prolongation
915  Array<BlockSend> sendUp; // send up charge to blocks on (level+1)
916  Array<BlockSend> sendAcross; // send across potential to blocks on (level)
917  Array<int> indexGridCutoff; // index of MsmGridCutoff chares to calculate
918  // each charge -> potential block interaction
919  Array<int> recvGridCutoff; // index of MsmGridCutoff chares contributing
920  // back into my potential block
921  Array<BlockSend> sendDown; // send down potential to blocks on (level-1)
922  Array<PatchSend> sendPatch; // send my (level=0) potential block to patch
923  int numRecvsCharge; // number of expected receives of charge
924  int numRecvsPotential; // number of expected receives of potential
925 
926  void reset() {
927  nrange = IndexRange();
931  sendUp.resize(0);
932  sendAcross.resize(0);
933  sendDown.resize(0);
934  sendPatch.resize(0);
935  numRecvsCharge = 0;
936  numRecvsPotential = 0;
937  } // reset
938  };
939 
940 
941  struct Map {
942  Array<IndexRange> gridrange; // dimensions for each MSM grid level
943 
944  Array<Grid<Float> > gc; // grid constant weights for each level
945  Array<Grid<Float> > gvc; // virial grid weights for each level
946  Grid<Float> grespro; // restriction / prolongation nonzero stencil
947  // requires correct IndexOffset array
948 
949  Array<Grid<C1Matrix> > gc_c1hermite; // grid constant weights C1 Hermite
950  Array<Grid<C1Matrix> > gvc_c1hermite; // virial grid weights C1 Hermite
951  Array<Grid<C1Matrix> > gres_c1hermite; // restriction stencil C1 Hermite
952  Array<Grid<C1Matrix> > gpro_c1hermite; // prolongation stencil C1 Hermite
953  // requires index offsets
954 
957 
958  int ispx, ispy, ispz; // is periodic in x, y, z?
959 
960  Array<int> bsx, bsy, bsz; // block size in x, y, z for each level
961 
962  Array<FoldFactor> foldfactor; // for uppermost grid levels
963  // replicate periodic dimensions in order to fill up block size
964 
965  // clip index to grid level, using periodicity flags
966  Ivec clipIndexToLevel(const Ivec& n, int level) const {
967  ASSERT(level >= 0 && level < gridrange.len());
968  Ivec pn(n);
969  if ( ! ispx) {
970  int a = gridrange[level].ia();
971  int b = gridrange[level].ib();
972  if (pn.i < a) pn.i = a;
973  if (pn.i > b) pn.i = b;
974  }
975  if ( ! ispy) {
976  int a = gridrange[level].ja();
977  int b = gridrange[level].jb();
978  if (pn.j < a) pn.j = a;
979  if (pn.j > b) pn.j = b;
980  }
981  if ( ! ispz) {
982  int a = gridrange[level].ka();
983  int b = gridrange[level].kb();
984  if (pn.k < a) pn.k = a;
985  if (pn.k > b) pn.k = b;
986  }
987  return pn;
988  }
989 
990  // determine relative (unwrapped) block index for the given grid index
991  BlockIndex blockOfGridIndex(const Ivec& n, int level) const {
992  ASSERT(level >= 0 && level < gridrange.len());
993  BlockIndex bn;
994  // we want floor((i - ia) / bsx), etc.
995  // modify case i < ia to avoid integer division of negative numbers
996  int d = n.i - gridrange[level].ia();
997  bn.n.i = (d >= 0 ? d / bsx[level] : -((-d+bsx[level]-1) / bsx[level]));
998  d = n.j - gridrange[level].ja();
999  bn.n.j = (d >= 0 ? d / bsy[level] : -((-d+bsy[level]-1) / bsy[level]));
1000  d = n.k - gridrange[level].ka();
1001  bn.n.k = (d >= 0 ? d / bsz[level] : -((-d+bsz[level]-1) / bsz[level]));
1002  bn.level = level;
1003  return bn;
1004  }
1005 
1006  // determine relative (unwrapped) block index for the given grid index
1007  // for unfolded replication of image charges
1008  BlockIndex blockOfGridIndexFold(const Ivec& n, int level) const {
1009  ASSERT(level >= 0 && level < gridrange.len());
1010  BlockIndex bn;
1011  int bsi = foldfactor[level].numrep.i * bsx[level];
1012  int bsj = foldfactor[level].numrep.j * bsy[level];
1013  int bsk = foldfactor[level].numrep.k * bsz[level];
1014  // we want floor((i - ia) / bsx), etc.
1015  // modify case i < ia to avoid integer division of negative numbers
1016  int d = n.i - gridrange[level].ia();
1017  bn.n.i = (d >= 0 ? d / bsi : -((-d+bsi-1) / bsi));
1018  d = n.j - gridrange[level].ja();
1019  bn.n.j = (d >= 0 ? d / bsj : -((-d+bsj-1) / bsj));
1020  d = n.k - gridrange[level].ka();
1021  bn.n.k = (d >= 0 ? d / bsk : -((-d+bsk-1) / bsk));
1022  bn.level = level;
1023  return bn;
1024  }
1025 
1026  // find the natural index range of the given relative block number
1028  ASSERT(nb.level >= 0 && nb.level < gridrange.len());
1029  IndexRange nr;
1030  int ia = nb.n.i * bsx[nb.level] + gridrange[nb.level].ia();
1031  int ja = nb.n.j * bsy[nb.level] + gridrange[nb.level].ja();
1032  int ka = nb.n.k * bsz[nb.level] + gridrange[nb.level].ka();
1033  nr.set(ia, bsx[nb.level], ja, bsy[nb.level], ka, bsz[nb.level]);
1034  return nr;
1035  }
1036 
1037  // find the natural index range of the given relative block number
1038  // for unfolded replication of image charges
1040  ASSERT(nb.level >= 0 && nb.level < gridrange.len());
1041  int bsi = foldfactor[nb.level].numrep.i * bsx[nb.level];
1042  int bsj = foldfactor[nb.level].numrep.j * bsy[nb.level];
1043  int bsk = foldfactor[nb.level].numrep.k * bsz[nb.level];
1044  IndexRange nr;
1045  int ia = nb.n.i * bsi + gridrange[nb.level].ia();
1046  int ja = nb.n.j * bsj + gridrange[nb.level].ja();
1047  int ka = nb.n.k * bsk + gridrange[nb.level].ka();
1048  nr.set(ia, bsi, ja, bsj, ka, bsk);
1049  return nr;
1050  }
1051 
1052  // clip the natural block index range to not exceed the given index range
1054  const IndexRange& nrange) const {
1055  IndexRange nr = indexRangeOfBlock(nb);
1056  int nia = nrange.ia();
1057  int nib = nrange.ib();
1058  int nja = nrange.ja();
1059  int njb = nrange.jb();
1060  int nka = nrange.ka();
1061  int nkb = nrange.kb();
1062  int ia = nr.ia();
1063  if (ia < nia) ia = nia;
1064  int ib = nr.ib();
1065  if (ib > nib) ib = nib;
1066  int ja = nr.ja();
1067  if (ja < nja) ja = nja;
1068  int jb = nr.jb();
1069  if (jb > njb) jb = njb;
1070  int ka = nr.ka();
1071  if (ka < nka) ka = nka;
1072  int kb = nr.kb();
1073  if (kb > nkb) kb = nkb;
1074  nr.setbounds(ia, ib, ja, jb, ka, kb);
1075  return nr;
1076  }
1077 
1078  // clip the natural block index range to not exceed the given index range
1079  // for unfolded replication of image charges
1081  const IndexRange& nrange) const {
1083  int nia = nrange.ia();
1084  int nib = nrange.ib();
1085  int nja = nrange.ja();
1086  int njb = nrange.jb();
1087  int nka = nrange.ka();
1088  int nkb = nrange.kb();
1089  int ia = nr.ia();
1090  if (ia < nia) ia = nia;
1091  int ib = nr.ib();
1092  if (ib > nib) ib = nib;
1093  int ja = nr.ja();
1094  if (ja < nja) ja = nja;
1095  int jb = nr.jb();
1096  if (jb > njb) jb = njb;
1097  int ka = nr.ka();
1098  if (ka < nka) ka = nka;
1099  int kb = nr.kb();
1100  if (kb > nkb) kb = nkb;
1101  nr.setbounds(ia, ib, ja, jb, ka, kb);
1102  return nr;
1103  }
1104 
1105  // set the nblock_wrap and nrange_wrap fields based on periodicity
1106  void wrapBlockSend(BlockSend& bs) const {
1107  BlockIndex nb = bs.nblock;
1108  IndexRange nr = bs.nrange;
1109  int level = bs.nblock.level;
1110  ASSERT(level >= 0 && level < blockLevel.len());
1111  int ni = blockLevel[level].ni();
1112  int nj = blockLevel[level].nj();
1113  int nk = blockLevel[level].nk();
1114  int di=0, dj=0, dk=0;
1115  if (ispx) {
1116  while (nb.n.i < 0) {
1117  nb.n.i += ni;
1118  di += ni * bsx[level];
1119  }
1120  while (nb.n.i >= ni) {
1121  nb.n.i -= ni;
1122  di -= ni * bsx[level];
1123  }
1124  }
1125  if (ispy) {
1126  while (nb.n.j < 0) {
1127  nb.n.j += nj;
1128  dj += nj * bsy[level];
1129  }
1130  while (nb.n.j >= nj) {
1131  nb.n.j -= nj;
1132  dj -= nj * bsy[level];
1133  }
1134  }
1135  if (ispz) {
1136  while (nb.n.k < 0) {
1137  nb.n.k += nk;
1138  dk += nk * bsz[level];
1139  }
1140  while (nb.n.k >= nk) {
1141  nb.n.k -= nk;
1142  dk -= nk * bsz[level];
1143  }
1144  }
1145  int ia = nr.ia();
1146  int ib = nr.ib();
1147  int ja = nr.ja();
1148  int jb = nr.jb();
1149  int ka = nr.ka();
1150  int kb = nr.kb();
1151  nr.setbounds(ia + di, ib + di, ja + dj, jb + dj, ka + dk, kb + dk);
1152  bs.nblock_wrap = nb;
1153  bs.nrange_wrap = nr;
1154  }
1155 
1156  // set the nblock_wrap and nrange_wrap fields based on periodicity
1157  // for unfolded replication of image charges
1158  void wrapBlockSendFold(BlockSend& bs) const {
1159  BlockIndex nb = bs.nblock;
1160  IndexRange nr = bs.nrange;
1161  int level = bs.nblock.level;
1162  ASSERT(level >= 0 && level < blockLevel.len());
1163  int foldi = foldfactor[level].numrep.i;
1164  int foldj = foldfactor[level].numrep.j;
1165  int foldk = foldfactor[level].numrep.k;
1166  int ni = blockLevel[level].ni();
1167  int nj = blockLevel[level].nj();
1168  int nk = blockLevel[level].nk();
1169  int bsi = foldi * bsx[level];
1170  int bsj = foldj * bsy[level];
1171  int bsk = foldk * bsz[level];
1172  int di=0, dj=0, dk=0;
1173  if (ispx) {
1174  while (nb.n.i < 0) {
1175  nb.n.i += ni;
1176  di += ni * bsi;
1177  }
1178  while (nb.n.i >= ni) {
1179  nb.n.i -= ni;
1180  di -= ni * bsi;
1181  }
1182  }
1183  if (ispy) {
1184  while (nb.n.j < 0) {
1185  nb.n.j += nj;
1186  dj += nj * bsj;
1187  }
1188  while (nb.n.j >= nj) {
1189  nb.n.j -= nj;
1190  dj -= nj * bsj;
1191  }
1192  }
1193  if (ispz) {
1194  while (nb.n.k < 0) {
1195  nb.n.k += nk;
1196  dk += nk * bsk;
1197  }
1198  while (nb.n.k >= nk) {
1199  nb.n.k -= nk;
1200  dk -= nk * bsk;
1201  }
1202  }
1203  int ia = nr.ia();
1204  int ib = nr.ib();
1205  int ja = nr.ja();
1206  int jb = nr.jb();
1207  int ka = nr.ka();
1208  int kb = nr.kb();
1209  nr.setbounds(ia + di, ib + di, ja + dj, jb + dj, ka + dk, kb + dk);
1210  bs.nblock_wrap = nb;
1211  bs.nrange_wrap = nr;
1212  }
1213 
1214  void wrapBlockIndex(BlockIndex& bn) const {
1215  int level = bn.level;
1216  int ni = blockLevel[level].ni();
1217  int nj = blockLevel[level].nj();
1218  int nk = blockLevel[level].nk();
1219  if (ispx) {
1220  while (bn.n.i < 0) {
1221  bn.n.i += ni;
1222  }
1223  while (bn.n.i >= ni) {
1224  bn.n.i -= ni;
1225  }
1226  }
1227  if (ispy) {
1228  while (bn.n.j < 0) {
1229  bn.n.j += nj;
1230  }
1231  while (bn.n.j >= nj) {
1232  bn.n.j -= nj;
1233  }
1234  }
1235  if (ispz) {
1236  while (bn.n.k < 0) {
1237  bn.n.k += nk;
1238  }
1239  while (bn.n.k >= nk) {
1240  bn.n.k -= nk;
1241  }
1242  }
1243  }
1244 
1245  }; // Map
1246 
1247 
1248  struct AtomCoord {
1251  int id;
1252  };
1253 
1256 
1257  struct PatchData;
1259 
1260 } // namespace msm
1261 
1262 #endif // MSMMAP_H
int amax
Definition: MsmMap.h:276
Array< int > bsz
Definition: MsmMap.h:960
IndexRange clipBlockToIndexRangeFold(const BlockIndex &nb, const IndexRange &nrange) const
Definition: MsmMap.h:1080
GridFixed< T, N > & operator+=(const GridFixed< T, N > &g)
Definition: MsmMap.h:545
Array< Grid< C1Matrix > > gpro_c1hermite
Definition: MsmMap.h:952
Array< PatchDiagram > patchList
Definition: MsmMap.h:955
void reset(const T &t)
Definition: MsmMap.h:670
void resize(int n)
Definition: MsmMap.h:254
void init(int nelems)
Definition: MsmMap.h:342
Array< int > recvGridCutoff
Definition: MsmMap.h:919
int len() const
Definition: MsmMap.h:218
Definition: MsmMap.h:187
void reset()
Definition: MsmMap.h:926
BlockIndex nblock_wrap
Definition: MsmMap.h:869
Array(int n)
Definition: MsmMap.h:211
#define ASSERT(expr)
Definition: MsmMap.h:68
IndexRange nrange_unwrap
Definition: MsmMap.h:886
Definition: MsmMap.h:187
int k
Definition: MsmMap.h:411
void reset()
Definition: MsmMap.h:888
const Array< T > & data() const
Definition: MsmMap.h:666
Grid()
Definition: MsmMap.h:602
Grid< T > & operator+=(const GridFixed< T, N > &g)
Definition: MsmMap.h:774
friend C1Vector operator+(const C1Vector &u, const C1Vector &v)
Definition: MsmMap.h:101
Array< int > bsy
Definition: MsmMap.h:960
int ispx
Definition: MsmMap.h:958
int max() const
Definition: MsmMap.h:219
int ka() const
Definition: MsmMap.h:438
void wrapBlockSendFold(BlockSend &bs) const
Definition: MsmMap.h:1158
Array & operator=(const Array &a)
Definition: MsmMap.h:214
Definition: Vector.h:64
BlockIndex nblock
Definition: MsmMap.h:867
IndexRange indexRangeOfBlockFold(const BlockIndex &nb) const
Definition: MsmMap.h:1039
double Double
Definition: MsmMap.h:75
T & operator()(int i, int j, int k)
Definition: MsmMap.h:641
int i
Definition: MsmMap.h:411
int operator==(const Ivec &n)
Definition: MsmMap.h:414
const T & operator[](int i) const
Definition: MsmMap.h:220
int ispy
Definition: MsmMap.h:958
float Real
Definition: common.h:107
const T * buffer(int &n) const
Definition: MsmMap.h:261
int ia() const
Definition: MsmMap.h:434
BlockIndex blockOfGridIndexFold(const Ivec &n, int level) const
Definition: MsmMap.h:1008
Definition: MsmMap.h:187
BlockIndex blockOfGridIndex(const Ivec &n, int level) const
Definition: MsmMap.h:991
int ispz
Definition: MsmMap.h:958
Array< BlockSend > sendUp
Definition: MsmMap.h:915
void set(int pia, int pni, int pja, int pnj, int pka, int pnk)
Definition: MsmMap.h:426
void extract(Grid< T > &g)
Definition: MsmMap.h:748
void set(Float r)
Definition: MsmMap.h:89
void setbounds(int pia, int pib, int pja, int pjb, int pka, int pkb)
Definition: MsmMap.h:431
void swap(Array< T > &s, Array< T > &t)
Definition: MsmMap.h:319
T * buffer(int &n)
Definition: MsmMap.h:262
void setbounds(int pia, int pib, int pja, int pjb, int pka, int pkb)
Definition: MsmMap.h:480
T * abuffer
Definition: MsmMap.h:275
Ivec clipIndexToLevel(const Ivec &n, int level) const
Definition: MsmMap.h:966
const T & operator()(const Ivec &n) const
Definition: MsmMap.h:626
void set(int pia, int pni, int pja, int pnj, int pka, int pnk)
Definition: MsmMap.h:608
Array< Force > ForceArray
Definition: MsmMap.h:1255
T & operator()(int i, int j, int k)
Definition: MsmMap.h:506
Array< BlockSend > send
Definition: MsmMap.h:899
void setbounds(int pia, int pib, int pja, int pjb, int pka, int pkb)
Definition: MsmMap.h:612
const T & elem(int i, int j, int k) const
Definition: MsmMap.h:629
Array< Grid< BlockDiagram > > blockLevel
Definition: MsmMap.h:956
int flatindex(int i, int j, int k) const
Definition: MsmMap.h:528
T & operator()(const Ivec &n)
Definition: MsmMap.h:648
Ivec extent() const
Definition: MsmMap.h:445
int j
Definition: MsmMap.h:411
Definition: MsmMap.h:187
friend C1Vector operator*(const C1Matrix &m, const C1Vector &u)
Definition: MsmMap.h:116
void init(const IndexRange &n)
Definition: MsmMap.h:603
Array< BlockSend > sendAcross
Definition: MsmMap.h:916
void setmax(int m)
Definition: MsmMap.h:281
~Array()
Definition: MsmMap.h:213
Ivec nextent
Definition: MsmMap.h:459
Array< Grid< C1Matrix > > gres_c1hermite
Definition: MsmMap.h:951
Array< AtomCoord > AtomCoordArray
Definition: MsmMap.h:1254
C1Vector(Float r=0)
Definition: MsmMap.h:88
const T * buffer() const
Definition: MsmMap.h:531
PriorityQueue(int nelems=0)
Definition: MsmMap.h:339
int kb() const
Definition: MsmMap.h:439
Array< T > & data()
Definition: MsmMap.h:667
int nn() const
Definition: MsmMap.h:443
Array< PatchData * > PatchPtrArray
Definition: MsmMap.h:1257
void updateLower(const Ivec &n)
Definition: MsmMap.h:541
IndexRange nrangeRestricted
Definition: MsmMap.h:913
T & elem(int i)
Definition: MsmMap.h:242
int jb() const
Definition: MsmMap.h:437
void set(int i, int j, int k)
Definition: MsmMap.h:856
IndexRange nrange
Definition: MsmMap.h:868
int nj() const
Definition: MsmMap.h:441
void extract(GridFixed< T, N > &g)
Definition: MsmMap.h:801
BlockIndex(int ll, const Ivec &nn)
Definition: MsmMap.h:840
IndexRange nrange
Definition: MsmMap.h:911
Array< int > bsx
Definition: MsmMap.h:960
int alen
Definition: MsmMap.h:276
IndexRange nrangeCutoff
Definition: MsmMap.h:912
Position position
Definition: MsmMap.h:1249
void extract(GridFixed< T, N > &g)
Definition: MsmMap.h:571
Array< Grid< C1Matrix > > gvc_c1hermite
Definition: MsmMap.h:950
Array< IndexRange > gridrange
Definition: MsmMap.h:942
Array< PatchSend > sendPatch
Definition: MsmMap.h:922
void NAMD_die(const char *err_msg)
Definition: common.C:83
FoldFactor(int i, int j, int k)
Definition: MsmMap.h:855
int flatindex(int i, int j, int k) const
Definition: MsmMap.h:663
void wrapBlockIndex(BlockIndex &bn) const
Definition: MsmMap.h:1214
Definition: MsmMap.h:187
IndexRange clipBlockToIndexRange(const BlockIndex &nb, const IndexRange &nrange) const
Definition: MsmMap.h:1053
const T & elem(int i) const
Definition: MsmMap.h:227
float Float
Definition: MsmMap.h:74
void updateLower(const Ivec &n)
Definition: MsmMap.h:677
T * buffer()
Definition: MsmMap.h:532
int numRecvsPotential
Definition: MsmMap.h:924
const T & operator()(int i, int j, int k) const
Definition: MsmMap.h:619
Float melem[C1_MATRIX_SIZE]
Definition: MsmMap.h:111
Definition: MsmMap.h:187
void copy(const Array &a)
Definition: MsmMap.h:309
Definition: Array.h:10
const T & elem(int i, int j, int k) const
Definition: MsmMap.h:494
void wrapBlockSend(BlockSend &bs) const
Definition: MsmMap.h:1106
int nk() const
Definition: MsmMap.h:442
T & operator()(const Ivec &n)
Definition: MsmMap.h:513
IndexRange nrange_wrap
Definition: MsmMap.h:870
Array< FoldFactor > foldfactor
Definition: MsmMap.h:962
const T & operator()(const Ivec &n) const
Definition: MsmMap.h:491
Definition: MsmMap.h:187
void resize(int n)
Definition: MsmMap.h:616
void insert(const T &t)
Definition: MsmMap.h:349
Grid< T > & operator+=(const Grid< T > &g)
Definition: MsmMap.h:681
Array< int > indexGridCutoff
Definition: MsmMap.h:917
T & operator[](int i)
Definition: MsmMap.h:235
IndexRange nrange
Definition: MsmMap.h:898
void set(int pia, int pni, int pja, int pnj, int pka, int pnk)
Definition: MsmMap.h:476
Array< Grid< Float > > gc
Definition: MsmMap.h:944
int ib() const
Definition: MsmMap.h:435
const T & operator()(int i, int j, int k) const
Definition: MsmMap.h:484
Ivec lower() const
Definition: MsmMap.h:444
Array< Grid< Float > > gvc
Definition: MsmMap.h:945
void reset()
Definition: MsmMap.h:871
Definition: MsmMap.h:187
void set(Float r)
Definition: MsmMap.h:113
int ja() const
Definition: MsmMap.h:436
Grid< Float > grespro
Definition: MsmMap.h:946
T * buffer()
Definition: MsmMap.h:260
IndexRange indexRangeOfBlock(const BlockIndex &nb) const
Definition: MsmMap.h:1027
friend Float operator*(const C1Vector &u, const C1Vector &v)
Definition: MsmMap.h:96
int ni() const
Definition: MsmMap.h:440
void init(const IndexRange &n)
Definition: MsmMap.h:471
IndexRange nrange
Definition: MsmMap.h:885
Array< Grid< C1Matrix > > gc_c1hermite
Definition: MsmMap.h:949
int operator<=(const IndexRange &n)
Definition: MsmMap.h:446
C1Vector & operator+=(const C1Vector &v)
Definition: MsmMap.h:92
void reset()
Definition: MsmMap.h:901
Ivec(int n=0)
Definition: MsmMap.h:412
IndexRange nrangeProlongated
Definition: MsmMap.h:914
Ivec(int ni, int nj, int nk)
Definition: MsmMap.h:413
Float velem[C1_VECTOR_SIZE]
Definition: MsmMap.h:87
void reset(const T &t)
Definition: MsmMap.h:535
void reset(const T &t)
Definition: MsmMap.h:263
const T * buffer() const
Definition: MsmMap.h:259
Array< BlockSend > sendDown
Definition: MsmMap.h:921
int numRecvsCharge
Definition: MsmMap.h:923
T & elem(int i, int j, int k)
Definition: MsmMap.h:516
Array(const Array &a)
Definition: MsmMap.h:212
C1Matrix(Float r=0)
Definition: MsmMap.h:112
void append(const T &t)
Definition: MsmMap.h:250
T & elem(int i, int j, int k)
Definition: MsmMap.h:651