MPQC 3.0.0-alpha
Loading...
Searching...
No Matches
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
85namespace sc {
86
87namespace 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
1493 const blockhash_t &blockhash() const { return block_hash_; }
1494#endif
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::AVLMMap
Definition avlmmap.h:155
sc::ExEnv::outn
static std::ostream & outn()
Return an ostream that writes from all nodes.
Definition exenv.h:81
sc::ExEnv::out0
static std::ostream & out0()
Return an ostream that writes from node 0.
sc::ProgrammingError
This is thrown when a situations arises that should be impossible.
Definition scexception.h:92
sc::RefCount
The base class for all reference counted objects.
Definition ref.h:192
sc::Ref
A template class that maintains references counts.
Definition ref.h:361
sc::Ref::null
bool null() const
Return true if this is a reference to a null object.
Definition ref.h:423
sc::StateIn
Restores fundamental and user-defined types from images created with StateOut.
Definition statein.h:79
sc::StateIn::get
virtual int get(const ClassDesc **)
This restores ClassDesc's.
sc::StateIn::getstring
virtual int getstring(char *&)
This restores strings saved with StateOut::putstring.
sc::StateOut
Serializes fundamental and user-defined types.
Definition stateout.h:71
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::StateOut::put
virtual int put(const ClassDesc *)
Write out information about the given ClassDesc.
sc::auto_vec
The auto_vec class functions much like auto_ptr, except it contains references to arrays.
Definition autovec.h:59
sc::auto_vec::get
T * get() const MPQC__NOEXCEPT
Returns the pointer.
Definition autovec.h:82
sc::auto_vec::reset
void reset(T *d=0) MPQC__NOEXCEPT
Assign to a new value.
Definition autovec.h:98
sc::sma2::Array
Implements a block sparse tensor.
Definition sma.h:1247
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::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::Array::deallocate_blocks
void deallocate_blocks()
Deallocate storage for all blocks.
Definition sma.h:1801
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::Array::use_blockmap_cache
bool use_blockmap_cache() const
Return true if blockmap caches are in use.
Definition sma.h:2026
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::indices
const Range * indices() const
Gets the Range array.
Definition sma.h:1505
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:360
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::Array::clear
void clear()
Make this array store nothing.
Definition sma.h:1475
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::Array::initial_hint
blockmap_t::iterator initial_hint()
Returns an initial hint.
Definition sma.h:1684
sc::sma2::Array::parallel_accumulate
void parallel_accumulate(const sc::Ref< sc::MessageGrp > &grp)
Performs a reduce-broadcast on the data.
Definition parallel.h:51
sc::sma2::Array::block_size
int block_size(const BlockInfo< N > &bi) const
Return size of the given block.
Definition sma.h:1507
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::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::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 std::string &name="", double tol=DBL_EPSILON)
Initialize range 0 to i0, and the rest to i1.
Definition sma.h:1371
sc::sma2::Array::assign
void assign(double val)
Assigns the given value to all elements in the array.
Definition sma.h:1722
sc::sma2::Array::n_element
size_t n_element() const
Returns the number of elements contained in blocks.
Definition sma.h:1692
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::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::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::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::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::add_new_unallocated_block
blockmap_t::iterator add_new_unallocated_block(const BlockInfo< N > &b)
Adds block b.
Definition sma.h:1630
sc::sma2::Array::n_block
int n_block() const
Returns the number of blocks.
Definition sma.h:1688
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::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:310
sc::sma2::Array::remove_block
void remove_block(const BlockInfo< N > &b)
Removes a block.
Definition sma.h:1584
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::blockmap
const blockmap_t & blockmap() const
Return the block map.
Definition sma.h:1490
sc::sma2::Array::compute_bounds
void compute_bounds()
Computes the bound for each block.
Definition sma.h:1743
sc::sma2::Array::zero
void zero()
Zeros the array.
Definition sma.h:1739
sc::sma2::Array::index
const Range & index(int i) const
Gets the i'th Range.
Definition sma.h:1503
sc::sma2::Array::clear_blockmap_cache
void clear_blockmap_cache()
Get rid of cached blockmaps.
Definition sma.h:2008
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::Array::operator*=
void operator*=(double f)
Multiplication by a scalar.
Definition sma.h:1870
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:75
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::Array::max_abs_element
double max_abs_element()
Find maximum absolute value of elements.
Definition sma.h:1912
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::Array::value
double value()
If this is a scalar (N==0) return the value of the scalar.
Definition sma.h:2002
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::add_zeroed_block
blockmap_t::value_type & add_zeroed_block(const BlockInfo< N > &b)
Adds block b.
Definition sma.h:1561
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::Array::set_index
void set_index(int i, const Range &idx)
Sets the i'th Range to index.
Definition sma.h:1496
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::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::Array::nindex
static int nindex()
Return the number of indices.
Definition sma.h:1488
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::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::BlockDistrib
Provides information about how blocks are distributed onto processes.
Definition sma.h:1205
sc::sma2::BlockInfo< 0 >::zero
void zero()
Set all block indices to zero.
Definition sma.h:438
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< 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 >::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::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::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::BlockInfo
BlockInfo stores info about a block of data.
Definition sma.h:200
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::BlockInfo::zero
void zero()
Set all block indices to zero.
Definition sma.h:239
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::sma2::BlockInfo::size
unsigned int size(const Range *indices) const
Compute the size of this block.
Definition sma.h:247
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::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::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::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::BlockIter
BlockIter loops through the all the indices within a block.
Definition sma.h:2189
sc::sma2::BlockIter::ready
bool ready() const
Returns true if the iterator is valid.
Definition sma.h:2213
sc::sma2::BlockIter::operator++
void operator++(int)
Increments the iterator.
Definition sma.h:2215
sc::sma2::BlockIter::start
void start()
Initialize the iterator.
Definition sma.h:2211
sc::sma2::BlockIter::block_size
int block_size(int i) const
Returns the size of the given block.
Definition sma.h:2241
sc::sma2::ContractPart
Represents an array and symbolic indices in a contraction.
Definition sma.h:1095
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::ContractPart::value
double value()
Extract a value from the array.
Definition contractpartdef.h:705
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::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::ContractProd
Represents a pairs of contracted array and their symbolic indices.
Definition sma.h:2062
sc::sma2::ContractUnion
Definition sma.h:1191
sc::sma2::Data
Data holds the values for each block.
Definition sma.h:1026
sc::sma2::IndexListLess
Functor for determining if one IndexList is less than another.
Definition sma.h:977
sc::sma2::IndexListLess::IndexListLess
IndexListLess(const IndexList &il)
Use the IndexList il to sort the indices.
Definition sma.h:982
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::IndexList
An IndexList is a vector of indices.
Definition sma.h:160
sc::sma2::IndexList::IndexList
IndexList(const IndexList &il1, const IndexList &il2)
Append the indices from both IndexList objects, il1 and il2, to form this IndexList.
sc::sma2::IndexList::is_identity
bool is_identity() const
Returns true if this is 0, 1, ..., n-1.
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::Index
An Index is used in the symbolic notation for contractions.
Definition sma.h:1047
sc::sma2::Index::Index
Index(int value)
Allocate an index with a specific value.
Definition sma.h:1059
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::Index::operator==
bool operator==(const Index &i) const
Returns true if the indices are equivalent.
Definition sma.h:1082
sc::sma2::Index::Index
Index(const std::string &name)
Allocate a named index.
Definition sma.h:1053
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::Index::set_value
void set_value(int value)
Set the value of the index.
Definition sma.h:1085
sc::sma2::Index::value
int value() const
Return the value of the index.
Definition sma.h:1071
sc::sma2::Index::name
const std::string name() const
Return the name of the index.
Definition sma.h:1068
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::IndicesLess
Functor for comparing a block's indices.
Definition sma.h:333
sc::sma2::Range
An Range represent a set of integers, [0, N).
Definition sma.h:93
sc::sma2::Range::nindex
int nindex() const
Returns the dimension of this range.
Definition sma.h:125
sc::sma2::Range::block_size
int block_size(int i) const
Returns the size for the given block.
Definition sma.h:129
sc::sma2::Range::max_block_size
int max_block_size() const
Returns the maximum block size.
Definition sma.h:131
sc::sma2::Range::nblock
int nblock() const
Returns the number of blocks.
Definition sma.h:127
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::Range::all_size_one
bool all_size_one() const
Returns true if all blocks are size 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::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::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
Contains all MPQC code up to version 3.
Definition mpqcin.h:14
Generated at Wed Sep 25 2024 02:45:29 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.12.0.