MPQC  3.0.0-alpha
sma.h
1 
2 /*
3  * Copyright 2009 Sandia Corporation. Under the terms of Contract
4  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
5  * retains certain rights in this software.
6  *
7  * This file is a part of the MPQC LMP2 library.
8  *
9  * The MPQC LMP2 library is free software: you can redistribute it
10  * and/or modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation, either
12  * version 3 of the License, or (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful, but
15  * WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with this program. If not, see
21  * <http://www.gnu.org/licenses/>.
22  *
23  */
24 
25 #ifndef _chemistry_qc_lmp2_sma_h
26 #define _chemistry_qc_lmp2_sma_h
27 
28 //#define USE_HASH
29 //#define USE_STL_MULTIMAP 1
30 
31 #include <chemistry/qc/basis/basis.h>
32 #include <chemistry/qc/basis/integral.h>
33 #include <util/misc/autovec.h>
34 #include <util/misc/scint.h>
35 #include <util/misc/exenv.h>
36 #include <util/state/statein.h>
37 #include <util/state/stateout.h>
38 #include <util/misc/regtime.h>
39 #include <util/misc/scexception.h>
40 
41 #ifndef TWO_INDEX_SPECIALIZATIONS
42 # define THREE_INDEX_SPECIALIZATIONS 1
43 #endif
44 
45 #ifndef THREE_INDEX_SPECIALIZATIONS
46 # define THREE_INDEX_SPECIALIZATIONS 1
47 #endif
48 
49 #ifndef FOUR_INDEX_SPECIALIZATIONS
50 # define FOUR_INDEX_SPECIALIZATIONS 1
51 #endif
52 
53 #ifndef USE_STL_MULTIMAP
54 # define USE_STL_MULTIMAP 0
55 #endif
56 
57 #if ! USE_STL_MULTIMAP
58 # include <chemistry/qc/lmp2/avlmmap.h>
59 #endif
60 
61 #include <math.h>
62 #include <float.h>
63 #include <vector>
64 #include <map>
65 #ifdef USE_HASH
66 # include <ext/hash_map>
67 #endif
68 #include <set>
69 #include <stdexcept>
70 #include <string>
71 #ifdef WORDS_BIGENDIAN
72 # define r2(i) (i)
73 # define r4(i) (i)
74 #else
75 // For little endian machines, the order of indices in BlockInfo<4>
76 // must be reversed so a fast comparison on an uint64_t can be used
77 // to order blocks.
78 // BlockInfo<3> uses r4 and zeros out the last byte.
79 // BlockInfo<2> uses r2.
80 //# define r4(i) (3-(i))
81 # define r2(i) ((~(i))&1)
82 # define r4(i) ((~(i))&3)
83 #endif
84 
85 namespace sc {
86 
87 namespace sma2 {
88 
93  class Range {
94  private:
95  int nindex_;
96  std::vector<int> block_size_;
97  std::vector<int> block_offset_;
98  std::vector<int> index_to_block_;
99  std::vector<int> function_order_;
100  std::vector<int> range_order_;
101 
102  int max_block_size_;
103 
104  void init_offsets();
105  void init_extent_blocking(const sc::Ref<sc::GaussianBasisSet> &bs);
106  void init_order(int nbasis);
107  public:
108  enum BlockingMethod { AtomBlocking,
109  ShellBlocking,
110  FunctionBlocking,
111  ExtentBlocking };
112 
113  Range(const sc::Ref<sc::GaussianBasisSet> &,
114  BlockingMethod b = ShellBlocking, int blocksize = 1);
115  Range(int nindex, int block_size);
116  Range(const std::vector<int> & block_size);
117  Range();
118 
119  void init(const sc::Ref<sc::GaussianBasisSet> &,
120  BlockingMethod b = ShellBlocking, int blocksize = 1);
121  void init(int nindex, int block_size);
122  void clear();
123 
125  int nindex() const { return nindex_; }
127  int nblock() const { return block_size_.size(); }
129  int block_size(int i) const { return block_size_[i]; }
131  int max_block_size() const { return max_block_size_; }
133  int block_offset(int i) const { return block_offset_[i]; }
135  int index_to_block(int i) const { return index_to_block_[i]; }
138  int index_to_offset(int i) const {
139  return i - block_offset_[index_to_block_[i]];
140  }
145  int basis_index_to_range_index(int i) const;
146  int range_index_to_basis_index(int i) const;
147  bool operator == (const Range &r) const;
148  bool operator != (const Range &r) const { return ! operator == (r); }
150  bool all_size_one() const;
151  void print(std::ostream&o=sc::ExEnv::out0()) const;
152  void write(sc::StateOut&) const;
153  void read(sc::StateIn&);
154  };
155 
160  class IndexList {
161  private:
162  std::vector<int> indices_;
163  public:
164  IndexList();
165  IndexList(int);
166  IndexList(int,int);
167  IndexList(int,int,int);
168  IndexList(int,int,int,int);
169  IndexList(const IndexList &);
172  IndexList(const IndexList &il1, const IndexList &il2);
173  IndexList(const std::vector<int> &);
174  IndexList reverse_mapping() const;
175  void append_additional_indices(const IndexList &);
176  int &i(int a) { return indices_[a]; }
177  const int &i(int a) const { return indices_[a]; }
178  int n() const { return indices_.size(); }
179  void set_n(int n) { indices_.resize(n); }
180  static IndexList identity(int n);
182  bool is_identity() const;
184  bool is_identity_permutation() const;
185  void print(std::ostream &o=sc::ExEnv::outn()) const;
186  bool operator>(const IndexList&il)const{return indices_>il.indices_;}
187  bool operator<(const IndexList&il)const{return indices_<il.indices_;}
188  bool operator==(const IndexList&il)const{return indices_==il.indices_;}
189  bool operator!=(const IndexList&il)const{return indices_!=il.indices_;}
190  bool operator>=(const IndexList&il)const{return indices_>=il.indices_;}
191  bool operator<=(const IndexList&il)const{return indices_<=il.indices_;}
192  };
193  std::ostream &operator << (std::ostream &, const IndexList &);
194 
195  typedef unsigned short bi_t; // Block index type.
196 
199  template <int N>
200  class BlockInfo {
201  private:
202 #ifdef USE_BOUND
203  mutable double bound_;
204 #endif
205  bi_t block_[N];
206  public:
207  BlockInfo() {
208 #ifdef USE_BOUND
209  bound_ = DBL_MAX;
210 #endif
211  }
214  BlockInfo(const std::vector<bi_t> &v) {
215 #ifdef USE_BOUND
216  bound_ = DBL_MAX;
217 #endif
218  for (int i=0; i<N && i <v.size(); i++) block_[i] = v[i];
219  }
223  BlockInfo(const BlockInfo<N> &b, const IndexList &l) {
224 #ifdef USE_BOUND
225  bound_ = DBL_MAX;
226 #endif
227  for (int i=0; i<l.n(); i++) block_[i] = b.block_[l.i(i)];
228  }
229  template <int NB>
230  BlockInfo(const IndexList &l,
231  const BlockInfo<NB> &b,
232  const IndexList &lb) {
233 #ifdef USE_BOUND
234  bound_ = DBL_MAX;
235 #endif
236  for (int i=0; i<l.n(); i++) block_[l.i(i)] = b.block(lb.i(i));
237  }
239  void zero() { for (int i=0; i<N; i++) block_[i] = 0; }
240 #ifdef USE_BOUND
241  double bound() const { return bound_; }
242  void set_bound(double b) const { bound_ = b; }
243 #endif
244  bi_t &block(int i) { return block_[i]; }
245  const bi_t &block(int i) const { return block_[i]; }
247  unsigned int size(const Range *indices) const {
248  unsigned int r = 1;
249  for (int i=0; i<N; i++) r *= indices[i].block_size(block_[i]);
250  return r;
251  }
254  unsigned int subset_size(const Range *indices,
255  const IndexList &indexlist) const {
256  unsigned int r = 1;
257  for (int i=0; i<indexlist.n(); i++) {
258  int index = indexlist.i(i);
259  r *= indices[index].block_size(block_[index]);
260  }
261  return r;
262  }
267  template <int N2>
268  void assign_blocks(const IndexList &il,
269  const BlockInfo<N2> &bi2, const IndexList &il2) {
270  for (int i=0; i<il.n(); i++) {
271  block_[il.i(i)] = bi2.block(il2.i(i));
272  }
273  }
277  template <int N2>
278  void assign_blocks(const IndexList &il,
279  const BlockInfo<N2> &bi2) {
280  for (int i=0; i<il.n(); i++) {
281  block_[il.i(i)] = bi2.block(i);
282  }
283  }
287  template <int N2>
288  bool equiv_blocks(const IndexList &il,
289  const BlockInfo<N2> &bi2) {
290  for (int i=0; i<il.n(); i++) {
291  if (block_[il.i(i)] != bi2.block(i)) return false;
292  }
293  return true;
294  }
295  void print(std::ostream &o=sc::ExEnv::outn()) const;
296  void print_block_sizes(const Range *indices,
297  std::ostream &o=sc::ExEnv::outn()) const {
298  o << "{";
299  for (int i=0; i<N; i++) {
300  if (i) o << " ";
301  o << indices[i].block_size(block_[i]);
302  }
303  o << "}";
304  }
305  void write(sc::StateOut& so) const {
306 #ifdef USE_BOUND
307  so.put(bound_);
308 #endif
309  for (int i=0; i<N; i++) {
310  so.put(int(block_[i]));
311  }
312  }
313  void read(sc::StateIn& si) {
314 #ifdef USE_BOUND
315  si.get(bound_);
316 #endif
317  for (int i=0; i<N; i++) {
318  int b;
319  si.get(b);
320  block_[i] = b;
321  }
322  }
323  };
324  template <int N>
325  inline std::ostream& operator << (std::ostream&o, const BlockInfo<N> &b)
326  {
327  b.print(o);
328  return o;
329  }
330 
332  template <int N>
333  class IndicesLess {
334  public:
335  bool operator() (const BlockInfo<N>&b1, const BlockInfo<N>&b2) const {
336  int b1b = b1.block(0), b2b = b2.block(0);
337  if (b1b < b2b) return true;
338  if (b1b > b2b) return false;
339  for (int i=1; i<N; i++) {
340  b1b = b1.block(i); b2b = b2.block(i);
341  if (b1b < b2b) return true;
342  if (b1b > b2b) return false;
343  }
344  return false;
345  }
346  int compare(const BlockInfo<N>&b1, const BlockInfo<N>&b2) const {
347  int b1b = b1.block(0), b2b = b2.block(0);
348  if (b1b < b2b) return -1;
349  if (b1b > b2b) return 1;
350  for (int i=1; i<N; i++) {
351  b1b = b1.block(i); b2b = b2.block(i);
352  if (b1b < b2b) return -1;
353  if (b1b > b2b) return 1;
354  }
355  return 0;
356  }
357  };
358 
359  // Specializations for 0 indices
360  template <>
361  class BlockInfo<0> {
362  private:
363 #ifdef USE_BOUND
364  mutable double bound_;
365 #endif
366  public:
367  friend class IndicesLess<0>;
368  public:
369  BlockInfo() {
370 #ifdef USE_BOUND
371  bound_ = DBL_MAX;
372 #endif
373  }
374  BlockInfo(const std::vector<bi_t> &v) {
375 #ifdef USE_BOUND
376  bound_ = DBL_MAX;
377 #endif
378  }
379  BlockInfo(const BlockInfo<2> &b, const IndexList &l) {
380 #ifdef USE_BOUND
381  bound_ = DBL_MAX;
382 #endif
383  }
384  template <int NB>
385  BlockInfo(const IndexList &l,
386  const BlockInfo<NB> &b,
387  const IndexList &lb) {
388 #ifdef USE_BOUND
389  bound_ = DBL_MAX;
390 #endif
391  }
392 #ifdef USE_BOUND
393  double bound() const { return bound_; }
394  void set_bound(double b) const { bound_ = b; }
395 #endif
396  bi_t &block(int i) {
397  throw sc::ProgrammingError("BlockInfo<0>:block: can not be called",
398  __FILE__, __LINE__);
399  }
400  const bi_t &block(int i) const {
401  throw sc::ProgrammingError("BlockInfo<0>:block: can not be called",
402  __FILE__, __LINE__);
403  }
405  unsigned int size(const Range *indices) const {
406  return 1;
407  }
410  unsigned int subset_size(const Range *indices,
411  const IndexList &indexlist) const {
412  return 1;
413  }
418  template <int N2>
419  void assign_blocks(const IndexList &il,
420  const BlockInfo<N2> &bi2, const IndexList &il2) {
421  }
425  template <int N2>
426  void assign_blocks(const IndexList &il,
427  const BlockInfo<N2> &bi2) {
428  }
432  template <int N2>
433  bool equiv_blocks(const IndexList &il,
434  const BlockInfo<N2> &bi2) {
435  return true;
436  }
438  void zero() {}
439  void print(std::ostream &o=sc::ExEnv::outn()) const {
440  o << "{}";
441  }
442  void print_block_sizes(const Range *indices,
443  std::ostream &o=sc::ExEnv::outn()) const {
444  o << "{}";
445  }
446  void write(sc::StateOut& so) const {
447 #ifdef USE_BOUND
448  so.put(bound_);
449 #endif
450  }
451  void read(sc::StateIn& si) {
452 #ifdef USE_BOUND
453  si.get(bound_);
454 #endif
455  }
456  };
458  template <>
459  class IndicesLess<0> {
460  public:
461  bool operator() (const BlockInfo<0>&b1, const BlockInfo<0>&b2) const {
462  return false;
463  }
464  int compare(const BlockInfo<0>&b1, const BlockInfo<0>&b2) const {
465  return 0;
466  }
467  };
468 
469 #if TWO_INDEX_SPECIALIZATIONS
470  // Specializations for 2 indices
471  template <>
472  class BlockInfo<2> {
473  private:
474 #ifdef USE_BOUND
475  mutable double bound_;
476 #endif
477  public:
478  union {
479  sc::sc_uint32_t i;
480  sc::sc_uint16_t s[2];
481  sc::sc_uint8_t c[4];
482  } b_;
483  friend class IndicesLess<2>;
484  public:
485  BlockInfo() {
486 #ifdef USE_BOUND
487  bound_ = DBL_MAX;
488 #endif
489  }
490  BlockInfo(const std::vector<bi_t> &v) {
491 #ifdef USE_BOUND
492  bound_ = DBL_MAX;
493 #endif
494  for (int i=0; i<2 && i<v.size(); i++) b_.s[r2(i)] = v[i];
495  }
496  BlockInfo(const BlockInfo<2> &b, const IndexList &l) {
497 #ifdef USE_BOUND
498  bound_ = DBL_MAX;
499 #endif
500  for (int i=0; i<l.n(); i++) b_.s[r2(i)] = b.b_.s[r2(l.i(i))];
501  }
502  template <int NB>
503  BlockInfo(const IndexList &l,
504  const BlockInfo<NB> &b,
505  const IndexList &lb) {
506 #ifdef USE_BOUND
507  bound_ = DBL_MAX;
508 #endif
509  for (int i=0; i<l.n(); i++) b_.s[r2(l.i(i))] = b.block(lb.i(i));
510  }
511 #ifdef USE_BOUND
512  double bound() const { return bound_; }
513  void set_bound(double b) const { bound_ = b; }
514 #endif
515  bi_t &block(int i) { return b_.s[r2(i)]; }
516  const bi_t &block(int i) const { return b_.s[r2(i)]; }
518  unsigned int size(const Range *indices) const {
519  unsigned int r = 1;
520  for (int i=0; i<2; i++) r *= indices[i].block_size(b_.s[r2(i)]);
521  return r;
522  }
525  unsigned int subset_size(const Range *indices,
526  const IndexList &indexlist) const {
527  unsigned int r = 1;
528  for (int i=0; i<indexlist.n(); i++) {
529  int index = indexlist.i(i);
530  r *= indices[index].block_size(b_.s[r2(index)]);
531  }
532  return r;
533  }
538  template <int N2>
539  void assign_blocks(const IndexList &il,
540  const BlockInfo<N2> &bi2, const IndexList &il2) {
541  for (int i=0; i<il.n(); i++) {
542  b_.s[r2(il.i(i))] = bi2.block(il2.i(i));
543  }
544  }
548  template <int N2>
549  void assign_blocks(const IndexList &il,
550  const BlockInfo<N2> &bi2) {
551  for (int i=0; i<il.n(); i++) {
552  b_.s[r2(il.i(i))] = bi2.block(i);
553  }
554  }
558  template <int N2>
559  bool equiv_blocks(const IndexList &il,
560  const BlockInfo<N2> &bi2) {
561  for (int i=0; i<il.n(); i++) {
562  if (b_.s[r2(il.i(i))] != bi2.block(i)) return false;
563  }
564  return true;
565  }
567  void zero() { b_.i = 0; }
568  void print(std::ostream &o=sc::ExEnv::outn()) const {
569  o << "{";
570  for (int i=0; i<2; i++) {
571  if (i!=0) o << " ";
572  o << block(i);
573  }
574  o << "}";
575  }
576  void print_block_sizes(const Range *indices,
577  std::ostream &o=sc::ExEnv::outn()) const {
578  o << "{";
579  for (int i=0; i<2; i++) {
580  if (i) o << " ";
581  o << indices[i].block_size(b_.s[r2(i)]);
582  }
583  o << "}";
584  }
585  void write(sc::StateOut& so) const {
586 #ifdef USE_BOUND
587  so.put(bound_);
588 #endif
589  for (int i=0; i<2; i++) {
590  so.put(int(block(i)));
591  }
592  }
593  void read(sc::StateIn& si) {
594 #ifdef USE_BOUND
595  si.get(bound_);
596 #endif
597  for (int i=0; i<2; i++) {
598  int b;
599  si.get(b);
600  block(i) = b;
601  }
602  }
603  };
604 
605  template<> inline
606  bool IndicesLess<2>::operator() (const BlockInfo<2>&b1,
607  const BlockInfo<2>&b2) const
608  {
609  if (b1.b_.i < b2.b_.i) return true;
610  return false;
611  }
612 
613  template<> inline
614  int IndicesLess<2>::compare(const BlockInfo<2>&b1,
615  const BlockInfo<2>&b2) const
616  {
617  if (b1.b_.i < b2.b_.i) return -1;
618  else if (b1.b_.i > b2.b_.i) return 1;
619  return 0;
620  }
621 #endif // TWO_INDEX_SPECIALIZATIONS
622 
623 #if THREE_INDEX_SPECIALIZATIONS
624  // Specializations for 3 indices
625  template <>
626  class BlockInfo<3> {
627  private:
628 #ifdef USE_BOUND
629  mutable double bound_;
630 #endif
631  public:
632  union {
633  sc::sc_uint64_t l;
634  sc::sc_uint32_t i[2];
635  sc::sc_uint16_t s[4];
636  sc::sc_uint8_t c[8];
637  } b_;
638  friend class IndicesLess<3>;
639  public:
640  BlockInfo() {
641 #ifdef USE_BOUND
642  bound_ = DBL_MAX;
643 #endif
644  b_.l = 0;
645  }
646  BlockInfo(const std::vector<bi_t> &v) {
647 #ifdef USE_BOUND
648  bound_ = DBL_MAX;
649 #endif
650  b_.l = 0;
651  for (int i=0; i<3 && i<v.size(); i++) b_.s[r4(i)] = v[i];
652  }
653  BlockInfo(const BlockInfo<3> &b, const IndexList &l) {
654 #ifdef USE_BOUND
655  bound_ = DBL_MAX;
656 #endif
657  b_.l = 0;
658  for (int i=0; i<l.n(); i++) b_.s[r4(i)] = b.b_.s[r4(l.i(i))];
659  }
660  template <int NB>
661  BlockInfo(const IndexList &l,
662  const BlockInfo<NB> &b,
663  const IndexList &lb) {
664 #ifdef USE_BOUND
665  bound_ = DBL_MAX;
666 #endif
667  b_.l = 0;
668  for (int i=0; i<l.n(); i++) b_.s[r4(l.i(i))] = b.block(lb.i(i));
669  }
670 #ifdef USE_BOUND
671  double bound() const { return bound_; }
672  void set_bound(double b) const { bound_ = b; }
673 #endif
674  bi_t &block(int i) { return b_.s[r4(i)]; }
675  const bi_t &block(int i) const { return b_.s[r4(i)]; }
677  unsigned int size(const Range *indices) const {
678  unsigned int r = 1;
679  for (int i=0; i<3; i++) r *= indices[i].block_size(b_.s[r4(i)]);
680  return r;
681  }
684  unsigned int subset_size(const Range *indices,
685  const IndexList &indexlist) const {
686  unsigned int r = 1;
687  for (int i=0; i<indexlist.n(); i++) {
688  int index = indexlist.i(i);
689  r *= indices[index].block_size(b_.s[r4(index)]);
690  }
691  return r;
692  }
697  template <int N2>
698  void assign_blocks(const IndexList &il,
699  const BlockInfo<N2> &bi2, const IndexList &il2) {
700  for (int i=0; i<il.n(); i++) {
701  b_.s[r4(il.i(i))] = bi2.block(il2.i(i));
702  }
703  }
707  template <int N2>
708  void assign_blocks(const IndexList &il,
709  const BlockInfo<N2> &bi2) {
710  for (int i=0; i<il.n(); i++) {
711  b_.s[r4(il.i(i))] = bi2.block(i);
712  }
713  }
717  template <int N2>
718  bool equiv_blocks(const IndexList &il,
719  const BlockInfo<N2> &bi2) {
720  for (int i=0; i<il.n(); i++) {
721  if (b_.s[r4(il.i(i))] != bi2.block(i)) return false;
722  }
723  return true;
724  }
726  void zero() { b_.l = 0; }
727  void print(std::ostream &o=sc::ExEnv::outn()) const {
728  o << "{";
729  for (int i=0; i<3; i++) {
730  if (i!=0) o << " ";
731  o << block(i);
732  }
733  o << "}";
734  }
735  void print_block_sizes(const Range *indices,
736  std::ostream &o=sc::ExEnv::outn()) const {
737  o << "{";
738  for (int i=0; i<3; i++) {
739  if (i) o << " ";
740  o << indices[i].block_size(b_.s[r4(i)]);
741  }
742  o << "}";
743  }
744  void write(sc::StateOut& so) const {
745 #ifdef USE_BOUND
746  so.put(bound_);
747 #endif
748  for (int i=0; i<3; i++) {
749  so.put(int(block(i)));
750  }
751  }
752  void read(sc::StateIn& si) {
753 #ifdef USE_BOUND
754  si.get(bound_);
755 #endif
756  for (int i=0; i<3; i++) {
757  int b;
758  si.get(b);
759  block(i) = b;
760  }
761  }
762  };
763 
764  template<> inline
765  bool IndicesLess<3>::operator() (const BlockInfo<3>&b1,
766  const BlockInfo<3>&b2) const
767  {
768  if (b1.b_.l < b2.b_.l) return true;
769  return false;
770  }
771 
772  template<> inline
773  int IndicesLess<3>::compare(const BlockInfo<3>&b1,
774  const BlockInfo<3>&b2) const
775  {
776  if (b1.b_.l < b2.b_.l) return -1;
777  else if (b1.b_.l > b2.b_.l) return 1;
778  return 0;
779  }
780 #endif // THREE_INDEX_SPECIALIZATIONS
781 
782 #if FOUR_INDEX_SPECIALIZATIONS
783  // Specializations for 4 indices
784  template <>
785  class BlockInfo<4> {
786  private:
787 #ifdef USE_BOUND
788  mutable double bound_;
789 #endif
790  public:
791  union {
792  sc::sc_uint64_t l;
793  sc::sc_uint32_t i[2];
794  sc::sc_uint16_t s[4];
795  sc::sc_uint8_t c[8];
796  } b_;
797  friend class IndicesLess<4>;
798  public:
799  BlockInfo() {
800 #ifdef USE_BOUND
801  bound_ = DBL_MAX;
802 #endif
803  }
804  BlockInfo(const std::vector<bi_t> &v) {
805 #ifdef USE_BOUND
806  bound_ = DBL_MAX;
807 #endif
808  for (int i=0; i<4 && i<v.size(); i++) b_.s[r4(i)] = v[i];
809  }
810  BlockInfo(bi_t b0,bi_t b1,bi_t b2,bi_t b3) {
811 #ifdef USE_BOUND
812  bound_ = DBL_MAX;
813 #endif
814  b_.s[0] = b0;
815  b_.s[1] = b1;
816  b_.s[2] = b2;
817  b_.s[3] = b3;
818  }
819  BlockInfo(const BlockInfo<4> &b, const IndexList &l) {
820 #ifdef USE_BOUND
821  bound_ = DBL_MAX;
822 #endif
823  for (int i=0; i<l.n(); i++) b_.s[r4(i)] = b.b_.s[r4(l.i(i))];
824  }
825  template <int NB>
826  BlockInfo(const IndexList &l,
827  const BlockInfo<NB> &b,
828  const IndexList &lb) {
829 #ifdef USE_BOUND
830  bound_ = DBL_MAX;
831 #endif
832  for (int i=0; i<l.n(); i++) b_.s[r4(l.i(i))] = b.block(lb.i(i));
833  }
834 #ifdef USE_BOUND
835  double bound() const { return bound_; }
836  void set_bound(double b) const { bound_ = b; }
837 #endif
838  bi_t &block(int i) { return b_.s[r4(i)]; }
839  const bi_t &block(int i) const { return b_.s[r4(i)]; }
841  unsigned int size(const Range *indices) const {
842  unsigned int r = 1;
843  for (int i=0; i<4; i++) r *= indices[i].block_size(b_.s[r4(i)]);
844  return r;
845  }
848  unsigned int subset_size(const Range *indices,
849  const IndexList &indexlist) const {
850  unsigned int r = 1;
851  for (int i=0; i<indexlist.n(); i++) {
852  int index = indexlist.i(i);
853  r *= indices[index].block_size(b_.s[r4(index)]);
854  }
855  return r;
856  }
861  template <int N2>
862  void assign_blocks(const IndexList &il,
863  const BlockInfo<N2> &bi2, const IndexList &il2) {
864  for (int i=0; i<il.n(); i++) {
865  b_.s[r4(il.i(i))] = bi2.block(il2.i(i));
866  }
867  }
871  template <int N2>
872  void assign_blocks(const IndexList &il,
873  const BlockInfo<N2> &bi2) {
874  for (int i=0; i<il.n(); i++) {
875  b_.s[r4(il.i(i))] = bi2.block(i);
876  }
877  }
881  template <int N2>
882  bool equiv_blocks(const IndexList &il,
883  const BlockInfo<N2> &bi2) {
884  for (int i=0; i<il.n(); i++) {
885  if (b_.s[r4(il.i(i))] != bi2.block(i)) return false;
886  }
887  return true;
888  }
890  void zero() { b_.l = 0; }
891  void print(std::ostream &o=sc::ExEnv::outn()) const {
892  o << "{";
893  for (int i=0; i<4; i++) {
894  if (i!=0) o << " ";
895  o << block(i);
896  }
897  o << "}";
898  }
899  void print_block_sizes(const Range *indices,
900  std::ostream &o=sc::ExEnv::outn()) const {
901  o << "{";
902  for (int i=0; i<4; i++) {
903  if (i) o << " ";
904  o << indices[i].block_size(b_.s[r4(i)]);
905  }
906  o << "}";
907  }
908  void write(sc::StateOut& so) const {
909 #ifdef USE_BOUND
910  so.put(bound_);
911 #endif
912  for (int i=0; i<4; i++) {
913  so.put(int(block(i)));
914  }
915  }
916  void read(sc::StateIn& si) {
917 #ifdef USE_BOUND
918  si.get(bound_);
919 #endif
920  for (int i=0; i<4; i++) {
921  int b;
922  si.get(b);
923  block(i) = b;
924  }
925  }
926  };
927 
928  template<> inline
929  bool IndicesLess<4>::operator() (const BlockInfo<4>&b1,
930  const BlockInfo<4>&b2) const
931  {
932  if (b1.b_.l < b2.b_.l) return true;
933  return false;
934  }
935 
936  template<> inline
937  int IndicesLess<4>::compare(const BlockInfo<4>&b1,
938  const BlockInfo<4>&b2) const
939  {
940  if (b1.b_.l < b2.b_.l) return -1;
941  else if (b1.b_.l > b2.b_.l) return 1;
942  return 0;
943  }
944 #endif // FOUR_INDEX_SPECIALIZATIONS
945 
946 #ifdef USE_HASH
947  // Unary functor for generating hashes from BlockInfo's
948  template <int N>
949  class BlockInfoHash {
950  public:
951  size_t operator()(const BlockInfo<N> &b) const {
952  size_t r = 0;
953  for (int i=0; i<N; i++) r ^= b.block(i);
954  return r;
955  }
956  };
957 #endif
958 
959 #ifdef USE_HASH
960  // Binary functor for checking if two BlockInfo's are equal
961  template <int N>
962  class BlockInfoEqual {
963  public:
964  bool operator()(const BlockInfo<N> &b1,const BlockInfo<N> &b2) const {
965  for (int i=0; i<N; i++) if (b1.block(i) != b2.block(i)) return false;
966  return true;
967  }
968  };
969 #endif
970 
976  template <int N>
977  class IndexListLess {
978  private:
979  IndexList il_;
980  public:
982  IndexListLess(const IndexList &il): il_(il) {}
985  IndexListLess(const IndexList &il1, const IndexList &il2): il_(il1,il2) {}
986  bool operator() (const BlockInfo<N>&b1, const BlockInfo<N>&b2) const {
987  if (il_.n() == 0) return false;
988  int i0 = il_.i(0);
989  int b1b = b1.block(i0), b2b = b2.block(i0);
990  if (b1b < b2b) return true;
991  if (b1b > b2b) return false;
992  for (int l=1; l<il_.n(); l++) {
993  int i = il_.i(l);
994  b1b = b1.block(i); b2b = b2.block(i);
995  if (b1b < b2b) return true;
996  if (b1b > b2b) return false;
997  }
998 #ifdef USE_BOUND
999  // if indices are equal order according to descending bounds
1000  if (b1.bound() > b2.bound()) return true;
1001 #endif
1002  return false;
1003  }
1004  int compare(const BlockInfo<N>&b1, const BlockInfo<N>&b2) const {
1005  if (il_.n() == 0) return 0;
1006  int i0 = il_.i(0);
1007  int b1b = b1.block(i0), b2b = b2.block(i0);
1008  if (b1b < b2b) return -1;
1009  if (b1b > b2b) return 1;
1010  for (int l=1; l<il_.n(); l++) {
1011  int i = il_.i(l);
1012  b1b = b1.block(i); b2b = b2.block(i);
1013  if (b1b < b2b) return -1;
1014  if (b1b > b2b) return 1;
1015  }
1016 #ifdef USE_BOUND
1017  // if indices are equal order according to descending bounds
1018  if (b1.bound() > b2.bound()) return -1;
1019  else if (b1.bound() < b2.bound()) return 1;
1020 #endif
1021  return 0;
1022  }
1023  };
1024 
1026  class Data: public sc::RefCount {
1027  private:
1028 // double *data_;
1029  size_t ndata_;
1030 
1031  int default_chunksize_;
1032  // maps n_data available to the pointer and n_data used.
1033  typedef std::multimap<size_t, std::pair<double *, size_t> > memmap_t;
1034  memmap_t memmap_;
1035  typedef std::map<double*,memmap_t::iterator> memitermap_t;
1036  memitermap_t memitermap_;
1037  public:
1038  Data();
1039  ~Data();
1040  virtual double *allocate(long size);
1041  virtual void deallocate(double *);
1042  double *data() const;
1043  long ndata() const { return ndata_; }
1044  };
1045 
1047  class Index {
1048  std::string name_;
1049  int value_;
1050  public:
1053  Index(const std::string &name): name_(name), value_(-1) {}
1059  Index(int value): value_(value) {}
1065  Index(const std::string &name, int value): name_(name), value_(value) {}
1068  const std::string name() const { return name_; }
1071  int value() const { return value_; }
1074  bool has_value() const { return value_ >= 0; }
1079  bool symbolically_equivalent(const Index &i) const
1080  { return name_.size() > 0 && name_ == i.name(); }
1082  bool operator == (const Index &i) const
1083  { return name_ == i.name() && value_ == i.value_; }
1085  void set_value(int value) { value_ = value; }
1086  };
1087 
1088  template <int N> class Array;
1089  template <int Nl,int Nr> class ContractProd;
1090  template <int Nl,int Nr> class ContractUnion;
1094  template <int N>
1095  class ContractPart {
1096  Array<N> &array_;
1097  std::vector<Index> indices_;
1098  double factor_;
1099  bool clear_after_use_;
1100  bool skip_bounds_update_;
1101 
1102  template <int N2, class Op>
1103  void do_binary_op_with_fixed_indices(double f, const ContractPart<N2> &o,
1104  bool initarray, Op &op) const;
1105 
1106  template <class Op>
1107  void do_binary_op(double f, const ContractPart<N> &o,
1108  bool initarray, Op &op) const;
1109 
1110  template <int Nl,int Nr>
1111  void doprod(double f, const ContractProd<Nl,Nr> &o,
1112  bool initarray) const;
1113 
1114  template <int Nl,int Nr>
1115  void dounion(const ContractProd<Nl,Nr> &o) const;
1116 
1117  public:
1118  ContractPart(Array<N> &array);
1119  ContractPart(Array<N> &array,
1120  const Index &i1);
1121  ContractPart(Array<N> &array,
1122  const Index &i1,
1123  const Index &i2);
1124  ContractPart(Array<N> &array,
1125  const Index &i1,
1126  const Index &i2,
1127  const Index &i3);
1128  ContractPart(Array<N> &array,
1129  const Index &i1,
1130  const Index &i2,
1131  const Index &i3,
1132  const Index &i4);
1133  ContractPart(Array<N> &array,
1134  const Index &i1,
1135  const Index &i2,
1136  const Index &i3,
1137  const Index &i4,
1138  const Index &i5);
1139  ContractPart(Array<N> &array,
1140  const Index &i1,
1141  const Index &i2,
1142  const Index &i3,
1143  const Index &i4,
1144  const Index &i5,
1145  const Index &i6);
1146  void apply_factor(double f);
1147  double factor() const;
1148  Array<N>& array() const;
1149  const Index &index(int i) const;
1150  bool clear_after_use() const { return clear_after_use_; }
1151 
1152  void operator = (const ContractPart &o) const;
1153  template <int N2>
1154  void operator = (const ContractPart<N2> &o) const;
1155  void operator += (const ContractPart<N> &o) const;
1156  template <int N2>
1157  void operator += (const ContractPart<N2> &o) const;
1158  void operator /= (const ContractPart<N> &o) const;
1159  template <int N2>
1160  void operator /= (const ContractPart<N2> &o) const;
1161  void operator -= (const ContractPart<N> &o) const;
1162  template <int Nl,int Nr>
1163  void operator = (const ContractProd<Nl,Nr> &o) const;
1164  template <int Nl,int Nr>
1165  void operator += (const ContractProd<Nl,Nr> &o) const;
1166  template <int Nl,int Nr>
1167  void operator -= (const ContractProd<Nl,Nr> &o) const;
1168 
1171  template <int N2>
1172  void operator |= (const ContractPart<N2> &o) const;
1173 
1174  template <int Nl,int Nr>
1175  void operator |= (const ContractProd<Nl,Nr> &o) const;
1176 
1180  ContractPart<N> operator ~ () const;
1183  ContractPart<N> skip_bounds_update() const;
1185  double value();
1186  };
1187  template <int N>
1188  ContractPart<N> operator *(double, const ContractPart<N> &);
1189 
1190  template <int Nl, int Nr>
1191  class ContractUnion {
1192  public:
1193  ContractPart<Nl> l;
1194  ContractPart<Nr> r;
1195  ContractUnion(const ContractPart<Nl> &a1, const ContractPart<Nr> &a2):
1196  l(a1), r(a2) {}
1197  };
1198  template <int Nl, int Nr>
1199  ContractUnion<Nl,Nr> operator |(const ContractPart<Nl>&l,
1200  const ContractPart<Nr>&r)
1201  {
1202  return ContractUnion<Nl,Nr>(l,r);
1203  }
1204 
1205  template <int N> class BlockDistrib;
1206 
1207  // This checks to see that all elements of vec, vec[i], are
1208  // in the range [start[i], fence[i]).
1209  // Used with incr (below).
1210  template <class T>
1211  bool
1212  ready(std::vector<T> &vec,
1213  const std::vector<T> &start, const std::vector<T> &fence)
1214  {
1215  if (vec.size() == 0) return false;
1216  for (int i=0; i<vec.size(); i++) {
1217  if (vec[i] < start[i]) return false;
1218  if (vec[i] >= fence[i]) return false;
1219  }
1220  return true;
1221  }
1222 
1223 
1224  // This increments the elements of vec in such a way that all possible
1225  // vec's with vec[i] in the range [start[i], fence[i]). are obtained.
1226  // Used with ready (above).
1227  template <class T>
1228  void
1229  incr(std::vector<T> &vec,
1230  const std::vector<T> &start, const std::vector<T> &fence)
1231  {
1232  if (vec.size() == 0) return;
1233  for (int i=vec.size()-1; i>0; i--) {
1234  if (vec[i] >= fence[i]-1) vec[i] = start[i];
1235  else {
1236  vec[i]++;
1237  return;
1238  }
1239  }
1240  vec[0]++;
1241  }
1242 
1246  template <int N>
1247  class Array {
1248  template <int N1>
1249  friend
1250  void remap(typename Array<N1>::cached_blockmap_t& target,
1251  const Array<N1> &source,
1252  const IndexList &fixed, const BlockInfo<N1> &fixedvals);
1253  template <int N1>
1254  friend
1255  void remap(Array<N1>& target, const Array<N1> &source,
1256  const IndexList &il);
1257  public:
1258  typedef IndicesLess<N> Compare;
1259 #if USE_STL_MULTIMAP
1260  typedef typename std::multimap<BlockInfo<N>, double*, Compare > blockmap_t;
1261  typedef typename std::multimap<BlockInfo<N>, double*, IndexListLess<N> > cached_blockmap_t;
1262 #else
1263  typedef AVLMMap<BlockInfo<N>, double*, Compare > blockmap_t;
1264  typedef AVLMMap<BlockInfo<N>, double*, IndexListLess<N> > cached_blockmap_t;
1265 #endif
1266 #ifdef USE_HASH
1267  typedef typename __gnu_cxx::hash_map<BlockInfo<N>, double*, BlockInfoHash<N>, BlockInfoEqual<N> > blockhash_t;
1268 #endif
1269  private:
1270  sc::auto_vec<Range> indices_;
1271  blockmap_t blocks_;
1272 #ifdef USE_HASH
1273  blockhash_t block_hash_;
1274 #endif
1275  sc::Ref<Data> data_;
1276  bool allocated_;
1277  std::string name_;
1278  int debug_;
1279 #ifdef USE_BOUND
1280  double bound_;
1281 #endif
1282  double tolerance_;
1283 
1284  std::map<IndexList, cached_blockmap_t*> blockmap_cache_;
1285  bool use_blockmap_cache_;
1286 
1287  void init_indices(int n) {
1288  // This routine is used to reset the indices_ to
1289  // avoid a warning about newing a length 0 array
1290  // which occurs when N is used.
1291  if (n) indices_.reset(new Range[n]);
1292  else indices_.reset(0);
1293  }
1294 
1295  void basic_init(const std::string &name, double tol) {
1296  clear();
1297  use_blockmap_cache_ = true;
1298  init_indices(N);
1299  name_ = name;
1300  debug_ = 0;
1301  tolerance_ = tol;
1302  // In the special case that this is a scalar (N == 0),
1303  // add the block and zero it.
1304  if (N == 0) {
1305  add_all_unallocated_blocks();
1306  zero();
1307  }
1308  }
1309  public:
1310  ~Array() {
1311  clear();
1312  }
1313  Array(const Array<N> &a) {
1314  init_indices(N);
1315  debug_ = 0;
1316  operator = (a);
1317  }
1321  void assign_tol(const Array<N> &a,
1322  double tol) {
1323 #ifdef USE_BOUND
1324  bound_ = a.bound_;
1325 #endif
1326  tolerance_ = a.tolerance_;
1327  init_blocks(a, tol);
1328  if (a.allocated_) {
1329  allocate_blocks();
1330  for (typename blockmap_t::const_iterator i = a.blocks_.begin();
1331  i != a.blocks_.end();
1332  i++) {
1333  typename blockmap_t::iterator ithis = blocks_.find(i->first);
1334  if (ithis == blocks_.end()) continue;
1335  memcpy(ithis->second, i->second,
1336  sizeof(*ithis->second)*block_size(ithis->first));
1337  }
1338  allocated_ = true;
1339  }
1340  }
1342  void assign_all(const Array<N> &a) {
1343  assign_tol(a,0.0);
1344  }
1346  Array<N> &operator=(const Array<N> &a) {
1347  assign_tol(a, a.tolerance_);
1348  return *this;
1349  }
1352  Array(const std::string &name = "",
1353  double tol = DBL_EPSILON) { init(name,tol); }
1354  void init(const std::string &name = "",
1355  double tol = DBL_EPSILON) {
1356  basic_init(name,tol);
1357  }
1359  Array(const Range &index,
1360  const std::string & name= "",
1361  double tol = DBL_EPSILON) { init(index,name,tol); }
1362  void init(const Range &index,
1363  const std::string & name= "",
1364  double tol = DBL_EPSILON) {
1365  basic_init(name,tol);
1366  if (N < 1)
1367  throw std::invalid_argument("Array init given too many indices");
1368  for (int i=0; i<N; i++) set_index(i,index);
1369  }
1371  Array(const Range &i0,const Range &i1,
1372  const std::string & name= "",
1373  double tol = DBL_EPSILON) { init(i0,i1,name,tol); }
1374  void init(const Range &i0,const Range &i1,
1375  const std::string & name= "",
1376  double tol = DBL_EPSILON) {
1377  basic_init(name,tol);
1378  if (N < 2)
1379  throw std::invalid_argument("Array init given too many indices");
1380  set_index(0,i0);
1381  for (int i=1; i<N; i++) set_index(i,i1);
1382  }
1384  Array(const Range &i0,const Range &i1,const Range &i2,
1385  const std::string & name= "",
1386  double tol = DBL_EPSILON) { init(i0,i1,i2,name,tol); }
1387  void init(const Range &i0,const Range &i1,const Range &i2,
1388  const std::string & name= "",
1389  double tol = DBL_EPSILON) {
1390  basic_init(name,tol);
1391  if (N < 3)
1392  throw std::invalid_argument("Array init given too many indices");
1393  set_index(0,i0);
1394  set_index(1,i1);
1395  for (int i=2; i<N; i++) set_index(i,i2);
1396  }
1399  Array(const Range &i0,const Range &i1,
1400  const Range &i2,const Range &i3,
1401  const std::string & name= "",
1402  double tol = DBL_EPSILON) { init(i0,i1,i2,i3,name,tol); }
1403  void init(const Range &i0,const Range &i1,
1404  const Range &i2,const Range &i3,
1405  const std::string & name= "",
1406  double tol = DBL_EPSILON) {
1407  basic_init(name,tol);
1408  if (N < 4)
1409  throw std::invalid_argument("Array init given too many indices");
1410  set_index(0,i0);
1411  set_index(1,i1);
1412  set_index(2,i2);
1413  for (int i=3; i<N; i++) set_index(i,i3);
1414  }
1417  Array(const Range &i0,const Range &i1,
1418  const Range &i2,const Range &i3,
1419  const Range &i4,
1420  const std::string & name= "",
1421  double tol = DBL_EPSILON) { init(i0,i1,i2,i3,i4,name,tol); }
1422  void init(const Range &i0,const Range &i1,
1423  const Range &i2,const Range &i3,
1424  const Range &i4,
1425  const std::string & name= "",
1426  double tol = DBL_EPSILON) {
1427  basic_init(name,tol);
1428  if (N < 5)
1429  throw std::invalid_argument("Array init given too many indices");
1430  set_index(0,i0);
1431  set_index(1,i1);
1432  set_index(2,i2);
1433  set_index(3,i3);
1434  for (int i=4; i<N; i++) set_index(i,i4);
1435  }
1438  Array(const Range &i0,const Range &i1,
1439  const Range &i2,const Range &i3,
1440  const Range &i4,const Range &i5,
1441  const std::string & name= "",
1442  double tol = DBL_EPSILON) { init(i0,i1,i2,i3,i4,i5,name,tol); }
1443  void init(const Range &i0,const Range &i1,
1444  const Range &i2,const Range &i3,
1445  const Range &i4,const Range &i5,
1446  const std::string & name= "",
1447  double tol = DBL_EPSILON) {
1448  basic_init(name,tol);
1449  if (N < 6)
1450  throw std::invalid_argument("Array init given too many indices");
1451  set_index(0,i0);
1452  set_index(1,i1);
1453  set_index(2,i2);
1454  set_index(3,i3);
1455  set_index(4,i4);
1456  for (int i=5; i<N; i++) set_index(i,i5);
1457  }
1460  Array(const Range *ranges,
1461  double tol = DBL_EPSILON): data_(new Data), debug_(0),
1462  tolerance_(tol) {
1463 #ifdef USE_BOUND
1464  bound_ = DBL_MAX;
1465 #endif
1466  set_indices(ranges);
1467  }
1468 #ifdef USE_BOUND
1469  double bound() const { return bound_; }
1470  void set_bound(double b) { bound_ = b; }
1471 #endif
1472  double tolerance() const { return tolerance_; }
1473  void set_tolerance(double t) { tolerance_ = t; }
1475  void clear() {
1476  clear_blockmap_cache();
1477  allocated_ = false;
1478  blocks_.clear();
1479 #ifdef USE_HASH
1480  block_hash_.clear();
1481 #endif
1482 #ifdef USE_BOUND
1483  bound_ = DBL_MAX;
1484 #endif
1485  data_ = new Data;
1486  }
1488  static int nindex() { return N; }
1490  const blockmap_t &blockmap() const { return blocks_; }
1491 #ifdef USE_HASH
1492  const blockhash_t &blockhash() const { return block_hash_; }
1494 #endif
1495 
1496  void set_index(int i, const Range &idx) { indices_[i] = idx; }
1499  void set_indices(const Range *r) {
1500  for (int i=0; i<N; i++) set_index(i,r[i]);
1501  }
1503  const Range &index(int i) const { return indices_.get()[i]; }
1505  const Range *indices() const { return indices_.get(); }
1507  int block_size(const BlockInfo<N> &bi) const {
1508  int bs = 1;
1509  for (int i=0; i<N; i++) bs *= indices_.get()[i].block_size(bi.block(i));
1510  return bs;
1511  }
1515  typename blockmap_t::value_type &
1516  add_unallocated_block(const BlockInfo<N>&b) {
1517  if (allocated_) {
1518  throw std::runtime_error(
1519  "cannot add unalloced block to alloced array");
1520  }
1521  // check the indices
1522  for (int i=0; i<N; i++) {
1523  if (b.block(i) >= index(i).nblock()) {
1524  throw std::runtime_error(
1525  "attempted to add a block with invalid block indices"
1526  );
1527  }
1528  }
1529 #if USE_STL_MULTIMAP
1530  typename blockmap_t::iterator bi = blocks_.find(b);
1531  if (bi != blocks_.end()) return *bi;
1532  return *blocks_.insert(typename blockmap_t::value_type(b,(double*)0));
1533 #else
1534  return *blocks_.insert_unique(typename blockmap_t::value_type(b,(double*)0));
1535 #endif
1536  }
1540  typename blockmap_t::value_type &
1541  add_allocated_block(const BlockInfo<N>&b) {
1542  allocated_ = 1;
1543  // check the indices
1544  for (int i=0; i<N; i++) {
1545  if (b.block(i) >= index(i).nblock()) {
1546  throw std::runtime_error(
1547  "attempted to add a block with invalid block indices"
1548  );
1549  }
1550  }
1551  typename blockmap_t::iterator bi = blocks_.find(b);
1552  if (bi != blocks_.end()) return *bi;
1553 
1554  double *data = data_->allocate(block_size(b));
1555  return *blocks_.insert(typename blockmap_t::value_type(b,data));
1556  }
1560  typename blockmap_t::value_type &
1561  add_zeroed_block(const BlockInfo<N>&b) {
1562  allocated_ = 1;
1563  // check the indices
1564  for (int i=0; i<N; i++) {
1565  if (b.block(i) >= index(i).nblock()) {
1566  throw std::runtime_error(
1567  "attempted to add a block with invalid block indices"
1568  );
1569  }
1570  }
1571 
1572  typename blockmap_t::iterator bi = blocks_.find(b);
1573  if (bi != blocks_.end()) return *bi;
1574 
1575  int size = block_size(b);
1576  double *data = data_->allocate(size);
1577  memset(data, 0, sizeof(double)*size);
1578 
1579  return *blocks_.insert(typename blockmap_t::value_type(b,data));
1580  }
1584  void remove_block(const BlockInfo<N>&b) {
1585  typename blockmap_t::iterator bi = blocks_.find(b);
1586  if (bi != blocks_.end()) {
1587  data_->deallocate(bi->second);
1588  blocks_.erase(bi);
1589  }
1590  }
1591  void relocate_block(const BlockInfo<N>&b_old,const BlockInfo<N>&b_new) {
1592  typename blockmap_t::iterator bi = blocks_.find(b_old);
1593  if (bi != blocks_.end()) {
1594  double *data = bi->second;
1595  blocks_.erase(bi);
1596  blocks_.insert(std::make_pair(b_new,data));
1597  }
1598  }
1604  typename blockmap_t::iterator
1605  add_new_unallocated_block(const typename blockmap_t::iterator &hint,
1606  const BlockInfo<N>&b) {
1607  if (allocated_) {
1608  throw std::runtime_error(
1609  "cannot add unalloced block to alloced array");
1610  }
1611  // check the indices
1612  for (int i=0; i<N; i++) {
1613  if (b.block(i) >= index(i).nblock()) {
1614  throw std::runtime_error(
1615  "attempted to add a block with invalid block indices"
1616  );
1617  }
1618  }
1619 #if USE_STL_MULTIMAP
1620  return blocks_.insert(blockmap_t::value_type(b,(double*)0));
1621 #else
1622  return blocks_.insert_new(hint, typename blockmap_t::value_type(b,(double*)0));
1623 #endif
1624  }
1629  typename blockmap_t::iterator
1630  add_new_unallocated_block(const BlockInfo<N>&b) {
1631  if (allocated_) {
1632  throw std::runtime_error(
1633  "cannot add unalloced block to alloced array");
1634  }
1635  // check the indices
1636  for (int i=0; i<N; i++) {
1637  if (b.block(i) > index(i).nblock()) {
1638  throw std::runtime_error(
1639  "attempted to add a block with invalid block indices"
1640  );
1641  }
1642  }
1643 #if USE_STL_MULTIMAP
1644  return blocks_.insert(blockmap_t::value_type(b,(double*)0));
1645 #else
1646  return blocks_.insert_new(typename blockmap_t::value_type(b,(double*)0));
1647 #endif
1648  }
1652  void
1653  add_all_unallocated_blocks() {
1654 
1655  if (allocated_) {
1656  throw std::runtime_error(
1657  "cannot add unalloced block to alloced array");
1658  }
1659 
1660  clear();
1661 
1662  if (N == 0) {
1663  BlockInfo<N> bi;
1664  add_new_unallocated_block(bi);
1665  return;
1666  }
1667 
1668  std::vector<sma2::bi_t> b_indices(N);
1669  std::vector<sma2::bi_t> b_start(N);
1670  std::vector<sma2::bi_t> b_fence(N);
1671 
1672  std::fill(b_start.begin(), b_start.end(), 0);
1673  for (int i=0; i<N; i++) b_fence[i] = index(i).nblock();
1674 
1675  for (b_indices=b_start;
1676  ready(b_indices, b_start, b_fence);
1677  incr(b_indices, b_start, b_fence)) {
1678  sma2::BlockInfo<N> bi(b_indices);
1679  add_new_unallocated_block(bi);
1680  }
1681  }
1684  typename blockmap_t::iterator initial_hint() {
1685  return blocks_.end();
1686  }
1688  int n_block() const {
1689  return blocks_.size();
1690  }
1692  size_t n_element() const {
1693  size_t n = 0;
1694  for (typename blockmap_t::const_iterator i = blocks_.begin();
1695  i != blocks_.end();
1696  i++) {
1697  n += block_size(i->first);
1698  }
1699  return n;
1700  }
1702  size_t n_element_allocated() const {
1703  if (data_.null()) return 0;
1704  return data_->ndata();
1705  }
1706  std::string name() const { return name_; }
1707  double max_n_block() {
1708  double n_block_max = 1;
1709  for (int i=0; i<N; i++) {
1710  n_block_max *= index(i).nblock();
1711  }
1712  return n_block_max;
1713  }
1714  double max_n_element() {
1715  double n_element_max = 1;
1716  for (int i=0; i<N; i++) {
1717  n_element_max *= index(i).nindex();
1718  }
1719  return n_element_max;
1720  }
1722  void assign(double val) {
1723  allocate_blocks();
1724  for (typename blockmap_t::iterator i = blocks_.begin();
1725  i != blocks_.end();
1726  i++) {
1727  int n = block_size(i->first);
1728  double *d = i->second;
1729  for (int j=0; j<n; j++) d[j] = val;
1730 #ifdef USE_BOUND
1731  i->first.set_bound(val);
1732 #endif
1733  }
1734 #ifdef USE_BOUND
1735  bound_ = val;
1736 #endif
1737  }
1739  void zero() {
1740  assign(0.0);
1741  }
1743  void compute_bounds() {
1744  allocate_blocks();
1745 #ifdef USE_BOUND
1746  bound_ = 0.0;
1747  for (typename blockmap_t::iterator i = blocks_.begin();
1748  i != blocks_.end();
1749  i++) {
1750  int n = block_size(i->first);
1751  double *d = i->second;
1752  double maxval = 0.0;
1753  for (int j=0; j<n; j++) {
1754  double val = fabs(d[j]);
1755  if (val > maxval) maxval = val;
1756  }
1757  i->first.set_bound(maxval);
1758  if (maxval > bound_) bound_ = maxval;
1759  }
1760 #endif
1761  }
1764  void allocate_blocks() {
1765  if (allocated_) return;
1766  size_t nblock = 0;
1767  long n = n_element();
1768  // clj debug
1769  //if (n > 1000000) {
1770  // sc::Ref<sc::MessageGrp> msg
1771  // = sc::MessageGrp::get_default_messagegrp();
1772  // int me = msg->me();
1773  // sc::ExEnv::outn() << me << ": " << name()
1774  // << " allocating " << n << std::endl;
1775  // }
1776  // clj end debug
1777  double *dat = data_->allocate(n);
1778  for (typename blockmap_t::iterator i = blocks_.begin();
1779  i != blocks_.end();
1780  i++,nblock++) {
1781  if (i->second == 0) {
1782  int n = block_size(i->first);
1783  i->second = dat;
1784  dat += n;
1785  }
1786  else {
1787  throw std::runtime_error("allocate_blocks: duplicate alloc");
1788  }
1789  }
1790 #ifdef USE_HASH
1791  block_hash_.resize(nblock);
1792  for (typename blockmap_t::iterator i = blocks_.begin();
1793  i != blocks_.end();
1794  i++) {
1795  block_hash_.insert(*i);
1796  }
1797 #endif
1798  allocated_ = true;
1799  }
1801  void deallocate_blocks() {
1802  data_ = new Data;
1803  for (typename blockmap_t::iterator i = blocks_.begin();
1804  i != blocks_.end();
1805  i++) {
1806  i->second = 0;
1807  }
1808  allocated_ = false;
1809 #ifdef USE_BOUND
1810  bound_ = DBL_MAX;
1811 #endif
1812  }
1813  ContractPart<N> operator()() {
1814  return ContractPart<N>(*this);
1815  }
1816  ContractPart<N> operator()(const Index &i1) {
1817  return ContractPart<N>(*this,i1);
1818  }
1819  ContractPart<N> operator()(const Index &i1, const Index &i2) {
1820  return ContractPart<N>(*this,i1,i2);
1821  }
1822  ContractPart<N> operator()(const std::string &i1,
1823  const std::string &i2) {
1824  return ContractPart<N>(*this,i1,i2);
1825  }
1826  ContractPart<N> operator()(const Index &i1, const Index &i2, const Index &i3) {
1827  return ContractPart<N>(*this,i1,i2,i3);
1828  }
1829  ContractPart<N> operator()(const std::string &i1,
1830  const std::string &i2,
1831  const std::string &i3) {
1832  return ContractPart<N>(*this,i1,i2,i3);
1833  }
1834  ContractPart<N> operator()(const Index &i1, const Index &i2,
1835  const Index &i3, const Index &i4) {
1836  return ContractPart<N>(*this,i1,i2,i3,i4);
1837  }
1838  ContractPart<N> operator()(const std::string &i1,
1839  const std::string &i2,
1840  const std::string &i3,
1841  const std::string &i4) {
1842  return ContractPart<N>(*this,i1,i2,i3,i4);
1843  }
1844  ContractPart<N> operator()(const Index &i1, const Index &i2,
1845  const Index &i3, const Index &i4,
1846  const Index &i5) {
1847  return ContractPart<N>(*this,i1,i2,i3,i4,i5);
1848  }
1849  ContractPart<N> operator()(const std::string &i1,
1850  const std::string &i2,
1851  const std::string &i3,
1852  const std::string &i4,
1853  const std::string &i5) {
1854  return ContractPart<N>(*this,i1,i2,i3,i4,i5);
1855  }
1856  ContractPart<N> operator()(const Index &i1, const Index &i2,
1857  const Index &i3, const Index &i4,
1858  const Index &i5, const Index &i6) {
1859  return ContractPart<N>(*this,i1,i2,i3,i4,i5,i6);
1860  }
1861  ContractPart<N> operator()(const std::string &i1,
1862  const std::string &i2,
1863  const std::string &i3,
1864  const std::string &i4,
1865  const std::string &i5,
1866  const std::string &i6) {
1867  return ContractPart<N>(*this,i1,i2,i3,i4,i5,i6);
1868  }
1870  void operator *= (double f) {
1871  if (f == 0.0) { clear(); compute_bounds(); return; }
1872  for (typename blockmap_t::const_iterator i = blocks_.begin();
1873  i != blocks_.end();
1874  i++) {
1875  double fabs_f = fabs(f);
1876  int n = block_size(i->first);
1877  double *data = i->second;
1878  for (int j=0; j<n; j++) data[j] *= f;
1879 #ifdef USE_BOUND
1880  i->first.set_bound(fabs_f * i->first.bound());
1881 #endif
1882  }
1883  }
1885  void init_blocks(const Array<N> &a, double tol) {
1886  clear();
1887  set_indices(a.indices());
1888  for (typename blockmap_t::const_iterator i = a.blocks_.begin();
1889  i != a.blocks_.end();
1890  i++) {
1891  if (a.block_max_abs(i) > tol) {
1892  add_unallocated_block(i->first);
1893  }
1894  }
1895  }
1897  void init_blocks(const Array<N> &a) {
1898  init_blocks(a, tolerance());
1899  }
1901  double block_max_abs(typename blockmap_t::const_iterator &b) const {
1902  int s = block_size(b->first);
1903  const double *dat = b->second;
1904  double maxabs = 0.0;
1905  for (int i=0; i<s; i++) {
1906  double v = fabs(dat[i]);
1907  if (v > maxabs) maxabs = v;
1908  }
1909  return maxabs;
1910  }
1912  double max_abs_element() {
1913  double max_abs = 0.0;
1914  for (typename blockmap_t::const_iterator i = blocks_.begin();
1915  i != blocks_.end();
1916  i++) {
1917  max_abs = std::max(max_abs, this->block_max_abs(i));
1918  }
1919  return max_abs;
1920  }
1921 
1923  void print_local(std::ostream&o=sc::ExEnv::outn()) const;
1924  void print(const sc::Ref<sc::MessageGrp> &grp = 0,
1925  bool distributed = false, std::ostream&o=sc::ExEnv::outn()) const;
1926 
1928  int &debug() { return debug_; }
1929  const int &debug() const { return debug_; }
1930 
1931  void read(sc::StateIn&si) {
1932  clear();
1933  for (int i=0; i<N; i++) {
1934  indices_[i].read(si);
1935  indices_[i].print();
1936  }
1937  int allocatedval;
1938  si.get(allocatedval);
1939  char *str;
1940  si.getstring(str);
1941  name_ = str;
1942  delete[] str;
1943 #ifdef USE_BOUND
1944  si.get(bound_);
1945  sc::ExEnv::outn() << "bound = " << bound_ << std::endl;
1946 #endif
1947  si.get(tolerance_);
1948  sc::ExEnv::outn() << "tolerance = " << tolerance_ << std::endl;
1949  int nblock;
1950  si.get(nblock);
1951  BlockInfo<N> bi;
1952  for (int i=0; i<nblock; i++) {
1953  bi.read(si);
1954  add_unallocated_block(bi);
1955  }
1956  if (allocatedval) {
1957  allocate_blocks();
1958  si.get_array_double(data_->data(), data_->ndata());
1959  }
1960  }
1961  void write(sc::StateOut&so) const {
1962  for (int i=0; i<N; i++) indices_[i].write(so);
1963  int bval = allocated_;
1964  so.put(bval);
1965  so.putstring(name_.c_str());
1966 #ifdef USE_BOUND
1967  so.put(bound_);
1968 #endif
1969  so.put(tolerance_);
1970  int blocks_size = int(blocks_.size());
1971  if (blocks_.size() != blocks_size) {
1972  throw std::runtime_error("Array::write: blocks truncated");
1973  }
1974  so.put(blocks_size);
1975  for (typename blockmap_t::const_iterator i = blocks_.begin();
1976  i != blocks_.end();
1977  i++) {
1978  i->first.write(so);
1979  }
1980  if (allocated_) so.put_array_double(data_->data(), data_->ndata());
1981  }
1984  void parallel_accumulate(const sc::Ref<sc::MessageGrp> &grp);
1988  void parallel_union(const sc::Ref<sc::MessageGrp> &grp);
1991  void replicated_from_distributed(const sc::Ref<sc::MessageGrp>&,
1992  const Array<N> &);
1995  void distributed_from_distributed(const sc::Ref<sc::MessageGrp>&,
1996  const BlockDistrib<N> &,
1997  Array<N> &,
1998  bool clear_source_array = false,
1999  bool ignore_block_distrib_throws = false);
2000 
2002  double value() {
2003  if (N != 0) throw sc::ProgrammingError("Array:value: N!=0",
2004  __FILE__, __LINE__);
2005  return blockmap().begin()->second[0];
2006  }
2008  void clear_blockmap_cache() {
2009  for (typename std::map<IndexList,cached_blockmap_t*>::iterator
2010  iter = blockmap_cache_.begin();
2011  iter != blockmap_cache_.end();
2012  iter++) {
2013  delete iter->second;
2014  }
2015  blockmap_cache_.clear();
2016  }
2022  void set_use_blockmap_cache(bool ubmc) {
2023  use_blockmap_cache_ = ubmc;
2024  }
2026  bool use_blockmap_cache() const {
2027  return use_blockmap_cache_;
2028  }
2030  bool blockmap_cache_entry_exists(const IndexList &il) {
2031  return blockmap_cache_.find(il) != blockmap_cache_.end();
2032  }
2034  cached_blockmap_t &blockmap_cache_entry(const IndexList &il) {
2035  cached_blockmap_t *&entry = blockmap_cache_[il];
2036  if (entry == 0) {
2037  entry = new cached_blockmap_t(il);
2038  IndexList nullil;
2039  BlockInfo<N> nullbi;
2040  remap(*entry,*this,nullil,nullbi);
2041  }
2042  return *entry;
2043  }
2044  };
2045  template <int N>
2046  inline std::ostream& operator << (std::ostream&o, const Array<N> &a)
2047  {
2048  a.print(o);
2049  return o;
2050  }
2051 
2052  inline int offset2(int i, int ni, int j, int nj) {
2053  return i*nj+j;
2054  }
2055 
2056  template <int N> double scalar_contract(Array<N> &c, Array<N> &a, const IndexList &alist);
2057 
2061  template <int Nl, int Nr>
2062  class ContractProd {
2063  public:
2064  sc::Ref<sc::RegionTimer> regtimer;
2065  ContractPart<Nl> l;
2066  ContractPart<Nr> r;
2067  void apply_factor(double f) { l.apply_factor(f); }
2068  ContractProd(const ContractPart<Nl> &a1, const ContractPart<Nr> &a2):
2069  l(a1), r(a2) {}
2070  operator double() {
2071  std::vector<int> ivec;
2072  for (int i=0; i<Nr; i++) {
2073  for (int j=0; j<Nl; j++) {
2074  if (r.index(i) == l.index(j)) {
2075  ivec.push_back(j);
2076  break;
2077  }
2078  }
2079  }
2080  IndexList il(ivec);
2081  return l.factor() * r.factor()
2082  * sma2::scalar_contract(l.array(), r.array(), il);
2083  }
2084  ContractProd timer(const sc::Ref<sc::RegionTimer> &t) {
2085  ContractProd<Nl,Nr> ret(*this);
2086  ret.regtimer = t;
2087  return ret;
2088  }
2089  };
2090  template <int Nl, int Nr>
2091  ContractProd<Nl,Nr> operator *(const ContractPart<Nl>&l,
2092  const ContractPart<Nr>&r)
2093  {
2094  return ContractProd<Nl,Nr>(l,r);
2095  }
2096  template <int Nl, int Nr>
2097  ContractProd<Nl,Nr> operator *(double f, const ContractProd<Nl,Nr> &prod)
2098  {
2099  ContractProd<Nl,Nr> result(prod);
2100  result.apply_factor(f);
2101  return result;
2102  }
2103 
2108  template <int N>
2109  void
2110  remap(Array<N>& target, const Array<N> &source,
2111  const IndexList &il)
2112  {
2113  for (int i=0; i<N; i++) target.set_index(i,source.index(i));
2114  const typename Array<N>::blockmap_t &sourceblocks
2115  = source.blockmap();
2116  target.blocks_.clear();
2117  for (typename Array<N>::blockmap_t::const_iterator i = sourceblocks.begin();
2118  i != sourceblocks.end();
2119  i++) {
2120  const BlockInfo<N> &orig_bi(i->first);
2121  BlockInfo<N> new_bi(orig_bi,il);
2122  target.blocks_.insert(std::pair<BlockInfo<N>,double*>(new_bi,
2123  i->second));
2124  }
2125  target.data_ = source.data_;
2126  }
2127 
2133  template <int N>
2134  void
2135  remap(typename Array<N>::cached_blockmap_t& target,
2136  const Array<N> &source,
2137  const IndexList &fixed, const BlockInfo<N> &fixedvals)
2138  {
2139  if (!fixed.is_identity_permutation())
2140  throw std::invalid_argument("remap requires fixed indices first");
2141  const typename Array<N>::blockmap_t &sourceblocks
2142  = source.blockmap();
2143  target.clear();
2144  typename Array<N>::blockmap_t::const_iterator begin, end;
2145  if (fixed.n() == 0) {
2146  begin = sourceblocks.begin();
2147  end = sourceblocks.end();
2148  }
2149  else {
2150  BlockInfo<N> sbi;
2151 
2152  sbi.zero();
2153  sbi.assign_blocks(fixed, fixedvals);
2154 #ifdef USE_BOUND
2155  sbi.set_bound(DBL_MAX);
2156 #endif
2157  begin = sourceblocks.lower_bound(sbi);
2158 
2159  for (int i=0; i<N; i++) sbi.block(i) = source.index(i).nblock();
2160  sbi.assign_blocks(fixed, fixedvals);
2161 #ifdef USE_BOUND
2162  sbi.set_bound(0.0);
2163 #endif
2164  end = sourceblocks.upper_bound(sbi);
2165  }
2166  for (typename Array<N>::blockmap_t::const_iterator i = begin;
2167  i != end;
2168  i++) {
2169  target.insert(*i);
2170  }
2171  }
2172 
2174  void extract_diagonal(Array<2> &array, std::vector<double> &diag);
2175 
2177  void scale_diagonal(Array<2> &array, double factor);
2178 
2181  template <int N, class V>
2182  void
2183  apply_denominator(Array<N> &array, double denominator_offset,
2184  const std::vector<V> &eigenvalues);
2185 
2188  template <int N>
2189  class BlockIter {
2190  int bs_[N];
2191  int i_[N];
2192  bool ready_;
2193  void init(const Range *index,
2194  const BlockInfo<N> &bi,
2195  const IndexList &il) {
2196  for (int i=0; i<N; i++) {
2197  bs_[i] = index[il.i(i)].block_size(bi.block(il.i(i)));
2198  }
2199  }
2200  public:
2201  BlockIter(const Range *index,
2202  const BlockInfo<N> &bi) {
2203  init(index,bi,IndexList::identity(N));
2204  }
2205  BlockIter(const Range *index,
2206  const BlockInfo<N> &bi,
2207  const IndexList &il) {
2208  init(index,bi,il);
2209  }
2211  void start() { ready_ = true; for (int i=0; i<N; i++) i_[i] = 0; }
2213  bool ready() const { return ready_; }
2215  void operator++(int) {
2216  for (int i=N-1; i>=0; i--) {
2217  i_[i]++;
2218  if (i_[i] == bs_[i]) {
2219  i_[i] = 0;
2220  }
2221  else return;
2222  }
2223  ready_ = false;
2224  }
2225  int offset() {
2226  int r = 0;
2227  for (int i=0; i<N; i++) {
2228  r = r * bs_[i] + i_[i];
2229  }
2230  return r;
2231  }
2232  int subset_offset(const IndexList &subset_indexlist) {
2233  int r = 0;
2234  for (int i=0; i<subset_indexlist.n(); i++) {
2235  r = r * bs_[subset_indexlist.i(i)]
2236  + i_[subset_indexlist.i(i)];
2237  }
2238  return r;
2239  }
2241  int block_size(int i) const { return bs_[i]; }
2242  int &index(int i) { return i_[i]; }
2243  const int &index(int i) const { return i_[i]; }
2244  };
2245 
2247  template <>
2248  class BlockIter<0> {
2249  bool ready_;
2250  public:
2251  BlockIter(const Range *index,
2252  const BlockInfo<0> &bi) {}
2253  BlockIter(const Range *index,
2254  const BlockInfo<0> &bi,
2255  const IndexList &il) {}
2256  void start() { ready_ = true; }
2257  bool ready() const { return ready_; }
2258  void operator++(int) {
2259  ready_ = false;
2260  }
2261  int offset() {
2262  return 0;
2263  }
2264  int subset_offset(const IndexList &subset_indexlist) {
2265  return 0;
2266  }
2267  int block_size(int i) const {
2268  throw sc::ProgrammingError("BlockIter<0>:block_size: can not be called");
2269  }
2270  int &index(int i) {
2271  throw sc::ProgrammingError("BlockIter<0>:index: can not be called");
2272  }
2273  const int &index(int i) const {
2274  throw sc::ProgrammingError("BlockIter<0>:index: can not be called");
2275  }
2276  };
2277 
2279  template <int N, class Op>
2280  void binary_op(Array<N> &c,
2281  double alpha, const Array<N> &a, const IndexList &alist,
2282  bool skip_bounds_update, Op &op)
2283  {
2284  // Consistency checks
2285  if (alist.n() != N) {
2286  throw std::invalid_argument("sma2::binary_op: # of indices inconsistent");
2287  }
2288 // sc::ExEnv::outn() << "a's indices on c" << std::endl;
2289 // for (int i=0; i<N; i++) {
2290 // sc::ExEnv::outn() << " " << alist.i(i);
2291 // }
2292 // sc::ExEnv::outn() << std::endl;
2293  for (int i=0; i<N; i++) {
2294  if (c.index(alist.i(i)) != a.index(i))
2295  throw std::invalid_argument("sma2::binary_op: indices don't agree");
2296  }
2297 
2298  bool same_index_order = alist.is_identity();
2299 
2300  // if c has fewer blocks then it drives the summation
2301  // otherwise a drives
2302  if (c.n_block() < a.n_block()) {
2303  const typename Array<N>::blockmap_t &amap = a.blockmap();
2304  const typename Array<N>::blockmap_t &cmap = c.blockmap();
2305  IndexList clist = alist.reverse_mapping();
2306  for (typename Array<N>::blockmap_t::const_iterator citer = cmap.begin();
2307  citer != cmap.end();
2308  citer++) {
2309  const BlockInfo<N> &cbi = citer->first;
2310  BlockInfo<N> abi(cbi,alist);
2311  //sc::ExEnv::outn() << "summing into " << cbi << std::endl;
2312  typename Array<N>::blockmap_t::const_iterator ablock
2313  = amap.find(abi);
2314  if (ablock == amap.end()) continue;
2315  double *cdata = citer->second;
2316  double *adata = ablock->second;
2317  if (same_index_order) {
2318  int sz = c.block_size(cbi);
2319  for (int i=0; i<sz; i++) op(cdata[i],alpha * adata[i]);
2320  }
2321  else {
2322  BlockIter<N> cbiter(c.indices(),cbi);
2323  int coff = 0;
2324  for (cbiter.start(); cbiter.ready(); cbiter++,coff++) {
2325  op(cdata[coff],
2326  alpha * adata[cbiter.subset_offset(alist)]);
2327  }
2328  }
2329  }
2330  }
2331  else {
2332  const typename Array<N>::blockmap_t &amap = a.blockmap();
2333  const typename Array<N>::blockmap_t &cmap = c.blockmap();
2334  IndexList clist = alist.reverse_mapping();
2335  for (typename Array<N>::blockmap_t::const_iterator aiter = amap.begin();
2336  aiter != amap.end();
2337  aiter++) {
2338  const BlockInfo<N> &abi = aiter->first;
2339  BlockInfo<N> cbi(abi,clist);
2340  typename Array<N>::blockmap_t::const_iterator cblock
2341  = cmap.find(cbi);
2342  if (cblock == cmap.end()) continue;
2343  double *adata = aiter->second;
2344  double *cdata = cblock->second;
2345  if (same_index_order) {
2346  int sz = a.block_size(abi);
2347  for (int i=0; i<sz; i++) op(cdata[i], alpha * adata[i]);
2348  }
2349  else {
2350  BlockIter<N> cbiter(c.indices(),cbi);
2351  int coff = 0;
2352  for (cbiter.start(); cbiter.ready(); cbiter++,coff++) {
2353  op(cdata[coff],
2354  alpha * adata[cbiter.subset_offset(alist)]);
2355  }
2356  }
2357  }
2358  }
2359  if (!skip_bounds_update) c.compute_bounds();
2360  }
2361 
2365  template <class T>
2366  void pack_matrix_into_array(const T &m, Array<2> &smaa) {
2367  smaa.allocate_blocks();
2368  const Array<2>::blockmap_t &Abm = smaa.blockmap();
2369  for (Array<2>::blockmap_t::const_iterator i = Abm.begin();
2370  i != Abm.end();
2371  i++) {
2372  const BlockInfo<2> &Abi = i->first;
2373  int off0 = smaa.index(0).block_offset(Abi.block(0));
2374  int off1 = smaa.index(1).block_offset(Abi.block(1));
2375  BlockIter<2> j(smaa.indices(), Abi);
2376  double *dat = i->second;
2377  for (j.start(); j.ready(); j++) {
2378  dat[j.offset()] = m(off0+j.index(0), off1+j.index(1));
2379  }
2380  }
2381  smaa.compute_bounds();
2382  }
2383 
2387  template <class T>
2388  void pack_vector_into_array(const T &m, Array<1> &smaa) {
2389  smaa.allocate_blocks();
2390  const Array<1>::blockmap_t &Abm = smaa.blockmap();
2391  for (Array<1>::blockmap_t::const_iterator i = Abm.begin();
2392  i != Abm.end();
2393  i++) {
2394  const BlockInfo<1> &Abi = i->first;
2395  int off0 = smaa.index(0).block_offset(Abi.block(0));
2396  BlockIter<1> j(smaa.indices(), Abi);
2397  double *dat = i->second;
2398  for (j.start(); j.ready(); j++) {
2399  dat[j.offset()] = m(off0+j.index(0));
2400  }
2401  }
2402  smaa.compute_bounds();
2403  }
2404 
2412  template <class T>
2413  void pack_submatrix_into_array(const T &m, Array<2> &smaa,
2414  const std::vector<int> &index_map_i,
2415  const std::vector<int> &index_map_j) {
2416 
2417  std::vector<int> reverse_map_i(smaa.index(0).nindex());
2418  std::vector<int> reverse_map_j(smaa.index(1).nindex());
2419  std::fill(reverse_map_i.begin(), reverse_map_i.end(), -1);
2420  std::fill(reverse_map_j.begin(), reverse_map_j.end(), -1);
2421  for (int i=0; i<index_map_i.size(); i++) {
2422  reverse_map_i[index_map_i[i]] = i;
2423  }
2424  for (int i=0; i<index_map_j.size(); i++) {
2425  reverse_map_j[index_map_j[i]] = i;
2426  }
2427 
2428  smaa.allocate_blocks();
2429  const Array<2>::blockmap_t &Abm = smaa.blockmap();
2430  for (Array<2>::blockmap_t::const_iterator i = Abm.begin();
2431  i != Abm.end();
2432  i++) {
2433  const BlockInfo<2> &Abi = i->first;
2434  int off0 = smaa.index(0).block_offset(Abi.block(0));
2435  int off1 = smaa.index(1).block_offset(Abi.block(1));
2436  BlockIter<2> j(smaa.indices(), Abi);
2437  double *dat = i->second;
2438  for (j.start(); j.ready(); j++) {
2439  int m_i = reverse_map_i[off0+j.index(0)];
2440  int m_j = reverse_map_j[off1+j.index(1)];
2441  if (m_i == -1 || m_j == -1) continue;
2442  dat[j.offset()] = m(m_i, m_j);
2443  }
2444  }
2445  smaa.compute_bounds();
2446  }
2447 
2455  template <class T>
2456  void pack_subvector_into_array(const T &m, Array<1> &smaa,
2457  const std::vector<int> &index_map) {
2458 
2459  std::vector<int> reverse_map(smaa.index(0).nindex());
2460  std::fill(reverse_map.begin(), reverse_map.end(), -1);
2461  for (int i=0; i<index_map.size(); i++) {
2462  reverse_map[index_map[i]] = i;
2463  }
2464 
2465  smaa.allocate_blocks();
2466  const Array<1>::blockmap_t &Abm = smaa.blockmap();
2467  for (Array<1>::blockmap_t::const_iterator i = Abm.begin();
2468  i != Abm.end();
2469  i++) {
2470  const BlockInfo<1> &Abi = i->first;
2471  int off0 = smaa.index(0).block_offset(Abi.block(0));
2472  BlockIter<1> j(smaa.indices(), Abi);
2473  double *dat = i->second;
2474  for (j.start(); j.ready(); j++) {
2475  int m_i = reverse_map[off0+j.index(0)];
2476  if (m_i == -1) continue;
2477  dat[j.offset()] = m(m_i);
2478  }
2479  }
2480  smaa.compute_bounds();
2481  }
2482 
2487  template <class T>
2488  void pack_dense_matrix_into_array(const T &m, Array<2> &smaa,
2489  const std::set<int> &row_blocks,
2490  const std::set<int> &col_blocks) {
2491  // allocate the nonzero blocks of smaa
2492  smaa.clear();
2493  for (std::set<int>::const_iterator it1 = row_blocks.begin();
2494  it1 != row_blocks.end(); it1++) {
2495  for (std::set<int>::const_iterator it2 = col_blocks.begin();
2496  it2 != col_blocks.end(); it2++) {
2497  BlockInfo<2> bi;
2498  bi.block(0) = *it1; bi.block(1) = *it2;
2499  smaa.add_unallocated_block(bi);
2500  }
2501  }
2502  smaa.allocate_blocks();
2503 
2504  // initialize the offset arrays
2505  int off=0;
2506  std::vector<int> row_offset(smaa.index(0).nblock());
2507  for (std::set<int>::const_iterator i = row_blocks.begin();
2508  i!=row_blocks.end(); i++) {
2509  row_offset[*i] = off;
2510  off+=smaa.index(0).block_size(*i);
2511  }
2512  off=0;
2513  std::vector<int> col_offset(smaa.index(1).nblock());
2514  for (std::set<int>::const_iterator i = col_blocks.begin();
2515  i!=col_blocks.end(); i++) {
2516  col_offset[*i] = off;
2517  off+=smaa.index(1).block_size(*i);
2518  }
2519 
2520  // place the data in the array
2521  const Array<2>::blockmap_t &Abm = smaa.blockmap();
2522  for (Array<2>::blockmap_t::const_iterator i = Abm.begin();
2523  i != Abm.end();
2524  i++) {
2525  const BlockInfo<2> &Abi = i->first;
2526  int off0 = row_offset[Abi.block(0)];
2527  int off1 = col_offset[Abi.block(1)];
2528  BlockIter<2> j(smaa.indices(), Abi);
2529  double *dat = i->second;
2530  for (j.start(); j.ready(); j++) {
2531  dat[j.offset()] = m(off0+j.index(0), off1+j.index(1));
2532  }
2533  }
2534  smaa.compute_bounds();
2535  }
2536 
2540  template <class T>
2541  void pack_matrix_into_empty_array(const T &m, Array<2> &smaa,
2542  double bound) {
2543  smaa.clear();
2544  const Range &index0 = smaa.index(0);
2545  const Range &index1 = smaa.index(1);
2546  for (int i0=0; i0 < index0.nblock(); i0++) {
2547  int bs0 = index0.block_size(i0);
2548  int bo0 = index0.block_offset(i0);
2549  for (int i1=0; i1 < index1.nblock(); i1++) {
2550  int bs1 = index1.block_size(i1);
2551  int bo1 = index1.block_offset(i1);
2552  double maxval = 0.0;
2553  for (int j0=0; j0<bs0; j0++) {
2554  for (int j1=0; j1<bs1; j1++) {
2555  double val = fabs(m(bo0+j0,bo1+j1));
2556  if (val > maxval) maxval = val;
2557  }
2558  }
2559  if (maxval >= bound) {
2560  BlockInfo<2> bi;
2561  bi.block(0) = i0;
2562  bi.block(1) = i1;
2563  double *data = smaa.add_unallocated_block(bi).second;
2564  }
2565  }
2566  }
2567  pack_matrix_into_array(m,smaa);
2568  }
2569 
2573  template <class T>
2574  void pack_vector_into_empty_array(const T &m, Array<1> &smaa,
2575  double bound) {
2576  if (m.dim().n() != smaa.index(0).nindex()) {
2577  throw std::runtime_error("pack_vector_into_empty_array: dims don't match");
2578  }
2579  smaa.clear();
2580  const Range &index0 = smaa.index(0);
2581  for (int i0=0; i0 < index0.nblock(); i0++) {
2582  int bs0 = index0.block_size(i0);
2583  int bo0 = index0.block_offset(i0);
2584  double maxval = 0.0;
2585  for (int j0=0; j0<bs0; j0++) {
2586  double val = fabs(m(bo0+j0));
2587  if (val > maxval) maxval = val;
2588  }
2589  if (maxval >= bound) {
2590  BlockInfo<1> bi;
2591  bi.block(0) = i0;
2592  double *data = smaa.add_unallocated_block(bi).second;
2593  }
2594  }
2595  pack_vector_into_array(m,smaa);
2596  }
2597 
2603  template <class T>
2604  void pack_submatrix_into_empty_array(const T &m, Array<2> &smaa,
2605  double bound,
2606  const std::vector<int> &index_map_i,
2607  const std::vector<int> &index_map_j) {
2608  smaa.clear();
2609  const Range &index0 = smaa.index(0);
2610  const Range &index1 = smaa.index(1);
2611  std::vector<int> reverse_map_i(index0.nindex());
2612  std::vector<int> reverse_map_j(index1.nindex());
2613  std::fill(reverse_map_i.begin(), reverse_map_i.end(), -1);
2614  std::fill(reverse_map_j.begin(), reverse_map_j.end(), -1);
2615  for (int i=0; i<index_map_i.size(); i++) {
2616  reverse_map_i[index_map_i[i]] = i;
2617  }
2618  for (int i=0; i<index_map_j.size(); i++) {
2619  reverse_map_j[index_map_j[i]] = i;
2620  }
2621  for (int i0=0; i0 < index0.nblock(); i0++) {
2622  int bs0 = index0.block_size(i0);
2623  int bo0 = index0.block_offset(i0);
2624  for (int i1=0; i1 < index1.nblock(); i1++) {
2625  int bs1 = index1.block_size(i1);
2626  int bo1 = index1.block_offset(i1);
2627  double maxval = 0.0;
2628  for (int j0=0; j0<bs0; j0++) {
2629  int i_m = reverse_map_i[bo0+j0];
2630  if (i_m == -1) continue;
2631  for (int j1=0; j1<bs1; j1++) {
2632  int j_m = reverse_map_j[bo1+j1];
2633  if (j_m == -1) continue;
2634  double val = fabs(m(i_m, j_m));
2635  if (val > maxval) maxval = val;
2636  }
2637  }
2638  if (maxval >= bound) {
2639  BlockInfo<2> bi;
2640  bi.block(0) = i0;
2641  bi.block(1) = i1;
2642  double *data = smaa.add_unallocated_block(bi).second;
2643  }
2644  }
2645  }
2646  pack_submatrix_into_array(m,smaa,index_map_i,index_map_j);
2647  }
2648 
2654  template <class T>
2655  void pack_subvector_into_empty_array(const T &m, Array<1> &smaa,
2656  double bound,
2657  const std::vector<int> &index_map) {
2658  smaa.clear();
2659  const Range &index0 = smaa.index(0);
2660  std::vector<int> reverse_map(index0.nindex());
2661  std::fill(reverse_map.begin(), reverse_map.end(), -1);
2662  for (int i=0; i<index_map.size(); i++) {
2663  reverse_map[index_map[i]] = i;
2664  }
2665  for (int i0=0; i0 < index0.nblock(); i0++) {
2666  int bs0 = index0.block_size(i0);
2667  int bo0 = index0.block_offset(i0);
2668  double maxval = 0.0;
2669  for (int j0=0; j0<bs0; j0++) {
2670  int i_m = reverse_map[bo0+j0];
2671  if (i_m == -1) continue;
2672  double val = fabs(m(i_m));
2673  if (val > maxval) maxval = val;
2674  }
2675  if (maxval >= bound) {
2676  BlockInfo<1> bi;
2677  bi.block(0) = i0;
2678  double *data = smaa.add_unallocated_block(bi).second;
2679  }
2680  }
2681  pack_subvector_into_array(m,smaa,index_map);
2682  }
2683 
2687  template <class T>
2688  void pack_array_into_matrix(const Array<2> &smaa, T &m) {
2689  m.assign(0.0); // Initialize matrix to zero (not all blocks will be assigned below)
2690  const Array<2>::blockmap_t &Abm = smaa.blockmap();
2691  for (Array<2>::blockmap_t::const_iterator i = Abm.begin();
2692  i != Abm.end();
2693  i++) {
2694  const BlockInfo<2> &Abi = i->first;
2695  int off0 = smaa.index(0).block_offset(Abi.block(0));
2696  int off1 = smaa.index(1).block_offset(Abi.block(1));
2697  BlockIter<2> j(smaa.indices(), Abi);
2698  double *dat = i->second;
2699  for (j.start(); j.ready(); j++) {
2700  m(off0+j.index(0), off1+j.index(1)) = dat[j.offset()];
2701  }
2702  }
2703  }
2704 
2708  template <class T>
2709  void pack_array_into_vector(const Array<1> &smaa, T &m) {
2710  m.assign(0.0);
2711  const Array<1>::blockmap_t &Abm = smaa.blockmap();
2712  for (Array<1>::blockmap_t::const_iterator i = Abm.begin();
2713  i != Abm.end();
2714  i++) {
2715  const BlockInfo<1> &Abi = i->first;
2716  int off0 = smaa.index(0).block_offset(Abi.block(0));
2717  BlockIter<1> j(smaa.indices(), Abi);
2718  double *dat = i->second;
2719  for (j.start(); j.ready(); j++) {
2720  m(off0+j.index(0)) = dat[j.offset()];
2721  }
2722  }
2723  }
2724 
2730  template <class T>
2731  void pack_array_into_subvector(const Array<1> &smaa, T &m,
2732  const std::vector<int> &index_map) {
2733  m.assign(0.0);
2734  const Array<1>::blockmap_t &Abm = smaa.blockmap();
2735  for (Array<1>::blockmap_t::const_iterator i = Abm.begin();
2736  i != Abm.end();
2737  i++) {
2738  const BlockInfo<1> &Abi = i->first;
2739  int off0 = smaa.index(0).block_offset(Abi.block(0));
2740  BlockIter<1> j(smaa.indices(), Abi);
2741  double *dat = i->second;
2742  for (j.start(); j.ready(); j++) {
2743  m(index_map[off0+j.index(0)]) = dat[j.offset()];
2744  }
2745  }
2746  }
2747 
2753  template <class T>
2754  void pack_array_into_submatrix(const Array<2> &smaa, T &m,
2755  const std::vector<int> &index_map_i,
2756  const std::vector<int> &index_map_j) {
2757  m.assign(0.0); // Initialize matrix to zero (not all blocks will be assigned below)
2758  const Array<2>::blockmap_t &Abm = smaa.blockmap();
2759  for (Array<2>::blockmap_t::const_iterator i = Abm.begin();
2760  i != Abm.end();
2761  i++) {
2762  const BlockInfo<2> &Abi = i->first;
2763  int off0 = smaa.index(0).block_offset(Abi.block(0));
2764  int off1 = smaa.index(1).block_offset(Abi.block(1));
2765  BlockIter<2> j(smaa.indices(), Abi);
2766  double *dat = i->second;
2767  for (j.start(); j.ready(); j++) {
2768  m(index_map_i[off0+j.index(0)], index_map_j[off1+j.index(1)]) = dat[j.offset()];
2769  }
2770  }
2771  }
2772 
2775  void pack_2e_integrals_into_array(Array<4> & array,
2776  const sc::Ref<sc::TwoBodyInt> &tbint);
2777 
2781  void pack_2e_integrals_into_shell_blocked_array(Array<4> & array,
2782  const sc::Ref<sc::TwoBodyInt> &tbint);
2783 
2784 }
2785 }
2786 
2787 #include "arraydef.h"
2788 #include "blockinfodef.h"
2789 #include "contract.h"
2790 #include "contractpartdef.h"
2791 
2792 #endif
2793 
sc::sma2::ContractPart::value
double value()
Extract a value from the array.
Definition: contractpartdef.h:705
sc::sma2::Array::init_blocks
void init_blocks(const Array< N > &a)
Initialize the stored blocks to be the same as those in a.
Definition: sma.h:1897
sc::sma2::Array::set_use_blockmap_cache
void set_use_blockmap_cache(bool ubmc)
Used to indicate whether or not a blockmap cache is to be used.
Definition: sma.h:2022
sc::sma2::Array::n_element
size_t n_element() const
Returns the number of elements contained in blocks.
Definition: sma.h:1692
sc::sma2::BlockInfo< 4 >::assign_blocks
void assign_blocks(const IndexList &il, const BlockInfo< N2 > &bi2, const IndexList &il2)
Assign blocks to those in another BlockInfo, bi2, given an IndexList that specifies the index mapping...
Definition: sma.h:862
sc::sma2::Array::assign_all
void assign_all(const Array< N > &a)
Assigns this to 'a', and keeps all blocks.
Definition: sma.h:1342
sc::sma2::Range::index_to_offset
int index_to_offset(int i) const
Given an index, return its offset within the block in which that index is found.
Definition: sma.h:138
sc::sma2::Array::remove_block
void remove_block(const BlockInfo< N > &b)
Removes a block.
Definition: sma.h:1584
sc::auto_vec::get
T * get() const MPQC__NOEXCEPT
Returns the pointer.
Definition: autovec.h:82
sc::sma2::Range::max_block_size
int max_block_size() const
Returns the maximum block size.
Definition: sma.h:131
sc::sma2::Array::index
const Range & index(int i) const
Gets the i'th Range.
Definition: sma.h:1503
sc::sma2::BlockInfo< 3 >::subset_size
unsigned int subset_size(const Range *indices, const IndexList &indexlist) const
Compute the size of a block formed from this block by using some subset of the indices given by index...
Definition: sma.h:684
sc::sma2::Array::compute_bounds
void compute_bounds()
Computes the bound for each block.
Definition: sma.h:1743
sc::sma2::Range::basis_index_to_range_index
int basis_index_to_range_index(int i) const
Maps a basis function number to the number within this Range object.
sc::sma2::ContractPart::skip_bounds_update
ContractPart< N > skip_bounds_update() const
Causes the bounds to not be computed for the LHS operand for certain operations.
Definition: contractpartdef.h:696
sc::sma2::Array::add_new_unallocated_block
blockmap_t::iterator add_new_unallocated_block(const BlockInfo< N > &b)
Adds block b.
Definition: sma.h:1630
sc::sma2::BlockInfo< 3 >::assign_blocks
void assign_blocks(const IndexList &il, const BlockInfo< N2 > &bi2, const IndexList &il2)
Assign blocks to those in another BlockInfo, bi2, given an IndexList that specifies the index mapping...
Definition: sma.h:698
sc::sma2::Data
Data holds the values for each block.
Definition: sma.h:1026
sc::sma2::BlockInfo::equiv_blocks
bool equiv_blocks(const IndexList &il, const BlockInfo< N2 > &bi2)
Return true if blocks are the same as in another BlockInfo, bi2, given an IndexList that specifies th...
Definition: sma.h:288
sc::sma2::Array::operator*=
void operator*=(double f)
Multiplication by a scalar.
Definition: sma.h:1870
sc::Ref< sc::GaussianBasisSet >
sc::StateIn::getstring
virtual int getstring(char *&)
This restores strings saved with StateOut::putstring.
sc::sma2::Array::add_unallocated_block
blockmap_t::value_type & add_unallocated_block(const BlockInfo< N > &b)
Adds block b.
Definition: sma.h:1516
sc::sma2::Index::operator==
bool operator==(const Index &i) const
Returns true if the indices are equivalent.
Definition: sma.h:1082
sc::sma2::Array::Array
Array(const std::string &name="", double tol=DBL_EPSILON)
This does not initialize any indices.
Definition: sma.h:1352
sc::sma2::BlockIter::block_size
int block_size(int i) const
Returns the size of the given block.
Definition: sma.h:2241
sc::auto_vec::reset
void reset(T *d=0) MPQC__NOEXCEPT
Assign to a new value.
Definition: autovec.h:98
sc::sma2::Array::Array
Array(const Range &i0, const Range &i1, const Range &i2, const Range &i3, const Range &i4, const std::string &name="", double tol=DBL_EPSILON)
Initialize range 0, 1, 2, 3, and 4 to i0, i1, i2, i3, and i4, and the rest to i4.
Definition: sma.h:1417
sc::sma2::Array::blockmap_cache_entry_exists
bool blockmap_cache_entry_exists(const IndexList &il)
Return true if the needed blockmap cache entry exists.
Definition: sma.h:2030
sc::sma2::BlockInfo< 4 >::equiv_blocks
bool equiv_blocks(const IndexList &il, const BlockInfo< N2 > &bi2)
Return true if blocks are the same as in another BlockInfo, bi2, given an IndexList that specifies th...
Definition: sma.h:882
sc::sma2::ContractPart::operator|=
void operator|=(const ContractPart< N2 > &o) const
Add blocks to this corresponding the blocks already allocated in o.
Definition: contractpartdef.h:549
sc::sma2::Array::add_new_unallocated_block
blockmap_t::iterator add_new_unallocated_block(const typename blockmap_t::iterator &hint, const BlockInfo< N > &b)
Adds block b.
Definition: sma.h:1605
sc::sma2::Array::add_zeroed_block
blockmap_t::value_type & add_zeroed_block(const BlockInfo< N > &b)
Adds block b.
Definition: sma.h:1561
sc::sma2::BlockIter::start
void start()
Initialize the iterator.
Definition: sma.h:2211
sc::sma2::Index::symbolically_equivalent
bool symbolically_equivalent(const Index &i) const
Returns true if the indices have non-empty names which are the same.
Definition: sma.h:1079
sc::sma2::Range::all_size_one
bool all_size_one() const
Returns true if all blocks are size 1.
sc::sma2::Array::parallel_union
void parallel_union(const sc::Ref< sc::MessageGrp > &grp)
Makes the block distribution the same on all processes.
Definition: parallel.h:74
sc::auto_vec
The auto_vec class functions much like auto_ptr, except it contains references to arrays.
Definition: autovec.h:59
sc::sma2::Array::Array
Array(const Range &i0, const Range &i1, const Range &i2, const std::string &name="", double tol=DBL_EPSILON)
Initialize range 0 to i0, 1 to i1, and the rest to i2.
Definition: sma.h:1384
sc::sma2::Index::value
int value() const
Return the value of the index.
Definition: sma.h:1071
sc::sma2::BlockInfo< 4 >::zero
void zero()
Set all block indices to zero.
Definition: sma.h:890
sc::sma2::BlockInfo
BlockInfo stores info about a block of data.
Definition: sma.h:200
sc::sma2::Array::assign_tol
void assign_tol(const Array< N > &a, double tol)
Assigns this to 'a', but throws out blocks smaller than tol.
Definition: sma.h:1321
sc::sma2::BlockInfo::size
unsigned int size(const Range *indices) const
Compute the size of this block.
Definition: sma.h:247
sc::sma2::ContractProd
Represents a pairs of contracted array and their symbolic indices.
Definition: sma.h:1089
sc::sma2::BlockIter
BlockIter loops through the all the indices within a block.
Definition: sma.h:2189
sc::sma2::BlockInfo< 4 >
Definition: sma.h:785
sc::sma2::BlockInfo< 3 >::size
unsigned int size(const Range *indices) const
Compute the size of this block.
Definition: sma.h:677
sc::StateIn
Definition: statein.h:79
sc::sma2::Array::init_blocks
void init_blocks(const Array< N > &a, double tol)
Initialize the stored blocks to be the same as those in a.
Definition: sma.h:1885
sc::sma2::IndexListLess::IndexListLess
IndexListLess(const IndexList &il1, const IndexList &il2)
Use the indices in both IndexList objects, il1 and il2, to sort the indices.
Definition: sma.h:985
sc::sma2::Array::nindex
static int nindex()
Return the number of indices.
Definition: sma.h:1488
sc::sma2::Array
Implements a block sparse tensor.
Definition: sma.h:1088
sc::Ref::null
bool null() const
Return true if this is a reference to a null object.
Definition: ref.h:423
sc::sma2::BlockInfo::assign_blocks
void assign_blocks(const IndexList &il, const BlockInfo< N2 > &bi2)
Assign blocks to those in another BlockInfo, bi2, given an IndexList that specifies the index mapping...
Definition: sma.h:278
sc::sma2::Array::Array
Array(const Range &i0, const Range &i1, const std::string &name="", double tol=DBL_EPSILON)
Initialize range 0 to i0, and the rest to i1.
Definition: sma.h:1371
sc::sma2::ContractUnion
Definition: sma.h:1090
sc::sma2::Range
An Range represent a set of integers, [0, N).
Definition: sma.h:93
sc::sma2::Array::value
double value()
If this is a scalar (N==0) return the value of the scalar.
Definition: sma.h:2002
sc::sma2::Array::deallocate_blocks
void deallocate_blocks()
Deallocate storage for all blocks.
Definition: sma.h:1801
sc::sma2::Range::index_to_block
int index_to_block(int i) const
Given an index, return the block in which that index is found.
Definition: sma.h:135
sc::sma2::ContractPart::operator~
ContractPart< N > operator~() const
Instruct the contract routine to clear this array immediately after use.
Definition: contractpartdef.h:688
sc::sma2::Array::add_all_unallocated_blocks
void add_all_unallocated_blocks()
Adds all possible blocks to give a dense array.
Definition: sma.h:1653
sc::sma2::Array::use_blockmap_cache
bool use_blockmap_cache() const
Return true if blockmap caches are in use.
Definition: sma.h:2026
sc::sma2::Range::block_size
int block_size(int i) const
Returns the size for the given block.
Definition: sma.h:129
sc::sma2::Array::block_size
int block_size(const BlockInfo< N > &bi) const
Return size of the given block.
Definition: sma.h:1507
sc::operator<<
std::vector< unsigned int > operator<<(const GaussianBasisSet &B, const GaussianBasisSet &A)
computes a map from basis functions in A to the equivalent basis functions in B.
sc::sma2::Array::assign
void assign(double val)
Assigns the given value to all elements in the array.
Definition: sma.h:1722
sc::sma2::IndexList
An IndexList is a vector of indices.
Definition: sma.h:160
sc::sma2::Array::clear
void clear()
Make this array store nothing.
Definition: sma.h:1475
sc::sma2::BlockInfo< 0 >::subset_size
unsigned int subset_size(const Range *indices, const IndexList &indexlist) const
Compute the size of a block formed from this block by using some subset of the indices given by index...
Definition: sma.h:410
sc::sma2::Array::Array
Array(const Range &index, const std::string &name="", double tol=DBL_EPSILON)
Initialize all N indices to index.
Definition: sma.h:1359
sc::sma2::ContractPart
Represents an array and symbolic indices in a contraction.
Definition: sma.h:1095
mpqc::ci::write
void write(ci::Vector &V, File::Dataspace< double > F, const std::vector< mpqc::range > &local)
write local segments of V to F
Definition: direct.hpp:41
sc::sma2::BlockInfo< 0 >::zero
void zero()
Set all block indices to zero.
Definition: sma.h:438
sc::sma2::Array::distributed_from_distributed
void distributed_from_distributed(const sc::Ref< sc::MessageGrp > &, const BlockDistrib< N > &, Array< N > &, bool clear_source_array=false, bool ignore_block_distrib_throws=false)
Redistribute an array.
Definition: parallel.h:359
sc::sma2::Array::debug
int & debug()
Set/get the debug level.
Definition: sma.h:1928
sc::sma2::Array::Array
Array(const Range &i0, const Range &i1, const Range &i2, const Range &i3, const std::string &name="", double tol=DBL_EPSILON)
Initialize range 0, 1 and 2 to i0, i1, and i2, and the rest to i3.
Definition: sma.h:1399
sc::sma2::Range::nindex
int nindex() const
Returns the dimension of this range.
Definition: sma.h:125
sc::sma2::Array::zero
void zero()
Zeros the array.
Definition: sma.h:1739
sc::sma2::Array::set_index
void set_index(int i, const Range &idx)
Sets the i'th Range to index.
Definition: sma.h:1496
sc::sma2::Index::Index
Index(int value)
Allocate an index with a specific value.
Definition: sma.h:1059
sc::sma2::BlockInfo< 3 >::zero
void zero()
Set all block indices to zero.
Definition: sma.h:726
sc::StateOut
Definition: stateout.h:71
sc::sma2::IndexListLess::IndexListLess
IndexListLess(const IndexList &il)
Use the IndexList il to sort the indices.
Definition: sma.h:982
sc::sma2::Array::set_indices
void set_indices(const Range *r)
Initialize each range to the ranges in the given range array.
Definition: sma.h:1499
sc::AVLMMap< BlockInfo< N >, double *, Compare >
sc::sma2::BlockInfo< 0 >::assign_blocks
void assign_blocks(const IndexList &il, const BlockInfo< N2 > &bi2, const IndexList &il2)
Assign blocks to those in another BlockInfo, bi2, given an IndexList that specifies the index mapping...
Definition: sma.h:419
sc::sma2::BlockInfo< 0 >::size
unsigned int size(const Range *indices) const
Compute the size of this block.
Definition: sma.h:405
sc::sma2::BlockInfo< 4 >::assign_blocks
void assign_blocks(const IndexList &il, const BlockInfo< N2 > &bi2)
Assign blocks to those in another BlockInfo, bi2, given an IndexList that specifies the index mapping...
Definition: sma.h:872
sc::sma2::Array::n_element_allocated
size_t n_element_allocated() const
Returns the number of elements contained in blocks.
Definition: sma.h:1702
sc::sma2::IndicesLess
Functor for comparing a block's indices.
Definition: sma.h:333
sc::sma2::Array::Array
Array(const Range &i0, const Range &i1, const Range &i2, const Range &i3, const Range &i4, const Range &i5, const std::string &name="", double tol=DBL_EPSILON)
Initialize range 0, 1, 2, 3, 4, and 5 to i0, i1, i2, i3, i4, and i5, and the rest to i5.
Definition: sma.h:1438
sc::sma2::Array::block_max_abs
double block_max_abs(typename blockmap_t::const_iterator &b) const
Find maximum absolute value of elements in a block.
Definition: sma.h:1901
sc::sma2::BlockInfo< 3 >::assign_blocks
void assign_blocks(const IndexList &il, const BlockInfo< N2 > &bi2)
Assign blocks to those in another BlockInfo, bi2, given an IndexList that specifies the index mapping...
Definition: sma.h:708
sc::ProgrammingError
This is thrown when a situations arises that should be impossible.
Definition: scexception.h:92
sc::sma2::BlockInfo< 4 >::size
unsigned int size(const Range *indices) const
Compute the size of this block.
Definition: sma.h:841
sc::sma2::Index::set_value
void set_value(int value)
Set the value of the index.
Definition: sma.h:1085
sc::sma2::Array::blockmap_cache_entry
cached_blockmap_t & blockmap_cache_entry(const IndexList &il)
Return a cached blockmap.
Definition: sma.h:2034
sc::sma2::Array::replicated_from_distributed
void replicated_from_distributed(const sc::Ref< sc::MessageGrp > &, const Array< N > &)
Convert a distributed array to a replicated array.
Definition: parallel.h:309
sc::sma2::BlockInfo< 3 >
Definition: sma.h:626
sc::sma2::BlockIter::operator++
void operator++(int)
Increments the iterator.
Definition: sma.h:2215
sc::sma2::Array::allocate_blocks
void allocate_blocks()
Allocate storage for all blocks if storage has not already been allocated.
Definition: sma.h:1764
sc::sma2::Array::initial_hint
blockmap_t::iterator initial_hint()
Returns an initial hint.
Definition: sma.h:1684
mpqc::ci::read
void read(ci::Vector &V, File::Dataspace< double > F, const std::vector< mpqc::range > &local)
read local segments into V from F
Definition: direct.hpp:25
sc::sma2::Index::Index
Index(const std::string &name, int value)
Allocate an Index with a fixed value and a name.
Definition: sma.h:1065
sc::sma2::IndexListLess
Functor for determining if one IndexList is less than another.
Definition: sma.h:977
sc::ExEnv::out0
static std::ostream & out0()
Return an ostream that writes from node 0.
sc::sma2::BlockInfo::assign_blocks
void assign_blocks(const IndexList &il, const BlockInfo< N2 > &bi2, const IndexList &il2)
Assign blocks to those in another BlockInfo, bi2, given an IndexList that specifies the index mapping...
Definition: sma.h:268
sc::sma2::Index::has_value
bool has_value() const
Returns true if the index has a value, indicating that the index is to be fixed.
Definition: sma.h:1074
sc::sma2::Array::print_local
void print_local(std::ostream &o=sc::ExEnv::outn()) const
Print the array.
Definition: arraydef.h:34
sc::sma2::Array::Array
Array(const Range *ranges, double tol=DBL_EPSILON)
Initialize each range to the ranges in the given range array.
Definition: sma.h:1460
sc::sma2::BlockInfo::subset_size
unsigned int subset_size(const Range *indices, const IndexList &indexlist) const
Compute the size of a block formed from this block by using some subset of the indices given by index...
Definition: sma.h:254
sc::StateOut::put
virtual int put(const ClassDesc *)
Write out information about the given ClassDesc.
sc::sma2::BlockInfo< 3 >::equiv_blocks
bool equiv_blocks(const IndexList &il, const BlockInfo< N2 > &bi2)
Return true if blocks are the same as in another BlockInfo, bi2, given an IndexList that specifies th...
Definition: sma.h:718
sc::RefCount
The base class for all reference counted objects.
Definition: ref.h:192
sc::sma2::BlockInfo< 4 >::subset_size
unsigned int subset_size(const Range *indices, const IndexList &indexlist) const
Compute the size of a block formed from this block by using some subset of the indices given by index...
Definition: sma.h:848
sc::StateIn::get
virtual int get(const ClassDesc **)
This restores ClassDesc's.
sc::ExEnv::outn
static std::ostream & outn()
Return an ostream that writes from all nodes.
Definition: exenv.h:81
sc::sma2::Index::Index
Index(const std::string &name)
Allocate a named index.
Definition: sma.h:1053
sc::sma2::BlockInfo< 0 >::assign_blocks
void assign_blocks(const IndexList &il, const BlockInfo< N2 > &bi2)
Assign blocks to those in another BlockInfo, bi2, given an IndexList that specifies the index mapping...
Definition: sma.h:426
sc::StateOut::putstring
virtual int putstring(const char *)
This is like put except the length of the char array is determined by interpreting the character arra...
sc::sma2::Array::n_block
int n_block() const
Returns the number of blocks.
Definition: sma.h:1688
sc::sma2::BlockInfo::BlockInfo
BlockInfo(const BlockInfo< N > &b, const IndexList &l)
Initialize the blocks to the values in the argument, b.
Definition: sma.h:223
sc::sma2::Array::add_allocated_block
blockmap_t::value_type & add_allocated_block(const BlockInfo< N > &b)
Adds block b.
Definition: sma.h:1541
sc::sma2::Index::name
const std::string name() const
Return the name of the index.
Definition: sma.h:1068
sc::sma2::Array::indices
const Range * indices() const
Gets the Range array.
Definition: sma.h:1505
sc::sma2::BlockInfo< 0 >::equiv_blocks
bool equiv_blocks(const IndexList &il, const BlockInfo< N2 > &bi2)
Return true if blocks are the same as in another BlockInfo, bi2, given an IndexList that specifies th...
Definition: sma.h:433
sc::sma2::Array::blockmap
const blockmap_t & blockmap() const
Return the block map.
Definition: sma.h:1490
sc::sma2::BlockInfo::BlockInfo
BlockInfo(const std::vector< bi_t > &v)
Initialize the blocks to the values in the argument, v.
Definition: sma.h:214
sc::sma2::IndexList::is_identity
bool is_identity() const
Returns true if this is 0, 1, ..., n-1.
sc::sma2::Range::nblock
int nblock() const
Returns the number of blocks.
Definition: sma.h:127
sc::sma2::Index
An Index is used in the symbolic notation for contractions.
Definition: sma.h:1047
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
sc::sma2::Array::operator=
Array< N > & operator=(const Array< N > &a)
Assigns this to 'a', but throws out small blocks.
Definition: sma.h:1346
sc::sma2::BlockIter::ready
bool ready() const
Returns true if the iterator is valid.
Definition: sma.h:2213
sc::sma2::Array::parallel_accumulate
void parallel_accumulate(const sc::Ref< sc::MessageGrp > &grp)
Performs a reduce-broadcast on the data.
Definition: parallel.h:50
sc::sma2::Array::max_abs_element
double max_abs_element()
Find maximum absolute value of elements.
Definition: sma.h:1912
sc::sma2::BlockInfo::zero
void zero()
Set all block indices to zero.
Definition: sma.h:239
sc::sma2::IndexList::is_identity_permutation
bool is_identity_permutation() const
Returns true if this is a permutation of 0, 1, ..., n-1.
sc::sma2::Range::block_offset
int block_offset(int i) const
Returns the offset for the given block.
Definition: sma.h:133
sc::sma2::BlockInfo< 0 >
Definition: sma.h:361
sc::sma2::Array::clear_blockmap_cache
void clear_blockmap_cache()
Get rid of cached blockmaps.
Definition: sma.h:2008
Generated at Sun Jan 26 2020 23:23:58 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.