MPQC  3.0.0-alpha
mtensor.h
1 //
2 // mtensor.h
3 //
4 // Copyright (C) 2009 Edward Valeev
5 //
6 // Author: Edward Valeev <evaleev@vt.edu>
7 // Maintainer: EV
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #ifndef _ccr12_mtensor_h
29 #define _ccr12_mtensor_h
30 
31 #include <numeric>
32 #include <chemistry/qc/ccr12/tensor.h>
33 #include <chemistry/qc/ccr12/ccr12_info.h>
34 #include <math/scmat/matrix.h>
35 
36 namespace {
37  void print_tile(double* data, int n01, int n23) {
38  using namespace sc;
40  new SCDimension(n23));
41  mat.assign(data);
42  mat.print("tile");
43  }
44 
45 }
46 
47 namespace sc {
48 
50  template <typename I> inline std::pair<I,I> intersect(const std::pair<I,I>& range1, const std::pair<I,I>& range2) {
51  I start_max = std::max(range1.first, range2.first);
52  I fence_min = std::min(range1.second, range2.second);
53  if (start_max < fence_min) return make_pair(start_max, fence_min);
54  return make_pair(I(0),I(0));
55  }
56 
58  template <typename I> inline bool in(const std::pair<I,I>& r, const std::pair<I,I>& range) {
59  return r.first >= range.first && r.second <= range.second;
60  }
62  template <typename I> inline bool in(I i, const std::pair<I,I>& range) {
63  return i >= range.first && i < range.second;
64  }
65 
67 
69  template <size_t NDIM>
70  class MTensor {
71  public:
72  typedef CCR12_Info Info;
73  typedef long tile_index;
74  typedef long element_index;
75  typedef std::pair<tile_index, tile_index> tile_range; //< [start, fence)
76  typedef std::vector<tile_range> tile_ranges;
77  typedef std::pair<element_index, element_index> element_range; //< [start, fence)
78  typedef std::vector<element_range> element_ranges;
79  typedef std::vector<element_index> element_index_map;
80  static const size_t ndim = NDIM;
81 
85  MTensor(Info const* info,
86  Tensor* tensor,
87  const tile_ranges& range) :
88  tensor_(tensor), info_(info), range_(range)
89  {
90  MPQC_ASSERT(range.size() == NDIM);
91  }
92 
104  void convert(const Ref<DistArray4>& src, unsigned int tbtype,
105  const element_index_map& eimap0,
106  const element_index_map& eimap1,
107  const element_index_map& eimap2,
108  const element_index_map& eimap3,
109  element_ranges const* mapped_element_ranges = 0,
110  bool src_is_2301 = false);
111 
121  template <typename SrcType> void convert(const SrcType& src,
122  int n1,
123  int n3,
124  const element_index_map& eimap0,
125  const element_index_map& eimap1,
126  const element_index_map& eimap2,
127  const element_index_map& eimap3,
128  element_ranges const* mapped_element_ranges = 0) {
129 
130  MPQC_ASSERT(NDIM == 4);
131  MPQC_ASSERT(mapped_element_ranges != 0);
132 
133  // determine which tiles map to the allowed element ranges in src
134  const std::vector<bool> tile0_in_erange = tiles_in_erange(range_[0], eimap0, (*mapped_element_ranges)[0]);
135  const std::vector<bool> tile1_in_erange = tiles_in_erange(range_[1], eimap1, (*mapped_element_ranges)[1]);
136  const std::vector<bool> tile2_in_erange = tiles_in_erange(range_[2], eimap2, (*mapped_element_ranges)[2]);
137  const std::vector<bool> tile3_in_erange = tiles_in_erange(range_[3], eimap3, (*mapped_element_ranges)[3]);
138 
139  // split work over tasks assuming that all tasks can access src
140  const int nproc_with_ints = info_->mem()->n();
141  const int me = info_->mem()->me();
142 
143  const size_t ntiles0 = range_[0].second - range_[0].first;
144  const size_t ntiles1 = range_[1].second - range_[1].first;
145  const size_t ntiles2 = range_[2].second - range_[2].first;
146  const size_t ntiles3 = range_[3].second - range_[3].first;
147 
148  // if espace0 is equivalent to espace1, or espace2 equivalent to espace3,
149  // can compute antisymmetric integrals using only (espace0 espace1| espace2 espace3)
150  // else also need (espace0 espace1| espace3 espace2), or an equivalent
151  const bool espace0_eq_espace1 = ((*mapped_element_ranges)[0] == (*mapped_element_ranges)[1]);
152  const bool espace2_eq_espace3 = ((*mapped_element_ranges)[2] == (*mapped_element_ranges)[3]);
153  const bool can_always_antisymmetrize = (espace0_eq_espace1 || espace2_eq_espace3);
154  MPQC_ASSERT(can_always_antisymmetrize == true); // this is likely to break the logic downstream
155  // I clearly don't understand enough at the moment
156 
157  // TODO check if equivalence of MTensor spaces is matched by their relationship in src
158 
159  // assume that src is accessible from all nodes
160  if (1) {
161 
162  // how do I determine max size of the tiles?
163  const size_t maxsize1 = info_->maxtilesize();
164  const size_t maxtilesize = maxsize1 * maxsize1 * maxsize1 * maxsize1;
165  double* data = info_->mem()->malloc_local_double(maxtilesize);
166 
167  for (long t0 = range_[0].first; t0 < range_[0].second; ++t0) {
168  const long size0 = info_->get_range(t0);
169  const long offset0 = info_->get_offset(t0);
170  const long spin0 = info_->get_spin(t0);
171  const bool in_erange0 = tile0_in_erange[t0];
172 
173  for (long t1 = std::max(range_[1].first, t0); t1 < range_[1].second; ++t1) {
174  const long size1 = info_->get_range(t1);
175  const long offset1 = info_->get_offset(t1);
176  const long spin1 = info_->get_spin(t1);
177  const bool in_erange1 = tile1_in_erange[t1];
178 
179  const bool aaaa = (spin0 == spin1);
180 
181  for (long t2 = range_[2].first; t2 < range_[2].second; ++t2) {
182  const long size2 = info_->get_range(t2);
183  const long offset2 = info_->get_offset(t2);
184  const long spin2 = info_->get_spin(t2);
185  const bool in_erange2 = tile2_in_erange[t2];
186 
187  const bool abab = ((spin0 != spin1) && (spin0 == spin2));
188  const bool abba = ((spin0 != spin1) && (spin0 != spin2));
189 
190  for (long t3 = std::max(range_[3].first, t2); t3 < range_[3].second; ++t3) {
191  const long size3 = info_->get_range(t3);
192  const long offset3 = info_->get_offset(t3);
193  const long spin3 = info_->get_spin(t3);
194  const bool in_erange3 = tile3_in_erange[t3];
195 
196  // cartesian ordinal is used as a key for hashing tiles
197  const long tile_key = (t3-range_[3].first) +
198  ntiles3 * ( (t2-range_[2].first) +
199  ntiles2 * ( (t1-range_[1].first) +
200  ntiles1 * (t0-range_[0].first)
201  )
202  );
203 
204  // since all procs can access integrals, use locality
205  if (tensor_->exists(tile_key) && tensor_->is_this_local(tile_key)) {
206 
207  if ((!in_erange0 || !in_erange1 || !in_erange2 || !in_erange3))
208  continue;
209 
210  // if erange[0] == erange[1], clearly can antisymmetrize 0 and 1
211  // if erange[2] == erange[3], clearly can antisymmetrize 2 and 3
212  bool antisymmetrize01 = espace0_eq_espace1;
213  bool antisymmetrize23 = espace2_eq_espace3;
214  // prefer to antisymmetrize 23 since can do that with 1 block of DistArray4
215  if (antisymmetrize23) antisymmetrize01 = false;
216 
217  long size = size0 * size1 * size2 * size3;
218  std::fill(data, data+size, 0.0);
219 
220  long i0123 = 0;
221  for (int i0 = 0; i0 < size0; ++i0) {
222  const int ii0 = eimap0[offset0 + i0];
223  for (int i1 = 0; i1 < size1; ++i1) {
224  const int ii1 = eimap1[offset1 + i1];
225 
226  //if (ii0 == ii1 && aaaa) continue;
227 
228  // split the work in round robin fashion by processors who have access to integrals
229  long i01;
230  if (antisymmetrize01) {
231  const long iii0 = std::max(ii0, ii1);
232  const long iii1 = std::min(ii0, ii1);
233  i01 = iii0 * (iii0+1)/2 + iii1;
234  } else {
235  i01 = ii0 * n1 + ii1;
236  }
237  const long i01_proc = i01 % nproc_with_ints;
238  if (i01_proc != me)
239  continue;
240 
241  if (antisymmetrize23) {
242 
243  for (int i2 = 0; i2 < size2; ++i2) {
244  const int ii2 = eimap2[offset2 + i2];
245  for (int i3 = 0; i3 < size3; ++i3, ++i0123) {
246  const int ii3 = eimap3[offset3 + i3];
247 
248  const int i23 = ii2 * n3 + ii3;
249  const int i32 = ii3 * n3 + ii2;
250 
251  if (abab) {
252  const double integral_0123 = src(i01,i23);
253  data[i0123] = integral_0123;
254  } else if (abba) {
255  const double integral_0132 = src(i01,i32);
256  data[i0123] = -integral_0132;
257  } else { // aaaa
258  const double integral_0123 = src(i01,i23);
259  const double integral_0132 = src(i01,i32);
260  data[i0123] = integral_0123 - integral_0132;
261  }
262 
263  }
264  }
265 
266  } // antisymmetrize23
267 
268  else if (antisymmetrize01) { // means can't antisymmetrize 2 and 3
269 
270  const int i10 = i1 * n1 + i0;
271 
272  for (int i2 = 0; i2 < size2; ++i2) {
273  const int ii2 = eimap2[offset2 + i2];
274  for (int i3 = 0; i3 < size3; ++i3, ++i0123) {
275  const int ii3 = eimap3[offset3 + i3];
276 
277  const int i23 = ii2 * n3 + ii3;
278 
279  if (abab) {
280  const double integral_0123 = src(i01,i23);
281  data[i0123] = integral_0123;
282  } else if (abba) {
283  const double integral_1023 = src(i10,i23);
284  data[i0123] = -integral_1023;
285  } else { // aaaa
286  const double integral_0123 = src(i01,i23);
287  const double integral_1023 = src(i10,i23);
288  data[i0123] = integral_0123 - integral_1023;
289  }
290 
291  }
292  }
293 
294  } // antisymmetrize01
295 
296  else { // do not antisymmetrize
297  MPQC_ASSERT(false); // should not happen, but only SMITH knows :-)
298  }
299 
300  }
301  }
302 
303  tensor_->put_block(tile_key, data);
304 
305 #if 0
306  const double sum = std::accumulate(data, data+size, 0.0);
307  ExEnv::out0() << "tiles = (" << t0 << "," << t1 << "," << t2 << "," << t3
308  << ") key = " << tile_key << " sum = " << sum << std::endl;
309 
310  for(int i=0; i<size; ++i) {
311  ExEnv::out0() << "data[" << i << "] = " << data[i] << std::endl;
312  }
313 #endif
314  } // if this task will process this tile
315 
316  } // t3
317  } // t2
318  } // t1
319  } // t0
320 
321  info_->mem()->free_local_double(data);
322 
323  } // if this task can access the integrals
324 
325  } // MTensor<4>::convert() from RefSCMatrix or RefSymmSCMatrix
326 
336  template <typename SrcType> void convert(const SrcType& src0, const SrcType& src1,
337  int n0,
338  int n1,
339  int n2,
340  int n3,
341  const element_index_map& eimap0,
342  const element_index_map& eimap1,
343  const element_index_map& eimap2,
344  const element_index_map& eimap3,
345  element_ranges const* mapped_element_ranges = 0,
346  bool src1_is_1023 = false) {
347 
348  MPQC_ASSERT(NDIM == 4);
349  MPQC_ASSERT(mapped_element_ranges != 0);
350 
351  // determine which tiles map to the allowed element ranges in src
352  const std::vector<bool> tile0_in_erange = tiles_in_erange(range_[0], eimap0, (*mapped_element_ranges)[0]);
353  const std::vector<bool> tile1_in_erange = tiles_in_erange(range_[1], eimap1, (*mapped_element_ranges)[1]);
354  const std::vector<bool> tile2_in_erange = tiles_in_erange(range_[2], eimap2, (*mapped_element_ranges)[2]);
355  const std::vector<bool> tile3_in_erange = tiles_in_erange(range_[3], eimap3, (*mapped_element_ranges)[3]);
356 
357  // split work over tasks assuming that all tasks can access src
358  const int nproc_with_ints = info_->mem()->n();
359  const int me = info_->mem()->me();
360 
361  const size_t ntiles0 = range_[0].second - range_[0].first;
362  const size_t ntiles1 = range_[1].second - range_[1].first;
363  const size_t ntiles2 = range_[2].second - range_[2].first;
364  const size_t ntiles3 = range_[3].second - range_[3].first;
365 
366  // if espace0 is equivalent to espace1, or espace2 equivalent to espace3,
367  // can compute antisymmetric integrals using only (espace0 espace1| espace2 espace3)
368  // else also need (espace0 espace1| espace3 espace2), or an equivalent
369  const bool espace0_eq_espace1 = ((*mapped_element_ranges)[0] == (*mapped_element_ranges)[1]);
370  const bool espace2_eq_espace3 = ((*mapped_element_ranges)[2] == (*mapped_element_ranges)[3]);
371  const bool can_always_antisymmetrize = (espace0_eq_espace1 || espace2_eq_espace3);
372  MPQC_ASSERT(can_always_antisymmetrize == true); // this is likely to break the logic downstream
373  // I clearly don't understand enough at the moment
374 
375  // TODO check if equivalence of MTensor spaces is matched by their relationship in src
376 
377  // assume that src is accessible from all nodes
378  if (1) {
379 
380  // how do I determine max size of the tiles?
381  const size_t maxsize1 = info_->maxtilesize();
382  const size_t maxtilesize = maxsize1 * maxsize1 * maxsize1 * maxsize1;
383  double* data = info_->mem()->malloc_local_double(maxtilesize);
384 
385  for (long t0 = range_[0].first; t0 < range_[0].second; ++t0) {
386  const long size0 = info_->get_range(t0);
387  const long offset0 = info_->get_offset(t0);
388  const long spin0 = info_->get_spin(t0);
389  const bool in_erange0 = tile0_in_erange[t0];
390 
391  for (long t1 = std::max(range_[1].first, t0); t1 < range_[1].second; ++t1) {
392  const long size1 = info_->get_range(t1);
393  const long offset1 = info_->get_offset(t1);
394  const long spin1 = info_->get_spin(t1);
395  const bool in_erange1 = tile1_in_erange[t1];
396 
397  const bool aaaa = (spin0 == spin1);
398 
399  for (long t2 = range_[2].first; t2 < range_[2].second; ++t2) {
400  const long size2 = info_->get_range(t2);
401  const long offset2 = info_->get_offset(t2);
402  const long spin2 = info_->get_spin(t2);
403  const bool in_erange2 = tile2_in_erange[t2];
404 
405  const bool abab = ((spin0 != spin1) && (spin0 == spin2));
406  const bool abba = ((spin0 != spin1) && (spin0 != spin2));
407 
408  for (long t3 = std::max(range_[3].first, t2); t3 < range_[3].second; ++t3) {
409  const long size3 = info_->get_range(t3);
410  const long offset3 = info_->get_offset(t3);
411  const long spin3 = info_->get_spin(t3);
412  const bool in_erange3 = tile3_in_erange[t3];
413 
414  // cartesian ordinal is used as a key for hashing tiles
415  const long tile_key = (t3-range_[3].first) +
416  ntiles3 * ( (t2-range_[2].first) +
417  ntiles2 * ( (t1-range_[1].first) +
418  ntiles1 * (t0-range_[0].first)
419  )
420  );
421 
422  // since all procs can access integrals, use locality
423  if (tensor_->exists(tile_key) && tensor_->is_this_local(tile_key)) {
424 
425  if ((!in_erange0 || !in_erange1 || !in_erange2 || !in_erange3))
426  continue;
427 
428  // if erange[0] == erange[1], clearly can antisymmetrize 0 and 1
429  // if erange[2] == erange[3], clearly can antisymmetrize 2 and 3
430  bool antisymmetrize01 = espace0_eq_espace1;
431  bool antisymmetrize23 = espace2_eq_espace3;
432  // prefer to antisymmetrize 23 since can do that with 1 block of DistArray4
433  if (antisymmetrize23) antisymmetrize01 = false;
434 
435  long size = size0 * size1 * size2 * size3;
436  std::fill(data, data+size, 0.0);
437 
438  long i0123 = 0;
439  for (int i0 = 0; i0 < size0; ++i0) {
440  const int ii0 = eimap0[offset0 + i0];
441  for (int i1 = 0; i1 < size1; ++i1) {
442  const int ii1 = eimap1[offset1 + i1];
443 
444  //if (ii0 == ii1 && aaaa) continue;
445 
446  // split the work in round robin fashion by processors who have access to integrals
447  long i01;
448  if (antisymmetrize01) {
449  const long iii0 = std::max(ii0, ii1);
450  const long iii1 = std::min(ii0, ii1);
451  i01 = iii0 * (iii0+1)/2 + iii1;
452  } else {
453  i01 = ii0 * n1 + ii1;
454  }
455  const long i01_proc = i01 % nproc_with_ints;
456  if (i01_proc != me)
457  continue;
458 
459  if (antisymmetrize23) {
460 
461  for (int i2 = 0; i2 < size2; ++i2) {
462  const int ii2 = eimap2[offset2 + i2];
463  for (int i3 = 0; i3 < size3; ++i3, ++i0123) {
464  const int ii3 = eimap3[offset3 + i3];
465 
466  const int i23 = ii2 * n3 + ii3;
467  const int i32 = ii3 * n2 + ii2;
468 
469  if (abab) {
470  const double integral_0123 = src0(i01,i23);
471  data[i0123] = integral_0123;
472  } else if (abba) {
473  const double integral_0132 = src1(i01,i32);
474  data[i0123] = -integral_0132;
475  } else { // aaaa
476  const double integral_0123 = src0(i01,i23);
477  const double integral_0132 = src1(i01,i32);
478  data[i0123] = integral_0123 - integral_0132;
479  }
480 
481  }
482  }
483 
484  } // antisymmetrize23
485 
486  else if (antisymmetrize01) { // means can't antisymmetrize 2 and 3
487 
488  const int i10 = i1 * n0 + i0;
489 
490  for (int i2 = 0; i2 < size2; ++i2) {
491  const int ii2 = eimap2[offset2 + i2];
492  for (int i3 = 0; i3 < size3; ++i3, ++i0123) {
493  const int ii3 = eimap3[offset3 + i3];
494 
495  const int i23 = ii2 * n3 + ii3;
496 
497  if (abab) {
498  const double integral_0123 = src0(i01,i23);
499  data[i0123] = integral_0123;
500  } else if (abba) {
501  const double integral_1023 = src1(i10,i23);
502  data[i0123] = -integral_1023;
503  } else { // aaaa
504  const double integral_0123 = src0(i01,i23);
505  const double integral_1023 = src1(i10,i23);
506  data[i0123] = integral_0123 - integral_1023;
507  }
508 
509  }
510  }
511 
512  } // antisymmetrize01
513 
514  else { // do not antisymmetrize
515  MPQC_ASSERT(false); // should not happen, but only SMITH knows :-)
516  }
517 
518  }
519  }
520 
521  tensor_->put_block(tile_key, data);
522 
523 #if 0
524  const double sum = std::accumulate(data, data+size, 0.0);
525  ExEnv::out0() << "tiles = (" << t0 << "," << t1 << "," << t2 << "," << t3
526  << ") key = " << tile_key << " sum = " << sum << std::endl;
527 
528  for(int i=0; i<size; ++i) {
529  ExEnv::out0() << "data[" << i << "] = " << data[i] << std::endl;
530  }
531 #endif
532  } // if this task will process this tile
533 
534  } // t3
535  } // t2
536  } // t1
537  } // t0
538 
539  info_->mem()->free_local_double(data);
540 
541  } // if this task can access the integrals
542 
543  } // MTensor<4>::convert() from RefSCMatrix or RefSymmSCMatrix
544 
551  template <typename SrcType> void convert(const SrcType& src,
552  const element_index_map& eimap0,
553  const element_index_map& eimap1,
554  element_ranges const* mapped_element_ranges = 0) {
555 
556  MPQC_ASSERT(NDIM == 2);
557 
558  // determine which tiles map to the allowed element ranges in src
559  const std::vector<bool> tile0_in_erange = tiles_in_erange(range_[0], eimap0, (*mapped_element_ranges)[0]);
560  const std::vector<bool> tile1_in_erange = tiles_in_erange(range_[1], eimap1, (*mapped_element_ranges)[1]);
561 
562  // split work over tasks assuming that all tasks can access src
563  const int nproc_with_ints = info_->mem()->n();
564  const int me = info_->mem()->me();
565 
566  const size_t ntiles0 = range_[0].second - range_[0].first;
567  const size_t ntiles1 = range_[1].second - range_[1].first;
568 
569  // check that the needed index ranges are present in src
570  #if 0 // assume that the user knows what she's doing
571  element_range src_range[2];
572  src_range[0] = this->map_range(range_[0], eimap0);
573  src_range[1] = this->map_range(range_[1], eimap1);
574  MPQC_ASSERT(src_range[0].second <= src->ni());
575  MPQC_ASSERT(src_range[1].second <= src->nj());
576  #endif
577 
578  // TODO check if equivalence of MTensor spaces is matched by their relationship in src
579 
580  // assume that src is accessible from all nodes
581  if (1) {
582 
583  // how do I determine max size of the tiles?
584  const size_t maxsize1 = info_->maxtilesize();
585  const size_t maxtilesize = maxsize1 * maxsize1;
586  double* data = info_->mem()->malloc_local_double(maxtilesize);
587 
588  for (long t0 = range_[0].first; t0 < range_[0].second; ++t0) {
589  const long size0 = info_->get_range(t0);
590  const long offset0 = info_->get_offset(t0);
591  const long spin0 = info_->get_spin(t0);
592  const bool in_erange0 = tile0_in_erange[t0];
593 
594  for (long t1 = range_[1].first; t1 < range_[1].second; ++t1) {
595  const long size1 = info_->get_range(t1);
596  const long offset1 = info_->get_offset(t1);
597  const long spin1 = info_->get_spin(t1);
598  const bool in_erange1 = tile1_in_erange[t1];
599 
600  // cartesian ordinal is used as a key for hashing tiles
601  const long tile_key = (t1-range_[1].first) + ntiles1 * (t0-range_[0].first);
602 
603  // since all procs can access integrals, use locality
604  if (tensor_->exists(tile_key) && tensor_->is_this_local(tile_key)) {
605 
606  if ((!in_erange0 || !in_erange1))
607  continue;
608 
609  long size = size0 * size1;
610  std::fill(data, data+size, 0.0);
611 
612  long i01 = 0;
613  for (int i0 = 0; i0 < size0; ++i0) {
614  const int ii0 = eimap0[offset0 + i0];
615  for (int i1 = 0; i1 < size1; ++i1, ++i01) {
616  const int ii1 = eimap1[offset1 + i1];
617 
618  data[i01] = src(ii0, ii1);
619  }
620  }
621 
622  tensor_->add_block(tile_key, data);
623 
624 #if 0
625  const double sum = std::accumulate(data, data+size, 0.0);
626  ExEnv::out0() << "tiles = (" << t0 << "," << t1
627  << ") key = " << tile_key << " sum = " << sum << endl;
628 #endif
629  } // if this task will process this tile
630 
631  } // t1
632  } // t0
633 
634  info_->mem()->free_local_double(data);
635 
636  } // if this task can access the integrals
637 
638  } // MTensor<2>::convert() from RefSCMatrix or RefSymmSCMatrix
639 
640  void print(const std::string& label, std::ostream& os = ExEnv::out0()) const {
641  tensor_->print(label, os);
642  }
643 
644  private:
645  Info const* info_;
646  Tensor* tensor_;
647  tile_ranges range_; // range of tiles that defines this tensor. See info_ for tiling info
648 
649  // given a range of indices, return pair<mmin,mmax> where mmin and mmax are the minimum and maximum mapped indices
650  element_range map_range(const element_range& rng,
651  const element_index_map& eimap) {
652  element_range result(eimap[rng.first],
653  eimap[rng.first]);
654  for(element_index i=rng.first; i != rng.second; ++i) {
655  const element_index ii = eimap[i];
656  result.first = std::min(result.first, ii);
657  result.second = std::max(result.second, ii);
658  }
659  return result;
660  }
661 
662  // given tile range [start, fence), element map, and element_range, result[t] is true if all elements within
663  // the tile map to inside element_range
664  std::vector<bool>
665  tiles_in_erange(const tile_range& range,
666  const element_index_map& eimap,
667  const element_range& erange) {
668 
669  std::vector<bool> result(range.second, false);
670  for(long t=range.first; t!=range.second; ++t) {
671  const long start = info_->get_offset(t);
672  const long size = info_->get_range(t);
673  const long fence = start + size;
674  const element_range tile(start, fence);
675  // map elements
676  const element_range mapped_tile = map_range(tile, eimap);
677  result[t] = in(mapped_tile,erange);
678  }
679  return result;
680  }
681 
682  };
683 
684 } // end of namespace sc
685 
686 #include <chemistry/qc/ccr12/mtensor.timpl.h>
687 
688 #endif // end of header guard
689 
690 // Local Variables:
691 // mode: c++
692 // c-file-style: "CLJ-CONDENSED"
693 // End:
sc::Tensor::add_block
void add_block(long tag, double *data)
add a block to the distributed file (non-blocking); double* data will be destroyed
sc::Tensor::exists
bool exists(long tag) const
does this block exist?
sc::Tensor
Definition: tensor.h:49
sc::SCDimension
The SCDimension class is used to determine the size and blocking of matrices.
Definition: dim.h:105
mpqc::range
Definition: range.hpp:23
sc::RefSCMatrix
The RefSCMatrix class is a smart pointer to an SCMatrix specialization.
Definition: matrix.h:135
sc::SCMatrixKit::matrix
virtual SCMatrix * matrix(const RefSCDimension &, const RefSCDimension &)=0
Given the dimensions, create matrices or vectors.
sc::Ref
A template class that maintains references counts.
Definition: ref.h:361
sc::MTensor::convert
void convert(const SrcType &src0, const SrcType &src1, int n0, int n1, int n2, int n3, const element_index_map &eimap0, const element_index_map &eimap1, const element_index_map &eimap2, const element_index_map &eimap3, element_ranges const *mapped_element_ranges=0, bool src1_is_1023=false)
copies contents of (src0, src1) to this, where src0 = (01|23) and src1 = (01|32), if src1_is_1023 == ...
Definition: mtensor.h:336
sc::MTensor::convert
void convert(const SrcType &src, const element_index_map &eimap0, const element_index_map &eimap1, element_ranges const *mapped_element_ranges=0)
copies contents of src to this.
Definition: mtensor.h:551
sc::MTensor::convert
void convert(const Ref< DistArray4 > &src, unsigned int tbtype, const element_index_map &eimap0, const element_index_map &eimap1, const element_index_map &eimap2, const element_index_map &eimap3, element_ranges const *mapped_element_ranges=0, bool src_is_2301=false)
copies contents of src to this.
sc::in
bool in(const std::pair< I, I > &r, const std::pair< I, I > &range)
return true if r is contained in range defined as pair<start,fence>, i.e. [start, fence)
Definition: mtensor.h:58
sc::MTensor::convert
void convert(const SrcType &src, int n1, int n3, const element_index_map &eimap0, const element_index_map &eimap1, const element_index_map &eimap2, const element_index_map &eimap3, element_ranges const *mapped_element_ranges=0)
copies contents of src to this.
Definition: mtensor.h:121
sc::SCMatrixKit::default_matrixkit
static SCMatrixKit * default_matrixkit()
This returns the default matrix kit.
sc::Tensor::print
void print(const std::string &label, std::ostream &os=ExEnv::out0()) const
print
sc::Tensor::put_block
void put_block(long tag, double *data)
put a block to the distributed file (non-blocking)
sc::Tensor::is_this_local
bool is_this_local(long tag)
returns if a block associated with long tag resides in a local memory or not
sc::intersect
std::pair< I, I > intersect(const std::pair< I, I > &range1, const std::pair< I, I > &range2)
return intersect of two ranges defined as pair<start,fence>, i.e. [start, fence)
Definition: mtensor.h:50
sc::ExEnv::out0
static std::ostream & out0()
Return an ostream that writes from node 0.
sc::MTensor
Tensor metadata is implicit; MTensor is Tensor + metadata.
Definition: mtensor.h:70
sc::MTensor::MTensor
MTensor(Info const *info, Tensor *tensor, const tile_ranges &range)
Definition: mtensor.h:85
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
sc::CCR12_Info
CCR12_Info is the compilation of members that are used in CC and CC-R12 methods.
Definition: ccr12_info.h:50

Generated at Sun Jan 26 2020 23:23:57 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.