QDP++
qdp_scalarvec_specific.h
Go to the documentation of this file.
1// -*- C++ -*-
2
6
7#ifndef QDP_SCALARVEC_SPECIFIC_H
8#define QDP_SCALARVEC_SPECIFIC_H
9
10namespace QDP {
11
12//-----------------------------------------------------------------------------
13// Layout stuff specific to a scalarvec architecture
14namespace Layout
15{
18
21}
22
23
24//-----------------------------------------------------------------------------
25// Internal ops designed to look like those in parscalar
26// These dummy routines exist just to make code more portable
27namespace QDPInternal
28{
30 template<class T>
31 inline void globalSumArray(T* dest, int n) {}
32
34 template<class T>
35 inline void globalSumArray(multi1d<T>& dest) {}
36
38 template<class T>
39 inline void globalSumArray(multi2d<T>& dest) {}
40
42 template<class T>
43 inline void globalSum(T& dest) {}
44
46 template<class T>
47 inline void broadcast(T& dest) {}
48
49 template<>
50 inline void broadcast(std::string& dest) {}
51
53 inline void broadcast_str(std::string& dest) {}
54
56 inline void broadcast(void* dest, size_t nbytes) {}
57}
58
59#define QDP_NOT_IMPLEMENTED
60
61//-----------------------------------------------------------------------------
63
67template<class T, class T1, class Op, class RHS>
68//inline
69void evaluate(OLattice<T>& dest, const Op& op, const QDPExpr<RHS,OScalar<T1> >& rhs,
70 const Subset& s)
71{
72// cerr << "In evaluateUnorderedSubet(olattice,oscalar)\n";
73
74#if defined(QDP_USE_PROFILING)
75 static QDPProfile_t prof(dest, op, rhs);
76 prof.time -= getClockTime();
77#endif
78
79#if ! defined(QDP_NOT_IMPLEMENTED)
80 const int *tab = s.siteTable().slice();
81 for(int j=0; j < s.numSiteTable(); ++j)
82 {
83 int i = tab[j];
84// fprintf(stderr,"eval(olattice,oscalar): site %d\n",i);
85// op(dest.elem(i), forEach(rhs, ElemLeaf(), OpCombine()));
86 op(dest.elem(i), forEach(rhs, EvalLeaf1(0), OpCombine()));
87 }
88#else
89 QDP_error_exit("evaluateSubset not implemented");
90#endif
91
92#if defined(QDP_USE_PROFILING)
93 prof.time += getClockTime();
94 prof.count++;
95 prof.print();
96#endif
97}
98
100
104template<class T, class T1, class Op, class RHS>
105//inline
106void evaluate(OLattice<T>& dest, const Op& op, const QDPExpr<RHS,OLattice<T1> >& rhs,
107 const Subset& s)
108{
109// cerr << "In evaluateSubset(olattice,olattice)" << endl;
110
111#if defined(QDP_USE_PROFILING)
112 static QDPProfile_t prof(dest, op, rhs);
113 prof.time -= getClockTime();
114#endif
115
116#if ! defined(QDP_NOT_IMPLEMENTED)
117 // General form of loop structure
118 const int *tab = s.siteTable().slice();
119 for(int j=0; j < s.numSiteTable(); ++j)
120 {
121 int i = tab[j];
122// fprintf(stderr,"eval(olattice,olattice): site %d\n",i);
123 op(dest.elem(i), forEach(rhs, EvalLeaf1(i), OpCombine()));
124 }
125#else
126 QDP_error_exit("evaluateSubset not implemented");
127#endif
128
129#if defined(QDP_USE_PROFILING)
130 prof.time += getClockTime();
131 prof.count++;
132 prof.print();
133#endif
134}
135
136
137
138
139//-----------------------------------------------------------------------------
141template<class T1, class T2>
142void
143copymask(OSubLattice<T2> d, const OLattice<T1>& mask, const OLattice<T2>& s1)
144{
145 OLattice<T2>& dest = d.field();
146 const Subset& s = d.subset();
147
148#if ! defined(QDP_NOT_IMPLEMENTED)
149 const int *tab = s.siteTable().slice();
150 for(int j=0; j < s.numSiteTable(); ++j)
151 {
152 int i = tab[j];
153 copymask(dest.elem(i), mask.elem(i), s1.elem(i));
154 }
155#else
156 QDP_error_exit("copymask_Subset not implemented");
157#endif
158}
159
160
161
163template<class T1, class T2>
164void
165copymask(OLattice<T2>& dest, const OLattice<T1>& mask, const OLattice<T2>& s1)
166{
167 const int iend = Layout::outerSitesOnNode();
168 for(int i=0; i < iend; ++i)
169 copymask(dest.elem(i), mask.elem(i), s1.elem(i));
170}
171
172
173
174//-----------------------------------------------------------------------------
175// Random numbers
176namespace RNG
177{
178 extern Seed ran_seed;
179 extern Seed ran_mult;
180 extern Seed ran_mult_n;
182}
183
184
186
187template<class T>
188void
189random(OScalar<T>& d)
190{
191 Seed seed = RNG::ran_seed;
192 Seed skewed_seed = RNG::ran_seed * RNG::ran_mult;
193
194 fill_random(d.elem(), seed, skewed_seed, RNG::ran_mult);
195
196 RNG::ran_seed = seed; // The seed from any site is the same as the new global seed
197}
198
199
201template<class T>
202void
203random(OLattice<T>& d, const Subset& s)
204{
205 Seed seed;
206 Seed skewed_seed;
207
208#if ! defined(QDP_NOT_IMPLEMENTED)
209#error "random(unorderedsubset) broken"
210 const int *tab = s.siteTable().slice();
211 for(int j=0; j < s.numSiteTable(); ++j)
212 {
213 int i = tab[j];
214 seed = RNG::ran_seed;
215 skewed_seed.elem() = RNG::ran_seed.elem() * RNG::lattice_ran_mult->elem(i);
216 fill_random(d.elem(i), seed, skewed_seed, RNG::ran_mult_n);
217 }
218
219 RNG::ran_seed = seed; // The seed from any site is the same as the new global seed
220#else
221 QDP_error_exit("random_Subset not implemented");
222#endif
223}
224
225
226
227
229template<class T>
230void random(const OSubLattice<T>& dd)
231{
232 OLattice<T>& d = const_cast<OSubLattice<T>&>(dd).field();
233 const Subset& s = dd.subset();
234
235 random(d,s);
236}
237
238
240template<class T>
241void random(OLattice<T>& d)
242{
243 random(d,all);
244}
245
246
248template<class T>
249void gaussian(OLattice<T>& d, const Subset& s)
250{
251 OLattice<T> r1, r2;
252
253 random(r1,s);
254 random(r2,s);
255
256#if ! defined(QDP_NOT_IMPLEMENTED)
257 const int *tab = s.siteTable().slice();
258 for(int j=0; j < s.numSiteTable(); ++j)
259 {
260 int i = tab[j];
261 fill_gaussian(d.elem(i), r1.elem(i), r2.elem(i));
262 }
263#else
264 QDP_error_exit("gaussianSubset not implemented");
265#endif
266}
267
268
270template<class T>
272{
273 OLattice<T>& d = const_cast<OSubLattice<T>&>(dd).field();
274 const Subset& s = dd.subset();
275
276 gaussian(d,s);
277}
278
279
281template<class T>
282void gaussian(OLattice<T>& d)
283{
284 gaussian(d,all);
285}
286
287
288
289//-----------------------------------------------------------------------------
290// Broadcast operations
292template<class T>
293void zero_rep(OLattice<T>& dest, const Subset& s)
294{
295#if ! defined(QDP_NOT_IMPLEMENTED)
296 const int *tab = s.siteTable().slice();
297 for(int j=0; j < s.numSiteTable(); ++j)
298 {
299 int i = tab[j];
300 zero_rep(dest.elem(i));
301 }
302#else
303 QDP_error_exit("zero_rep_Subset not implemented");
304#endif
305}
306
307
308
310template<class T>
312{
313 OLattice<T>& d = dd.field();
314 const Subset& s = dd.subset();
315
316 zero_rep(d,s);
317}
318
319
321template<class T>
322void zero_rep(OLattice<T>& dest)
323{
324 const int iend = Layout::outerSitesOnNode();
325 for(int i=0; i < iend; ++i)
326 zero_rep(dest.elem(i));
327}
328
329
330
331//-----------------------------------------------
332// Global sums
334
338template<class RHS, class T>
339typename UnaryReturn<OScalar<T>, FnSum>::Type_t
340sum(const QDPExpr<RHS,OScalar<T> >& s1, const Subset& s)
341{
342 typename UnaryReturn<OScalar<T>, FnSum>::Type_t d;
343
344#if defined(QDP_USE_PROFILING)
345 static QDPProfile_t prof(d, OpAssign(), FnSum(), s1);
346 prof.time -= getClockTime();
347#endif
348
349 evaluate(d,OpAssign(),s1,all); // since OScalar, no global sum needed
350
351#if defined(QDP_USE_PROFILING)
352 prof.time += getClockTime();
353 prof.count++;
354 prof.print();
355#endif
356
357 return d;
358}
359
360
362
366template<class RHS, class T>
367typename UnaryReturn<OScalar<T>, FnSum>::Type_t
368sum(const QDPExpr<RHS,OScalar<T> >& s1)
369{
370 typename UnaryReturn<OScalar<T>, FnSum>::Type_t d;
371
372#if defined(QDP_USE_PROFILING)
373 static QDPProfile_t prof(d, OpAssign(), FnSum(), s1);
374 prof.time -= getClockTime();
375#endif
376
377 evaluate(d,OpAssign(),s1,all);
378
379#if defined(QDP_USE_PROFILING)
380 prof.time += getClockTime();
381 prof.count++;
382 prof.print();
383#endif
384
385 return d;
386}
387
388
389
391
400template<class RHS, class T>
401typename UnaryReturn<OLattice<T>, FnSum>::Type_t
402sum(const QDPExpr<RHS,OLattice<T> >& s1, const Subset& s)
403{
404 typename UnaryReturn<OLattice<T>, FnSum>::Type_t d;
405 OScalar<T> tmp; // Note, expect to have ILattice inner grid
406
407#if defined(QDP_USE_PROFILING)
408 static QDPProfile_t prof(d, OpAssign(), FnSum(), s1);
409 prof.time -= getClockTime();
410#endif
411
412 // Must initialize to zero since we do not know if the loop will be entered
413 zero_rep(d.elem());
414
415 if (s.hasOrderedRep())
416 {
417 const int istart = s.start() >> INNER_LOG;
418 const int iend = s.end() >> INNER_LOG;
419
420 for(int i=istart; i <= iend; ++i)
421 {
422 tmp.elem() = forEach(s1, EvalLeaf1(i), OpCombine()); // Evaluate to ILattice part
423 d.elem() += sum(tmp.elem()); // sum as well the ILattice part
424 }
425 }
426 else
427 {
428 const int *tab = s.siteTable().slice();
429 for(int j=0; j < s.numSiteTable(); ++j)
430 {
431 int i = tab[j];
432 int outersite = i >> INNER_LOG;
433 int innersite = i & ((1 << INNER_LOG)-1);
434
435 tmp.elem() = forEach(s1, EvalLeaf1(outersite), OpCombine()); // Evaluate to ILattice part
436 d.elem() += getSite(tmp.elem(),innersite); // wasteful - only extract a single site worth
437 }
438 }
439
440#if defined(QDP_USE_PROFILING)
441 prof.time += getClockTime();
442 prof.count++;
443 prof.print();
444#endif
445
446 return d;
447}
448
449
450
451
452
454
458template<class RHS, class T>
459typename UnaryReturn<OLattice<T>, FnSum>::Type_t
460sum(const QDPExpr<RHS,OLattice<T> >& s1)
461{
462 return sum(s1,all);
463}
464
465
466//-----------------------------------------------------------------------------
467// Multiple global sums
469
476template<class RHS, class T>
477multi1d<typename UnaryReturn<OScalar<T>, FnSum>::Type_t>
478sumMulti(const QDPExpr<RHS,OScalar<T> >& s1, const Set& ss)
479{
481
482#if defined(QDP_USE_PROFILING)
483 static QDPProfile_t prof(dest(0,0), OpAssign(), FnSum(), s1);
484 prof.time -= getClockTime();
485#endif
486
487 // lazy - evaluate repeatedly
488 for(int i=0; i < ss.numSubsets(); ++i)
489 dest[i] = sum(s1,ss[i]);
490
491#if defined(QDP_USE_PROFILING)
492 prof.time += getClockTime();
493 prof.count++;
494 prof.print();
495#endif
496
497 return dest;
498}
499
500
502
510template<class RHS, class T>
511typename UnaryReturn<OLattice<T>, FnSumMulti>::Type_t
512sumMulti(const QDPExpr<RHS,OLattice<T> >& s1, const Set& ss)
513{
514 typename UnaryReturn<OLattice<T>, FnSumMulti>::Type_t dest(ss.numSubsets());
515
516#if defined(QDP_USE_PROFILING)
517 static QDPProfile_t prof(dest[0], OpAssign(), FnSum(), s1);
518 prof.time -= getClockTime();
519#endif
520
521 // lazy - evaluate repeatedly
522 for(int i=0; i < ss.numSubsets(); ++i)
523 dest[i] = sum(s1,ss[i]);
524
525#if defined(QDP_USE_PROFILING)
526 prof.time += getClockTime();
527 prof.count++;
528 prof.print();
529#endif
530
531 return dest;
532}
533
534
535//-----------------------------------------------------------------------------
536// Multiple global sums
538
545template<class T>
546multi2d<typename UnaryReturn<OScalar<T>, FnSum>::Type_t>
547sumMulti(const multi1d< OScalar<T> >& s1, const Set& ss)
548{
549 multi2d<typename UnaryReturn<OScalar<T>, FnSum>::Type_t> dest(s1.size(), ss.numSubsets());
550
551#if defined(QDP_USE_PROFILING)
552 static QDPProfile_t prof(dest(0,0), OpAssign(), FnSum(), s1);
553 prof.time -= getClockTime();
554#endif
555
556 // lazy - evaluate repeatedly
557 for(int i=0; i < dest.size1(); ++i)
558 for(int j=0; j < dest.size2(); ++j)
559 dest(j,i) = s1[j];
560
561#if defined(QDP_USE_PROFILING)
562 prof.time += getClockTime();
563 prof.count++;
564 prof.print();
565#endif
566
567 return dest;
568}
569
570
572
580template<class T>
581multi2d<typename UnaryReturn<OLattice<T>, FnSum>::Type_t>
582sumMulti(const multi1d< OLattice<T> >& s1, const Set& ss)
583{
584 multi2d<typename UnaryReturn<OLattice<T>, FnSum>::Type_t> dest(s1.size(),ss.numSubsets());
585
586#if defined(QDP_USE_PROFILING)
587 static QDPProfile_t prof(dest(0,0), OpAssign(), FnSum(), s1);
588 prof.time -= getClockTime();
589#endif
590
591 // lazy - evaluate repeatedly
592 for(int k=0; k < s1.size(); ++k)
593 for(int i=0; i < ss.numSubsets(); ++i)
594 dest(k,i) = sum(s1[k],ss[i]);
595
596#if defined(QDP_USE_PROFILING)
597 prof.time += getClockTime();
598 prof.count++;
599 prof.print();
600#endif
601
602 return dest;
603}
604
605
606//-----------------------------------------------------------------------------
608
614template<class T>
615inline typename UnaryReturn<OScalar<T>, FnNorm2>::Type_t
616norm2(const multi1d< OScalar<T> >& s1)
617{
618 typename UnaryReturn<OScalar<T>, FnNorm2>::Type_t d;
619
620#if defined(QDP_USE_PROFILING)
621 static QDPProfile_t prof(d, OpAssign(), FnNorm2(), s1[0]);
622 prof.time -= getClockTime();
623#endif
624
625 // Possibly loop entered
626 zero_rep(d.elem());
627
628 for(int n=0; n < s1.size(); ++n)
629 {
630 OScalar<T>& ss1 = s1[n];
631 d.elem() += localNorm2(ss1.elem());
632 }
633
634#if defined(QDP_USE_PROFILING)
635 prof.time += getClockTime();
636 prof.count++;
637 prof.print();
638#endif
639
640 return d;
641}
642
644
645template<class T>
646inline typename UnaryReturn<OScalar<T>, FnNorm2>::Type_t
647norm2(const multi1d< OScalar<T> >& s1, const Subset& s)
648{
649 return norm2(s1);
650}
651
652
653
655
661template<class T>
662inline typename UnaryReturn<OLattice<T>, FnNorm2>::Type_t
663norm2(const multi1d< OLattice<T> >& s1, const Subset& s)
664{
665 typename UnaryReturn<OLattice<T>, FnNorm2>::Type_t d;
666
667#if defined(QDP_USE_PROFILING)
668 static QDPProfile_t prof(d, OpAssign(), FnNorm2(), s1[0]);
669 prof.time -= getClockTime();
670#endif
671
672 // Possibly loop entered
673 zero_rep(d.elem());
674
675#if ! defined(QDP_NOT_IMPLEMENTED)
676 const int *tab = s.siteTable().slice();
677 for(int n=0; n < s1.size(); ++n)
678 {
679 const OLattice<T>& ss1 = s1[n];
680 for(int j=0; j < s.numSiteTable(); ++j)
681 {
682 int i = tab[j];
683 d.elem() += localNorm2(ss1.elem(i));
684 }
685 }
686#else
687 QDP_error_exit("norm2-Subset not implemented");
688#endif
689
690#if defined(QDP_USE_PROFILING)
691 prof.time += getClockTime();
692 prof.count++;
693 prof.print();
694#endif
695
696 return d;
697}
698
699
701
707template<class T>
708inline typename UnaryReturn<OLattice<T>, FnNorm2>::Type_t
709norm2(const multi1d< OLattice<T> >& s1)
710{
711 return norm2(s1,all);
712}
713
714
715
716//-----------------------------------------------------------------------------
717// Peek and poke at individual sites. This is very architecture specific
718// NOTE: these two routines assume there is no underlying inner grid
719
721
727template<class T1>
728inline typename UnaryReturn<OScalar<T1>, FnPeekSite>::Type_t
729peekSite(const OScalar<T1>& l, const multi1d<int>& coord)
730{
731 return l;
732}
733
735
741template<class RHS, class T1>
742inline OScalar<T1>
743peekSite(const QDPExpr<RHS,OScalar<T1> > & l, const multi1d<int>& coord)
744{
745 // For now, simply evaluate the expression and then call the function
746 typedef OScalar<T1> C1;
747
748 return peekSite(C1(l), coord);
749}
750
751
753
759template<class T1>
761peekSite(const OLattice<T1>& l, const multi1d<int>& coord)
762{
764
765 int i = Layout::linearSiteIndex(coord);
766 int iouter = i >> INNER_LOG;
767 int iinner = i & ((1 << INNER_LOG)-1);
768
769 dest.elem() = getSite(l.elem(iouter), iinner);
770 return dest;
771}
772
774
780template<class RHS, class T1>
781inline OScalar<T1>
782peekSite(const QDPExpr<RHS,OLattice<T1> > & l, const multi1d<int>& coord)
783{
784 // For now, simply evaluate the expression and then call the function
785 typedef OLattice<T1> C1;
786
787 return peekSite(C1(l), coord);
788}
789
790
792
799template<class T1, class T2>
800inline OLattice<T1>&
802{
803 int i = Layout::linearSiteIndex(coord);
804 int iouter = i >> INNER_LOG;
805 int iinner = i & ((1 << INNER_LOG)-1);
806 copy_site(l.elem(iouter), iinner, r.elem());
807 return l;
808}
809
810
812
818template<class T>
819inline void
821 const OLattice<T>& src, const Subset& s)
822{
823 const int *tab = s.siteTable().slice();
824 for(int j=0; j < s.numSiteTable(); ++j)
825 {
826 int i = tab[j];
827 int iouter = i >> INNER_LOG;
828 int iinner = i & ((1 << INNER_LOG)-1);
829
830 dest[i].elem() = getSite(src.elem(iouter),iinner);
831 }
832}
833
834
836
842template<class T>
843inline void
845 const multi1d<OScalar<typename UnaryReturn<T, FnGetSite>::Type_t> >& src,
846 const Subset& s)
847{
848 const int *tab = s.siteTable().slice();
849 for(int j=0; j < s.numSiteTable(); ++j)
850 {
851 int i = tab[j];
852 int iouter = i >> INNER_LOG;
853 int iinner = i & ((1 << INNER_LOG)-1);
854 copy_site(dest.elem(iouter), iinner, src[i].elem());
855 }
856}
857
858
859//-----------------------------------------------------------------------------
860// Map
861//
862// Empty map
863struct FnMap
864{
866};
867
868#if defined(QDP_USE_PROFILING)
869template <>
870struct TagVisitor<FnMap, PrintTag> : public ParenPrinter<FnMap>
871{
872 static void visit(FnMap op, PrintTag t)
873 { t.os_m << "shift"; }
874};
875#endif
876
877
879class Map
880{
881public:
883 Map() {}
884
886 ~Map() {}
887
889 Map(const MapFunc& fn) {make(fn);}
890
892
893 void make(const MapFunc& func);
894
896
905 template<class T1>
908 {
909 OLattice<T1> d;
910
911#if QDP_DEBUG >= 3
912 QDP_info("Map()");
913#endif
914
915 // *** SHOULD IMPROVE THIS - JUST GET IT TO WORK FIRST ***
916 // For now, use the all subset
917#if INNER_LOG == 2
918 const int vvol = Layout::sitesOnNode();
919 for(int i=0; i < vvol; i+= INNER_LEN)
920 {
921 int ii = i >> INNER_LOG;
922 int o0 = goffsets[i+0] >> INNER_LOG;
923 int i0 = goffsets[i+0] & (INNER_LEN - 1);
924
925 int o1 = goffsets[i+1] >> INNER_LOG;
926 int i1 = goffsets[i+1] & (INNER_LEN - 1);
927
928 int o2 = goffsets[i+2] >> INNER_LOG;
929 int i2 = goffsets[i+2] & (INNER_LEN - 1);
930
931 int o3 = goffsets[i+3] >> INNER_LOG;
932 int i3 = goffsets[i+3] & (INNER_LEN - 1);
933
934#if QDP_DEBUG >= 3
935 QDP_info("Map(lattice[%d]=lattice([%d,%d],[%d,%d],[%d,%d],[%d,%d])",
936 ii,o0,i0,o1,i1,o2,i2,o3,i3);
937#endif
938
939 // Gather 4 inner-grid sites together
940 gather_sites(d.elem(ii),
941 l.elem(o0),i0,
942 l.elem(o1),i1,
943 l.elem(o2),i2,
944 l.elem(o3),i3);
945 }
946
947#else
948#error "Map: this inner grid length is not supported - easy to fix"
949#endif
950
951#if QDP_DEBUG >= 3
952 QDP_info("exiting Map()");
953#endif
954
955 return d;
956 }
957
958
959 template<class T1>
962 {
963 return l;
964 }
965
966 template<class RHS, class T1>
969 {
970 // For now, simply evaluate the expression and then do the map
971 typedef OScalar<T1> C1;
972
973// fprintf(stderr,"map(QDPExpr<OScalar>)\n");
974 OScalar<T1> d = this->operator()(C1(l));
975
976 return d;
977 }
978
979 template<class RHS, class T1>
982 {
983 // For now, simply evaluate the expression and then do the map
984 typedef OLattice<T1> C1;
985
986// fprintf(stderr,"map(QDPExpr<OLattice>)\n");
987 OLattice<T1> d = this->operator()(C1(l));
988
989 return d;
990 }
991
992
993public:
995 const multi1d<int>& goffset() const {return goffsets;}
996
997private:
999 Map(const Map&) {}
1000
1002 void operator=(const Map&) {}
1003
1004private:
1009 multi1d<int> goffsets;
1010};
1011
1012
1013//-----------------------------------------------------------------------------
1015class ArrayMap
1016{
1017public:
1020
1023
1025 ArrayMap(const ArrayMapFunc& fn) {make(fn);}
1026
1028
1029 void make(const ArrayMapFunc& func);
1030
1032
1041 template<class T1>
1043 operator()(const OLattice<T1> & l, int dir)
1044 {
1045#if QDP_DEBUG >= 3
1046 QDP_info("ArrayMap(OLattice,%d)",dir);
1047#endif
1048
1049 return mapsa[dir](l);
1050 }
1051
1052 template<class T1>
1054 operator()(const OScalar<T1> & l, int dir)
1055 {
1056#if QDP_DEBUG >= 3
1057 QDP_info("ArrayMap(OScalar,%d)",dir);
1058#endif
1059
1060 return mapsa[dir](l);
1061 }
1062
1063
1064 template<class RHS, class T1>
1066 operator()(const QDPExpr<RHS,OScalar<T1> > & l, int dir)
1067 {
1068// fprintf(stderr,"ArrayMap(QDPExpr<OScalar>,%d)\n",dir);
1069
1070 // For now, simply evaluate the expression and then do the map
1071 return mapsa[dir](l);
1072 }
1073
1074 template<class RHS, class T1>
1076 operator()(const QDPExpr<RHS,OLattice<T1> > & l, int dir)
1077 {
1078// fprintf(stderr,"ArrayMap(QDPExpr<OLattice>,%d)\n",dir);
1079
1080 // For now, simply evaluate the expression and then do the map
1081 return mapsa[dir](l);
1082 }
1083
1084
1085private:
1087 ArrayMap(const ArrayMap&) {}
1088
1090 void operator=(const ArrayMap&) {}
1091
1092private:
1093 multi1d<Map> mapsa;
1094
1095};
1096
1097//-----------------------------------------------------------------------------
1099class BiDirectionalMap
1100{
1101public:
1104
1107
1110
1112
1113 void make(const MapFunc& func);
1114
1116
1125 template<class T1>
1127 operator()(const OLattice<T1> & l, int isign)
1128 {
1129#if QDP_DEBUG >= 3
1130 QDP_info("BiDirectionalMap(OLattice,%d)",isign);
1131#endif
1132
1133 return bimaps[(isign+1)>>1](l);
1134 }
1135
1136
1137 template<class T1>
1139 operator()(const OScalar<T1> & l, int isign)
1140 {
1141#if QDP_DEBUG >= 3
1142 QDP_info("BiDirectionalMap(OScalar,%d)",isign);
1143#endif
1144
1145 return bimaps[(isign+1)>>1](l);
1146 }
1147
1148
1149 template<class RHS, class T1>
1151 operator()(const QDPExpr<RHS,OScalar<T1> > & l, int isign)
1152 {
1153// fprintf(stderr,"BiDirectionalMap(QDPExpr<OScalar>,%d)\n",isign);
1154
1155 // For now, simply evaluate the expression and then do the map
1156 return bimaps[(isign+1)>>1](l);
1157 }
1158
1159 template<class RHS, class T1>
1161 operator()(const QDPExpr<RHS,OLattice<T1> > & l, int isign)
1162 {
1163// fprintf(stderr,"BiDirectionalMap(QDPExpr<OLattice>,%d)\n",isign);
1164
1165 // For now, simply evaluate the expression and then do the map
1166 return bimaps[(isign+1)>>1](l);
1167 }
1168
1169
1170private:
1173
1175 void operator=(const BiDirectionalMap&) {}
1176
1177private:
1178 multi1d<Map> bimaps;
1179
1180};
1181
1182
1183//-----------------------------------------------------------------------------
1185class ArrayBiDirectionalMap
1186{
1187public:
1190
1193
1196
1198
1199 void make(const ArrayMapFunc& func);
1200
1202
1219 template<class T1>
1221 operator()(const OLattice<T1> & l, int isign, int dir)
1222 {
1223#if QDP_DEBUG >= 3
1224 QDP_info("ArrayBiDirectionalMap(OLattice,%d,%d)",isign,dir);
1225#endif
1226
1227 return bimapsa((isign+1)>>1,dir)(l);
1228 }
1229
1230 template<class T1>
1232 operator()(const OScalar<T1> & l, int isign, int dir)
1233 {
1234#if QDP_DEBUG >= 3
1235 QDP_info("ArrayBiDirectionalMap(OScalar,%d,%d)",isign,dir);
1236#endif
1237
1238 return bimapsa((isign+1)>>1,dir)(l);
1239 }
1240
1241
1242 template<class RHS, class T1>
1244 operator()(const QDPExpr<RHS,OScalar<T1> > & l, int isign, int dir)
1245 {
1246// fprintf(stderr,"ArrayBiDirectionalMap(QDPExpr<OScalar>,%d,%d)\n",isign,dir);
1247
1248 // For now, simply evaluate the expression and then do the map
1249 return bimapsa((isign+1)>>1,dir)(l);
1250 }
1251
1252 template<class RHS, class T1>
1254 operator()(const QDPExpr<RHS,OLattice<T1> > & l, int isign, int dir)
1255 {
1256// fprintf(stderr,"ArrayBiDirectionalMap(QDPExpr<OLattice>,%d,%d)\n",isign,dir);
1257
1258 // For now, simply evaluate the expression and then do the map
1259 return bimapsa((isign+1)>>1,dir)(l);
1260 }
1261
1262
1263private:
1266
1268 void operator=(const ArrayBiDirectionalMap&) {}
1269
1270private:
1271 multi2d<Map> bimapsa;
1272
1273};
1274
1275
1276//-----------------------------------------------------------------------------
1277// Input and output of various flavors that are architecture specific
1278
1280multi1d<int> crtesn(int ipos, const multi1d<int>& latt_size);
1281
1283template<class T>
1284XMLWriter& operator<<(XMLWriter& xml, const OLattice<T>& d)
1285{
1286 xml.openTag("OLattice");
1287
1288 XMLWriterAPI::AttributeList alist;
1289
1290 const int iend = Layout::vol();
1291 for(int site=0; site < iend; ++site)
1292 {
1293 multi1d<int> coord = crtesn(site, Layout::lattSize());
1294 std::ostringstream os;
1295 os << coord[0];
1296 for(int i=1; i < coord.size(); ++i)
1297 os << " " << coord[i];
1298
1299 int i = Layout::linearSiteIndex(site);
1300 int outersite = i >> INNER_LOG;
1301 int innersite = i & ((1 << INNER_LOG)-1);
1302
1303 alist.clear();
1304 alist.push_back(XMLWriterAPI::Attribute("site", site));
1305 alist.push_back(XMLWriterAPI::Attribute("coord", os.str()));
1306
1307 xml.openTag("elem", alist);
1308 xml << getSite(d.elem(outersite),innersite);
1309 xml.closeTag();
1310 }
1311
1312 xml.closeTag(); // OLattice
1313
1314 return xml;
1315}
1316
1317
1318
1320
1321template<class T>
1322inline
1323void write(BinaryWriter& bin, const OScalar<T>& d)
1324{
1325 bin.writeArray((const char *)&(d.elem()),
1326 sizeof(typename WordType<T>::Type_t),
1327 sizeof(T) / sizeof(typename WordType<T>::Type_t));
1328}
1329
1331
1332template<class T>
1333void write(BinaryWriter& bin, const OLattice<T>& d)
1334{
1335 const int iend = Layout::vol();
1336 for(int site=0; site < iend; ++site)
1337 {
1338 int i = Layout::linearSiteIndex(site);
1339 int outersite = i >> INNER_LOG;
1340 int innersite = i & ((1 << INNER_LOG)-1);
1341
1342 typedef typename UnaryReturn<T, FnGetSite>::Type_t Site_t;
1343 Site_t this_site = getSite(d.elem(outersite),innersite);
1344
1345 bin.writeArray((const char*)&this_site,
1346 sizeof(typename WordType<Site_t>::Type_t),
1347 sizeof(Site_t) / sizeof(typename WordType<Site_t>::Type_t));
1348 }
1349}
1350
1352
1353template<class T>
1354void write(BinaryWriter& bin, const OLattice<T>& d, const multi1d<int>& coord)
1355{
1356 int i = Layout::linearSiteIndex(coord);
1357 int outersite = i >> INNER_LOG;
1358 int innersite = i & ((1 << INNER_LOG)-1);
1359
1360 typedef typename UnaryReturn<T, FnGetSite>::Type_t Site_t;
1361 Site_t this_site = getSite(d.elem(outersite),innersite);
1362
1363 bin.writeArray((const char*)&this_site,
1364 sizeof(typename WordType<Site_t>::Type_t),
1365 sizeof(Site_t) / sizeof(typename WordType<Site_t>::Type_t));
1366}
1367
1369
1370template<class T>
1371void read(BinaryReader& bin, OScalar<T>& d)
1372{
1373 bin.readArray((char*)&(d.elem()),
1374 sizeof(typename WordType<T>::Type_t),
1375 sizeof(T) / sizeof(typename WordType<T>::Type_t));
1376}
1377
1379
1380template<class T>
1381void read(BinaryReader& bin, OLattice<T>& d)
1382{
1383 const int iend = Layout::vol();
1384 for(int site=0; site < iend; ++site)
1385 {
1386 int i = Layout::linearSiteIndex(site);
1387 int outersite = i >> INNER_LOG;
1388 int innersite = i & ((1 << INNER_LOG)-1);
1389
1390 typedef typename UnaryReturn<T, FnGetSite>::Type_t Site_t;
1391 Site_t this_site;
1392
1393 bin.readArray((char*)&this_site,
1394 sizeof(typename WordType<Site_t>::Type_t),
1395 sizeof(Site_t) / sizeof(typename WordType<Site_t>::Type_t));
1396
1397 copy_site(d.elem(outersite), innersite, this_site);
1398 }
1399}
1400
1402
1403template<class T>
1404void read(BinaryReader& bin, OLattice<T>& d, const multi1d<int>& coord)
1405{
1406 int i = Layout::linearSiteIndex(coord);
1407 int outersite = i >> INNER_LOG;
1408 int innersite = i & ((1 << INNER_LOG)-1);
1409
1410 typedef typename UnaryReturn<T, FnGetSite>::Type_t Site_t;
1411 Site_t this_site;
1412
1413 bin.readArray((char*)&this_site,
1414 sizeof(typename WordType<Site_t>::Type_t),
1415 sizeof(Site_t) / sizeof(typename WordType<Site_t>::Type_t));
1416
1417 copy_site(d.elem(outersite), innersite, this_site);
1418}
1419
1420} // namespace QDP
1421
1422#endif
ForEach< Expr, FTag, CTag >::Type_t forEach(const Expr &e, const FTag &f, const CTag &c)
Definition ForEach.h:87
#define PETE_EMPTY_CONSTRUCTORS(CLASS)
Definition PETE.h:58
ArrayBiDirectional of general permutation map class for communications.
OLattice< T1 > operator()(const OLattice< T1 > &l, int isign, int dir)
Function call operator for a shift.
void make(const ArrayMapFunc &func)
Actual constructor from a function object.
Definition qdp_map.cc:153
ArrayBiDirectionalMap(const ArrayMapFunc &fn)
Constructor from a function object.
ArrayBiDirectionalMap()
Constructor - does nothing really.
OScalar< T1 > operator()(const QDPExpr< RHS, OScalar< T1 > > &l, int isign, int dir)
OScalar< T1 > operator()(const OScalar< T1 > &l, int isign, int dir)
OLattice< T1 > operator()(const QDPExpr< RHS, OLattice< T1 > > &l, int isign, int dir)
ArrayMapFunc.
Definition qdp_map.h:45
Array of general permutation map class for communications.
OLattice< T1 > operator()(const OLattice< T1 > &l, int dir)
Function call operator for a shift.
ArrayMap()
Constructor - does nothing really.
OScalar< T1 > operator()(const OScalar< T1 > &l, int dir)
OLattice< T1 > operator()(const QDPExpr< RHS, OLattice< T1 > > &l, int dir)
ArrayMap(const ArrayMapFunc &fn)
Constructor from a function object.
OScalar< T1 > operator()(const QDPExpr< RHS, OScalar< T1 > > &l, int dir)
void make(const ArrayMapFunc &func)
Actual constructor from a function object.
Definition qdp_map.cc:74
BiDirectional of general permutation map class for communications.
OScalar< T1 > operator()(const QDPExpr< RHS, OScalar< T1 > > &l, int isign)
void make(const MapFunc &func)
Actual constructor from a function object.
Definition qdp_map.cc:114
OLattice< T1 > operator()(const QDPExpr< RHS, OLattice< T1 > > &l, int isign)
BiDirectionalMap(const MapFunc &fn)
Constructor from a function object.
OScalar< T1 > operator()(const OScalar< T1 > &l, int isign)
OLattice< T1 > operator()(const OLattice< T1 > &l, int isign)
Function call operator for a shift.
BiDirectionalMap()
Constructor - does nothing really.
MapFunc.
Definition qdp_map.h:32
General permutation map class for communications.
OLattice< T1 > operator()(const QDPExpr< RHS, OLattice< T1 > > &l)
Map(const MapFunc &fn)
Constructor from a function object.
~Map()
Destructor.
OScalar< T1 > operator()(const QDPExpr< RHS, OScalar< T1 > > &l)
OScalar< T1 > operator()(const OScalar< T1 > &l)
const multi1d< int > & goffset() const
Accessor to offsets.
Map()
Constructor - does nothing really.
OLattice< T1 > operator()(const OLattice< T1 > &l)
Function call operator for a shift.
void make(const MapFunc &func)
Actual constructor from a function object.
Outer grid Lattice type.
Definition qdp_outer.h:264
T & elem(int i)
Definition qdp_outer.h:400
Outer grid Scalar class *‍/.
Definition qdp_outer.h:37
OLattice class narrowed to a subset.
const Subset & subset() const
Expression class for QDP.
Definition qdp_qdpexpr.h:16
OLattice< T1 > & pokeSite(OLattice< T1 > &l, const OScalar< T2 > &r, const multi1d< int > &coord)
Insert site element.
OScalar< T1 > peekSite(const QDPExpr< RHS, OScalar< T1 > > &l, const multi1d< int > &coord)
Extract site element.
UnaryReturn< OScalar< T1 >, FnPeekSite >::Type_t peekSite(const OScalar< T1 > &l, const multi1d< int > &coord)
Extract site element.
UnaryReturn< OLattice< T1 >, FnPeekSite >::Type_t peekSite(const OLattice< T1 > &l, const multi1d< int > &coord)
Extract site element.
OScalar< T1 > peekSite(const QDPExpr< RHS, OLattice< T1 > > &l, const multi1d< int > &coord)
Extract site element.
void QDP_extract(multi1d< OScalar< typename UnaryReturn< T, FnGetSite >::Type_t > > &dest, const OLattice< T > &src, const Subset &s)
Copy data values from field src to array dest.
void QDP_insert(OLattice< T > &dest, const multi1d< OScalar< typename UnaryReturn< T, FnGetSite >::Type_t > > &src, const Subset &s)
Inserts data values from site array src.
Set - collection of subsets controlling which sites are involved in an operation.
Definition qdp_subset.h:96
int numSubsets() const
Return number of subsets.
Definition qdp_subset.h:111
Subsets - controls how lattices are looped.
Definition qdp_subset.h:39
const multi1d< int > & siteTable() const
Definition qdp_subset.h:83
int numSiteTable() const
Definition qdp_subset.h:84
Container for a multi-dimensional 1D array.
Definition qdp_multi.h:25
const T * slice() const
Return ref to a column slice.
Definition qdp_multi.h:225
Container for a multi-dimensional 2D array.
Definition qdp_multi.h:640
#define INNER_LEN
OLattice< PScalar< PSeed< RScalar< INTEGER32 > > > > LatticeSeed
OScalar< PScalar< PSeed< RScalar< INTEGER32 > > > > Seed
#define INNER_LOG
OLattice< PScalar< PScalar< RScalar< INTEGER32 > > > > LatticeInteger
OScalar< T1 > peekSite(const OScalar< T1 > &l, const multi1d< int > &coord)
Extract site element.
UnaryReturn< C, FnNorm2 >::Type_t norm2(const QDPType< T, C > &s1)
OScalar = norm2(trace(adj(source)*source)).
UnaryReturn< C, FnSumMulti >::Type_t sumMulti(const QDPType< T, C > &s1, const Set &ss)
dest = sumMulti(source1,Set)
UnaryReturn< C, FnSum >::Type_t sum(const QDPType< T, C > &s1)
OScalar = sum(source).
void gather_sites(ILattice< T, 2 > &d, const ILattice< T1, 2 > &s0, int i0, const ILattice< T1, 2 > &s1, int i1)
gather several inner sites together
Definition qdp_inner.h:3307
TextWriter & operator<<(TextWriter &txt, const std::string &output)
Definition qdp_io.cc:283
void write(BinaryWriter &bin, const std::string &output)
Definition qdp_io.cc:1204
void read(BinaryReader &bin, std::string &input, size_t maxBytes)
Definition qdp_io.cc:778
void fill_gaussian(IScalar< T > &d, IScalar< T > &r1, IScalar< T > &r2)
dest = gaussian
Definition qdp_inner.h:1861
void copy_site(IScalar< T > &d, int isite, const IScalar< T1 > &s1)
dest [some type] = source [some type]
Definition qdp_inner.h:1682
void zero_rep(IScalar< T > &dest)
dest = 0
Definition qdp_inner.h:1841
void copymask(IScalar< T > &d, const IScalar< T1 > &mask, const IScalar< T > &s1)
dest = (mask) ? s1 : dest
Definition qdp_inner.h:1656
UnaryReturn< IScalar< T >, FnGetSite >::Type_t getSite(const IScalar< T > &s1, int innersite)
Definition qdp_inner.h:1605
void evaluate(OLattice< DCol > &d, const OpAssign &op, const QDPExpr< BinaryNode< OpMultiply, Reference< QDPType< DCol, OLattice< DCol > > >, Reference< QDPType< DCol, OLattice< DCol > > > >, OLattice< DCol > > &rhs, const Subset &s)
void gaussian(OSubScalar< T > &d)
dest = gaussian
Definition qdp_outer.h:1462
void random(OScalar< T > &d)
dest = random
Subset all
Default all subset.
Definition qdp_subset.cc:16
Layout namespace holding info on problem size and machine info.
Definition qdp_layout.cc:16
int linearSiteIndex(int site)
The linearized site index for the corresponding lexicographic site.
int sitesOnNode()
Subgrid lattice volume.
int outerSitesOnNode()
Subgrid lattice volume.
const multi1d< int > & lattSize()
Virtual grid (problem grid) lattice size.
LatticeInteger latticeCoordinate(int mu)
coord[mu] <- mu : fill with lattice coord in mu direction
int vol()
Total lattice volume.
void broadcast(T &dest)
Broadcast from primary node to all other nodes.
void globalSum(T &dest)
Sum across all nodes.
void globalSumArray(unsigned int *dest, int len)
Wrapper to get a functional unsigned global sum.
Seed ran_mult_n
RNG multiplier raised to the volume+1.
Definition qdp_random.cc:24
Seed ran_mult
RNG multiplier.
Definition qdp_random.cc:22
LatticeSeed * lattice_ran_mult
The lattice of skewed RNG multipliers.
Definition qdp_random.cc:26
Seed ran_seed
Global (current) seed.
Definition qdp_random.cc:20
Yet another random number generator.
ForEach< Expr, FTag, CTag >::Type_t forEach(const Expr &e, const FTag &f, const CTag &c)
Definition qdp.h:88
void QDP_error_exit(const char *format,...)
Simple error display and abort routine.
Definition qdp_util.cc:93
QDPTime_t getClockTime()
Get the wallclock time.
multi1d< int > crtesn(int ipos, const multi1d< int > &latt_size)
Decompose a lexicographic site into coordinates.
void fill_random(float &d, T1 &seed, T2 &skewed_seed, const T1 &seed_mult)
dest = random
Definition qdp_random.h:54
MakeReturn< UnaryNode< FnLocalNorm2, typenameCreateLeaf< QDPExpr< T1, C1 > >::Leaf_t >, typenameUnaryReturn< C1, FnLocalNorm2 >::Type_t >::Expression_t localNorm2(const QDPExpr< T1, C1 > &l)
Definition qdp.h:4998
int QDP_info(const char *format,...)
Simple information display routine.
Definition qdp_util.cc:43
#define QDP_CONST
Definition qdp.h:63