QDP++
qdp_scalarsite_sse.cc
Go to the documentation of this file.
1
7
8
9#include "qdp.h"
10
11
12// These SSE asm instructions are only supported under GCC/G++
13#if defined(__GNUC__)
14#include "qdp_sse_intrin.h"
15namespace QDP {
16
17
18
19#if 1
20//-------------------------------------------------------------------
21// Specialization to optimize the case
22// LatticeColorMatrix[ Subset] = LatticeColorMatrix * LatticeColorMatrix
23template<>
24void evaluate(OLattice< TCol >& d,
25 const OpAssign& op,
26 const QDPExpr<BinaryNode<OpMultiply,
27 Reference<QDPType< TCol, OLattice< TCol > > >,
28 Reference<QDPType< TCol, OLattice< TCol > > > >,
29 OLattice< TCol > >& rhs,
30 const Subset& s)
31{
32// cout << "call single site QDP_M_eq_M_times_M" << endl;
33
34 typedef OLattice< TCol > C;
35
36 const C& l = static_cast<const C&>(rhs.expression().left());
37 const C& r = static_cast<const C&>(rhs.expression().right());
38
39 if( s.hasOrderedRep() ) {
40 for(int i=s.start(); i <= s.end(); i++) {
41
42 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
43 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
44 su3_matrixf *dm = (su3_matrixf *)&(d.elem(i).elem().elem(0,0).real());
45
46 intrin_sse_mult_su3_nn(lm, rm, dm);
47
48 }
49 }
50 else {
51 const int *tab = s.siteTable().slice();
52 for(int j=0; j < s.numSiteTable(); ++j) {
53 int i = tab[j];
54 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
55 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
56 su3_matrixf *dm = (su3_matrixf *)&(d.elem(i).elem().elem(0,0).real());
57
58 intrin_sse_mult_su3_nn(lm, rm, dm);
59
60 }
61 }
62}
63
64
65// Specialization to optimize the case
66// LatticeColorMatrix[Subset] = adj(LatticeColorMatrix) * LatticeColorMatrix
67template<>
68void evaluate(OLattice< TCol >& d,
69 const OpAssign& op,
70 const QDPExpr<BinaryNode<OpAdjMultiply,
71 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > >,
72 Reference<QDPType< TCol, OLattice< TCol > > > >,
73 OLattice< TCol > >& rhs,
74 const Subset& s)
75{
76// cout << "call single site QDP_M_eq_aM_times_M" << endl;
77
78 typedef OLattice< TCol > C;
79
80 const C& l = static_cast<const C&>(rhs.expression().left().child());
81 const C& r = static_cast<const C&>(rhs.expression().right());
82
83 if( s.hasOrderedRep() ) {
84 for(int i=s.start(); i <= s.end(); i++) {
85 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
86 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
87 su3_matrixf *dm = (su3_matrixf *)&(d.elem(i).elem().elem(0,0).real());
88
89 intrin_sse_mult_su3_an(lm, rm, dm);
90
91 }
92 }
93 else {
94 const int *tab = s.siteTable().slice();
95 for(int j=0; j < s.numSiteTable(); ++j) {
96
97 int i = tab[j];
98 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
99 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
100 su3_matrixf *dm = (su3_matrixf *)&(d.elem(i).elem().elem(0,0).real());
101
102 intrin_sse_mult_su3_an(lm, rm, dm);
103
104 }
105 }
106
107}
108
109
110// Specialization to optimize the case
111// LatticeColorMatrix[Subset] = LatticeColorMatrix * adj(LatticeColorMatrix)
112template<>
113void evaluate(OLattice< TCol >& d,
114 const OpAssign& op,
115 const QDPExpr<BinaryNode<OpMultiplyAdj,
116 Reference<QDPType< TCol, OLattice< TCol > > >,
117 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > > >,
118 OLattice< TCol > >& rhs,
119 const Subset& s)
120{
121// cout << "call single site QDP_M_eq_M_times_aM" << endl;
122
123 typedef OLattice< TCol > C;
124
125 const C& l = static_cast<const C&>(rhs.expression().left());
126 const C& r = static_cast<const C&>(rhs.expression().right().child());
127
128 if( s.hasOrderedRep() ) {
129 for(int i=s.start(); i <= s.end(); i++) {
130 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
131 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
132 su3_matrixf *dm = (su3_matrixf *)&(d.elem(i).elem().elem(0,0).real());
133
134 intrin_sse_mult_su3_na(lm, rm, dm);
135
136 }
137 }
138 else {
139
140 const int *tab = s.siteTable().slice();
141 for(int j=0; j < s.numSiteTable(); ++j) {
142 int i = tab[j];
143 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
144 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
145 su3_matrixf *dm = (su3_matrixf *)&(d.elem(i).elem().elem(0,0).real());
146
147 intrin_sse_mult_su3_na(lm, rm, dm);
148
149 }
150 }
151}
152
153
154// Specialization to optimize the case
155// LatticeColorMatrix[Subset] = adj(LatticeColorMatrix) * adj(LatticeColorMatrix)
156template<>
157void evaluate(OLattice< TCol >& d,
158 const OpAssign& op,
159 const QDPExpr<BinaryNode<OpAdjMultiplyAdj,
160 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > >,
161 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > > >,
162 OLattice< TCol > >& rhs,
163 const Subset& s)
164{
165// cout << "call single site QDP_M_eq_Ma_times_Ma" << endl;
166
167 typedef OLattice< TCol > C;
168
169 const C& l = static_cast<const C&>(rhs.expression().left().child());
170 const C& r = static_cast<const C&>(rhs.expression().right().child());
171
172 PColorMatrix<RComplexFloat,3> tmp;
173
174 if( s.hasOrderedRep() ) {
175 for(int i=s.start(); i <= s.end(); i++) {
176 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
177 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
178 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
179
180 intrin_sse_mult_su3_nn(rm, lm, tmpm);
181
182 // Take the adj(r*l) = adj(l)*adj(r)
183 d.elem(i).elem().elem(0,0).real() = tmp.elem(0,0).real();
184 d.elem(i).elem().elem(0,0).imag() = -tmp.elem(0,0).imag();
185 d.elem(i).elem().elem(0,1).real() = tmp.elem(1,0).real();
186 d.elem(i).elem().elem(0,1).imag() = -tmp.elem(1,0).imag();
187 d.elem(i).elem().elem(0,2).real() = tmp.elem(2,0).real();
188 d.elem(i).elem().elem(0,2).imag() = -tmp.elem(2,0).imag();
189
190 d.elem(i).elem().elem(1,0).real() = tmp.elem(0,1).real();
191 d.elem(i).elem().elem(1,0).imag() = -tmp.elem(0,1).imag();
192 d.elem(i).elem().elem(1,1).real() = tmp.elem(1,1).real();
193 d.elem(i).elem().elem(1,1).imag() = -tmp.elem(1,1).imag();
194 d.elem(i).elem().elem(1,2).real() = tmp.elem(2,1).real();
195 d.elem(i).elem().elem(1,2).imag() = -tmp.elem(2,1).imag();
196
197 d.elem(i).elem().elem(2,0).real() = tmp.elem(0,2).real();
198 d.elem(i).elem().elem(2,0).imag() = -tmp.elem(0,2).imag();
199 d.elem(i).elem().elem(2,1).real() = tmp.elem(1,2).real();
200 d.elem(i).elem().elem(2,1).imag() = -tmp.elem(1,2).imag();
201 d.elem(i).elem().elem(2,2).real() = tmp.elem(2,2).real();
202 d.elem(i).elem().elem(2,2).imag() = -tmp.elem(2,2).imag();
203 }
204 }
205 else {
206 const int *tab = s.siteTable().slice();
207 for(int j=0; j < s.numSiteTable(); ++j) {
208 int i = tab[j];
209 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
210 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
211 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
212
213 intrin_sse_mult_su3_nn(rm, lm, tmpm);
214
215 // Take the adj(r*l) = adj(l)*adj(r)
216 d.elem(i).elem().elem(0,0).real() = tmp.elem(0,0).real();
217 d.elem(i).elem().elem(0,0).imag() = -tmp.elem(0,0).imag();
218 d.elem(i).elem().elem(0,1).real() = tmp.elem(1,0).real();
219 d.elem(i).elem().elem(0,1).imag() = -tmp.elem(1,0).imag();
220 d.elem(i).elem().elem(0,2).real() = tmp.elem(2,0).real();
221 d.elem(i).elem().elem(0,2).imag() = -tmp.elem(2,0).imag();
222
223 d.elem(i).elem().elem(1,0).real() = tmp.elem(0,1).real();
224 d.elem(i).elem().elem(1,0).imag() = -tmp.elem(0,1).imag();
225 d.elem(i).elem().elem(1,1).real() = tmp.elem(1,1).real();
226 d.elem(i).elem().elem(1,1).imag() = -tmp.elem(1,1).imag();
227 d.elem(i).elem().elem(1,2).real() = tmp.elem(2,1).real();
228 d.elem(i).elem().elem(1,2).imag() = -tmp.elem(2,1).imag();
229
230 d.elem(i).elem().elem(2,0).real() = tmp.elem(0,2).real();
231 d.elem(i).elem().elem(2,0).imag() = -tmp.elem(0,2).imag();
232 d.elem(i).elem().elem(2,1).real() = tmp.elem(1,2).real();
233 d.elem(i).elem().elem(2,1).imag() = -tmp.elem(1,2).imag();
234 d.elem(i).elem().elem(2,2).real() = tmp.elem(2,2).real();
235 d.elem(i).elem().elem(2,2).imag() = -tmp.elem(2,2).imag();
236 }
237 }
238}
239
240//-------------------------------------------------------------------
241
242// Specialization to optimize the case
243// LatticeColorMatrix[Subset] += LatticeColorMatrix * LatticeColorMatrix
244template<>
245void evaluate(OLattice< TCol >& d,
246 const OpAddAssign& op,
247 const QDPExpr<BinaryNode<OpMultiply,
248 Reference<QDPType< TCol, OLattice< TCol > > >,
249 Reference<QDPType< TCol, OLattice< TCol > > > >,
250 OLattice< TCol > >& rhs,
251 const Subset& s)
252{
253// cout << "call single site QDP_M_peq_M_times_M" << endl;
254
255 typedef OLattice< TCol > C;
256
257 const C& l = static_cast<const C&>(rhs.expression().left());
258 const C& r = static_cast<const C&>(rhs.expression().right());
259
260 PColorMatrix<RComplexFloat,3> tmp;
261
262 if( s.hasOrderedRep() ) {
263 for(int i=s.start(); i <= s.end(); i++) {
264 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
265 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
266 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
267
268 intrin_sse_mult_su3_nn(lm, rm, tmpm);
269
270 d.elem(i).elem().elem(0,0).real() += tmp.elem(0,0).real();
271 d.elem(i).elem().elem(0,0).imag() += tmp.elem(0,0).imag();
272 d.elem(i).elem().elem(0,1).real() += tmp.elem(0,1).real();
273 d.elem(i).elem().elem(0,1).imag() += tmp.elem(0,1).imag();
274 d.elem(i).elem().elem(0,2).real() += tmp.elem(0,2).real();
275 d.elem(i).elem().elem(0,2).imag() += tmp.elem(0,2).imag();
276
277 d.elem(i).elem().elem(1,0).real() += tmp.elem(1,0).real();
278 d.elem(i).elem().elem(1,0).imag() += tmp.elem(1,0).imag();
279 d.elem(i).elem().elem(1,1).real() += tmp.elem(1,1).real();
280 d.elem(i).elem().elem(1,1).imag() += tmp.elem(1,1).imag();
281 d.elem(i).elem().elem(1,2).real() += tmp.elem(1,2).real();
282 d.elem(i).elem().elem(1,2).imag() += tmp.elem(1,2).imag();
283
284 d.elem(i).elem().elem(2,0).real() += tmp.elem(2,0).real();
285 d.elem(i).elem().elem(2,0).imag() += tmp.elem(2,0).imag();
286 d.elem(i).elem().elem(2,1).real() += tmp.elem(2,1).real();
287 d.elem(i).elem().elem(2,1).imag() += tmp.elem(2,1).imag();
288 d.elem(i).elem().elem(2,2).real() += tmp.elem(2,2).real();
289 d.elem(i).elem().elem(2,2).imag() += tmp.elem(2,2).imag();
290 }
291 }
292 else {
293 const int *tab = s.siteTable().slice();
294 for(int j=0; j < s.numSiteTable(); ++j) {
295 int i = tab[j];
296 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
297 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
298 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
299
300 intrin_sse_mult_su3_nn(lm, rm, tmpm);
301
302 d.elem(i).elem().elem(0,0).real() += tmp.elem(0,0).real();
303 d.elem(i).elem().elem(0,0).imag() += tmp.elem(0,0).imag();
304 d.elem(i).elem().elem(0,1).real() += tmp.elem(0,1).real();
305 d.elem(i).elem().elem(0,1).imag() += tmp.elem(0,1).imag();
306 d.elem(i).elem().elem(0,2).real() += tmp.elem(0,2).real();
307 d.elem(i).elem().elem(0,2).imag() += tmp.elem(0,2).imag();
308
309 d.elem(i).elem().elem(1,0).real() += tmp.elem(1,0).real();
310 d.elem(i).elem().elem(1,0).imag() += tmp.elem(1,0).imag();
311 d.elem(i).elem().elem(1,1).real() += tmp.elem(1,1).real();
312 d.elem(i).elem().elem(1,1).imag() += tmp.elem(1,1).imag();
313 d.elem(i).elem().elem(1,2).real() += tmp.elem(1,2).real();
314 d.elem(i).elem().elem(1,2).imag() += tmp.elem(1,2).imag();
315
316 d.elem(i).elem().elem(2,0).real() += tmp.elem(2,0).real();
317 d.elem(i).elem().elem(2,0).imag() += tmp.elem(2,0).imag();
318 d.elem(i).elem().elem(2,1).real() += tmp.elem(2,1).real();
319 d.elem(i).elem().elem(2,1).imag() += tmp.elem(2,1).imag();
320 d.elem(i).elem().elem(2,2).real() += tmp.elem(2,2).real();
321 d.elem(i).elem().elem(2,2).imag() += tmp.elem(2,2).imag();
322 }
323 }
324}
325
326
327// Specialization to optimize the case
328// LatticeColorMatrix[Subset] += adj(LatticeColorMatrix) * LatticeColorMatrix
329template<>
330void evaluate(OLattice< TCol >& d,
331 const OpAddAssign& op,
332 const QDPExpr<BinaryNode<OpAdjMultiply,
333 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > >,
334 Reference<QDPType< TCol, OLattice< TCol > > > >,
335 OLattice< TCol > >& rhs,
336 const Subset& s)
337{
338// cout << "call single site QDP_M_peq_aM_times_M" << endl;
339
340 typedef OLattice< TCol > C;
341
342 const C& l = static_cast<const C&>(rhs.expression().left().child());
343 const C& r = static_cast<const C&>(rhs.expression().right());
344
345 PColorMatrix<RComplexFloat,3> tmp;
346
347 if( s.hasOrderedRep() ) {
348 for(int i=s.start(); i <= s.end(); i++) {
349 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
350 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
351 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
352
353 intrin_sse_mult_su3_an(lm, rm, tmpm);
354
355 d.elem(i).elem().elem(0,0).real() += tmp.elem(0,0).real();
356 d.elem(i).elem().elem(0,0).imag() += tmp.elem(0,0).imag();
357 d.elem(i).elem().elem(0,1).real() += tmp.elem(0,1).real();
358 d.elem(i).elem().elem(0,1).imag() += tmp.elem(0,1).imag();
359 d.elem(i).elem().elem(0,2).real() += tmp.elem(0,2).real();
360 d.elem(i).elem().elem(0,2).imag() += tmp.elem(0,2).imag();
361
362 d.elem(i).elem().elem(1,0).real() += tmp.elem(1,0).real();
363 d.elem(i).elem().elem(1,0).imag() += tmp.elem(1,0).imag();
364 d.elem(i).elem().elem(1,1).real() += tmp.elem(1,1).real();
365 d.elem(i).elem().elem(1,1).imag() += tmp.elem(1,1).imag();
366 d.elem(i).elem().elem(1,2).real() += tmp.elem(1,2).real();
367 d.elem(i).elem().elem(1,2).imag() += tmp.elem(1,2).imag();
368
369 d.elem(i).elem().elem(2,0).real() += tmp.elem(2,0).real();
370 d.elem(i).elem().elem(2,0).imag() += tmp.elem(2,0).imag();
371 d.elem(i).elem().elem(2,1).real() += tmp.elem(2,1).real();
372 d.elem(i).elem().elem(2,1).imag() += tmp.elem(2,1).imag();
373 d.elem(i).elem().elem(2,2).real() += tmp.elem(2,2).real();
374 d.elem(i).elem().elem(2,2).imag() += tmp.elem(2,2).imag();
375 }
376 }
377 else {
378 const int *tab = s.siteTable().slice();
379 for(int j=0; j < s.numSiteTable(); ++j) {
380 int i = tab[j];
381 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
382 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
383 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
384
385 intrin_sse_mult_su3_an(lm, rm, tmpm);
386
387 d.elem(i).elem().elem(0,0).real() += tmp.elem(0,0).real();
388 d.elem(i).elem().elem(0,0).imag() += tmp.elem(0,0).imag();
389 d.elem(i).elem().elem(0,1).real() += tmp.elem(0,1).real();
390 d.elem(i).elem().elem(0,1).imag() += tmp.elem(0,1).imag();
391 d.elem(i).elem().elem(0,2).real() += tmp.elem(0,2).real();
392 d.elem(i).elem().elem(0,2).imag() += tmp.elem(0,2).imag();
393
394 d.elem(i).elem().elem(1,0).real() += tmp.elem(1,0).real();
395 d.elem(i).elem().elem(1,0).imag() += tmp.elem(1,0).imag();
396 d.elem(i).elem().elem(1,1).real() += tmp.elem(1,1).real();
397 d.elem(i).elem().elem(1,1).imag() += tmp.elem(1,1).imag();
398 d.elem(i).elem().elem(1,2).real() += tmp.elem(1,2).real();
399 d.elem(i).elem().elem(1,2).imag() += tmp.elem(1,2).imag();
400
401 d.elem(i).elem().elem(2,0).real() += tmp.elem(2,0).real();
402 d.elem(i).elem().elem(2,0).imag() += tmp.elem(2,0).imag();
403 d.elem(i).elem().elem(2,1).real() += tmp.elem(2,1).real();
404 d.elem(i).elem().elem(2,1).imag() += tmp.elem(2,1).imag();
405 d.elem(i).elem().elem(2,2).real() += tmp.elem(2,2).real();
406 d.elem(i).elem().elem(2,2).imag() += tmp.elem(2,2).imag();
407 }
408 }
409}
410
411
412// Specialization to optimize the case
413// LatticeColorMatrix[Subset] += LatticeColorMatrix * adj(LatticeColorMatrix)
414template<>
415void evaluate(OLattice< TCol >& d,
416 const OpAddAssign& op,
417 const QDPExpr<BinaryNode<OpMultiplyAdj,
418 Reference<QDPType< TCol, OLattice< TCol > > >,
419 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > > >,
420 OLattice< TCol > >& rhs,
421 const Subset& s)
422{
423// cout << "call single site QDP_M_peq_M_times_aM" << endl;
424
425 typedef OLattice< TCol > C;
426
427 const C& l = static_cast<const C&>(rhs.expression().left());
428 const C& r = static_cast<const C&>(rhs.expression().right().child());
429
430 PColorMatrix<RComplexFloat,3> tmp;
431
432 if( s.hasOrderedRep() ) {
433 for(int i=s.start(); i <= s.end(); i++) {
434 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
435 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
436 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
437
438 intrin_sse_mult_su3_na(lm, rm, tmpm);
439
440
441 d.elem(i).elem().elem(0,0).real() += tmp.elem(0,0).real();
442 d.elem(i).elem().elem(0,0).imag() += tmp.elem(0,0).imag();
443 d.elem(i).elem().elem(0,1).real() += tmp.elem(0,1).real();
444 d.elem(i).elem().elem(0,1).imag() += tmp.elem(0,1).imag();
445 d.elem(i).elem().elem(0,2).real() += tmp.elem(0,2).real();
446 d.elem(i).elem().elem(0,2).imag() += tmp.elem(0,2).imag();
447
448 d.elem(i).elem().elem(1,0).real() += tmp.elem(1,0).real();
449 d.elem(i).elem().elem(1,0).imag() += tmp.elem(1,0).imag();
450 d.elem(i).elem().elem(1,1).real() += tmp.elem(1,1).real();
451 d.elem(i).elem().elem(1,1).imag() += tmp.elem(1,1).imag();
452 d.elem(i).elem().elem(1,2).real() += tmp.elem(1,2).real();
453 d.elem(i).elem().elem(1,2).imag() += tmp.elem(1,2).imag();
454
455 d.elem(i).elem().elem(2,0).real() += tmp.elem(2,0).real();
456 d.elem(i).elem().elem(2,0).imag() += tmp.elem(2,0).imag();
457 d.elem(i).elem().elem(2,1).real() += tmp.elem(2,1).real();
458 d.elem(i).elem().elem(2,1).imag() += tmp.elem(2,1).imag();
459 d.elem(i).elem().elem(2,2).real() += tmp.elem(2,2).real();
460 d.elem(i).elem().elem(2,2).imag() += tmp.elem(2,2).imag();
461 }
462 }
463 else {
464
465 const int *tab = s.siteTable().slice();
466 for(int j=0; j < s.numSiteTable(); ++j) {
467 int i = tab[j];
468 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
469 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
470 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
471
472 intrin_sse_mult_su3_na(lm, rm, tmpm);
473
474 d.elem(i).elem().elem(0,0).real() += tmp.elem(0,0).real();
475 d.elem(i).elem().elem(0,0).imag() += tmp.elem(0,0).imag();
476 d.elem(i).elem().elem(0,1).real() += tmp.elem(0,1).real();
477 d.elem(i).elem().elem(0,1).imag() += tmp.elem(0,1).imag();
478 d.elem(i).elem().elem(0,2).real() += tmp.elem(0,2).real();
479 d.elem(i).elem().elem(0,2).imag() += tmp.elem(0,2).imag();
480
481 d.elem(i).elem().elem(1,0).real() += tmp.elem(1,0).real();
482 d.elem(i).elem().elem(1,0).imag() += tmp.elem(1,0).imag();
483 d.elem(i).elem().elem(1,1).real() += tmp.elem(1,1).real();
484 d.elem(i).elem().elem(1,1).imag() += tmp.elem(1,1).imag();
485 d.elem(i).elem().elem(1,2).real() += tmp.elem(1,2).real();
486 d.elem(i).elem().elem(1,2).imag() += tmp.elem(1,2).imag();
487
488 d.elem(i).elem().elem(2,0).real() += tmp.elem(2,0).real();
489 d.elem(i).elem().elem(2,0).imag() += tmp.elem(2,0).imag();
490 d.elem(i).elem().elem(2,1).real() += tmp.elem(2,1).real();
491 d.elem(i).elem().elem(2,1).imag() += tmp.elem(2,1).imag();
492 d.elem(i).elem().elem(2,2).real() += tmp.elem(2,2).real();
493 d.elem(i).elem().elem(2,2).imag() += tmp.elem(2,2).imag();
494 }
495 }
496}
497
498
499// Specialization to optimize the case
500// LatticeColorMatrix[Subset] += adj(LatticeColorMatrix) * adj(LatticeColorMatrix)
501template<>
502void evaluate(OLattice< TCol >& d,
503 const OpAddAssign& op,
504 const QDPExpr<BinaryNode<OpAdjMultiplyAdj,
505 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > >,
506 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > > >,
507 OLattice< TCol > >& rhs,
508 const Subset& s)
509{
510// cout << "call single site QDP_M_peq_Ma_times_Ma" << endl;
511
512 typedef OLattice< TCol > C;
513
514 const C& l = static_cast<const C&>(rhs.expression().left().child());
515 const C& r = static_cast<const C&>(rhs.expression().right().child());
516
517 PColorMatrix<RComplexFloat,3> tmp;
518
519 if( s.hasOrderedRep() ) {
520 for(int i=s.start(); i <= s.end(); i++) {
521
522 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
523 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
524 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
525
526 intrin_sse_mult_su3_nn(rm, lm, tmpm);
527
528 // Take the adj(r*l) = adj(l)*adj(r)
529 d.elem(i).elem().elem(0,0).real() += tmp.elem(0,0).real();
530 d.elem(i).elem().elem(0,0).imag() -= tmp.elem(0,0).imag();
531 d.elem(i).elem().elem(0,1).real() += tmp.elem(1,0).real();
532 d.elem(i).elem().elem(0,1).imag() -= tmp.elem(1,0).imag();
533 d.elem(i).elem().elem(0,2).real() += tmp.elem(2,0).real();
534 d.elem(i).elem().elem(0,2).imag() -= tmp.elem(2,0).imag();
535
536 d.elem(i).elem().elem(1,0).real() += tmp.elem(0,1).real();
537 d.elem(i).elem().elem(1,0).imag() -= tmp.elem(0,1).imag();
538 d.elem(i).elem().elem(1,1).real() += tmp.elem(1,1).real();
539 d.elem(i).elem().elem(1,1).imag() -= tmp.elem(1,1).imag();
540 d.elem(i).elem().elem(1,2).real() += tmp.elem(2,1).real();
541 d.elem(i).elem().elem(1,2).imag() -= tmp.elem(2,1).imag();
542
543 d.elem(i).elem().elem(2,0).real() += tmp.elem(0,2).real();
544 d.elem(i).elem().elem(2,0).imag() -= tmp.elem(0,2).imag();
545 d.elem(i).elem().elem(2,1).real() += tmp.elem(1,2).real();
546 d.elem(i).elem().elem(2,1).imag() -= tmp.elem(1,2).imag();
547 d.elem(i).elem().elem(2,2).real() += tmp.elem(2,2).real();
548 d.elem(i).elem().elem(2,2).imag() -= tmp.elem(2,2).imag();
549
550 }
551 }
552 else {
553
554 const int *tab = s.siteTable().slice();
555 for(int j=0; j < s.numSiteTable(); ++j) {
556 int i = tab[j];
557 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
558 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
559 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
560
561 intrin_sse_mult_su3_nn(rm, lm, tmpm);
562
563 // Take the adj(r*l) = adj(l)*adj(r)
564 d.elem(i).elem().elem(0,0).real() += tmp.elem(0,0).real();
565 d.elem(i).elem().elem(0,0).imag() -= tmp.elem(0,0).imag();
566 d.elem(i).elem().elem(0,1).real() += tmp.elem(1,0).real();
567 d.elem(i).elem().elem(0,1).imag() -= tmp.elem(1,0).imag();
568 d.elem(i).elem().elem(0,2).real() += tmp.elem(2,0).real();
569 d.elem(i).elem().elem(0,2).imag() -= tmp.elem(2,0).imag();
570
571 d.elem(i).elem().elem(1,0).real() += tmp.elem(0,1).real();
572 d.elem(i).elem().elem(1,0).imag() -= tmp.elem(0,1).imag();
573 d.elem(i).elem().elem(1,1).real() += tmp.elem(1,1).real();
574 d.elem(i).elem().elem(1,1).imag() -= tmp.elem(1,1).imag();
575 d.elem(i).elem().elem(1,2).real() += tmp.elem(2,1).real();
576 d.elem(i).elem().elem(1,2).imag() -= tmp.elem(2,1).imag();
577
578 d.elem(i).elem().elem(2,0).real() += tmp.elem(0,2).real();
579 d.elem(i).elem().elem(2,0).imag() -= tmp.elem(0,2).imag();
580 d.elem(i).elem().elem(2,1).real() += tmp.elem(1,2).real();
581 d.elem(i).elem().elem(2,1).imag() -= tmp.elem(1,2).imag();
582 d.elem(i).elem().elem(2,2).real() += tmp.elem(2,2).real();
583 d.elem(i).elem().elem(2,2).imag() -= tmp.elem(2,2).imag();
584 }
585 }
586}
587
588//-------------------------------------------------------------------
589// Specialization to optimize the case
590// LatticeColorMatrix = LatticeColorMatrix
591template<>
592void evaluate(OLattice< TCol >& d,
593 const OpAssign& op,
594 const QDPExpr<
595 UnaryNode<OpIdentity, Reference< QDPType< TCol, OLattice< TCol > > > >,
596 OLattice< TCol > >& rhs,
597 const Subset& s)
598{
599 typedef OLattice<TCol> C;
600 const C& l = static_cast<const C&>(rhs.expression().child());
601
602
603 if( s.hasOrderedRep() ) {
604
605 const int start = s.start();
606 const int end = s.end();
607
608 REAL32* d_ptr =&(d.elem(start).elem().elem(0,0).real());
609 const REAL32* r_ptr =&(l.elem(start).elem().elem(0,0).real());
610
611 const unsigned int total_reals = (end-start+1)*3*3*2;
612 const unsigned int total_v4sf = total_reals/4;
613 const unsigned int remainder = total_reals%4;
614
615 float* d_ptr_v4sf = (float *)d_ptr;
616 float* r_ptr_v4sf = (float *)r_ptr;
617
618
619 for(unsigned int i = 0 ; i < total_v4sf; i++, d_ptr_v4sf +=4, r_ptr_v4sf+=4 ) {
620 _mm_store_ps( d_ptr_v4sf, _mm_load_ps(r_ptr_v4sf));
621 }
622
623
624 r_ptr = (REAL32 *)r_ptr_v4sf;
625 d_ptr = (REAL32 *)d_ptr_v4sf;
626 for(unsigned int i=0; i < remainder; i++, r_ptr++, d_ptr++) {
627 *d_ptr = *r_ptr;
628 }
629 }
630 else {
631 // Unordered case
632 const int* tab = s.siteTable().slice();
633
634 // Loop through the sites
635 for(int j=0; j < s.numSiteTable(); j++) {
636 int i = tab[j];
637
638 // Do the copy in the dumb way -- this could become quite complex
639 // Depending on whether the individual matrices are aligned or not.
640 d.elem(i).elem() = l.elem(i).elem();
641
642 }
643
644 }
645}
646
647//-------------------------------------------------------------------
648// Specialization to optimize the case
649// LatticeColorMatrix += LatticeColorMatrix
650template<>
651void evaluate(OLattice< TCol >& d,
652 const OpAddAssign& op,
653 const QDPExpr<
654 UnaryNode<OpIdentity, Reference< QDPType< TCol, OLattice< TCol > > > >,
655 OLattice< TCol > >& rhs,
656 const Subset& s)
657{
658 typedef OLattice<TCol> C;
659 const C& l = static_cast<const C&>(rhs.expression().child());
660
661 if( s.hasOrderedRep() ) {
662 const int start = s.start();
663 const int end = s.end();
664
665 REAL32* d_ptr =&(d.elem(start).elem().elem(0,0).real());
666 const REAL32* r_ptr =&(l.elem(start).elem().elem(0,0).real());
667
668 const unsigned int total_reals = (end-start+1)*3*3*2;
669 const unsigned int total_v4sf = total_reals/4;
670 const unsigned int remainder = total_reals%4;
671
672 float* d_ptr_v4sf = (float *)d_ptr;
673 float* r_ptr_v4sf = (float *)r_ptr;
674
675
676 for(unsigned int i = 0 ; i < total_v4sf; i++, d_ptr_v4sf +=4, r_ptr_v4sf+=4 ) {
677 _mm_store_ps( d_ptr_v4sf, _mm_add_ps( _mm_load_ps(d_ptr_v4sf),
678 _mm_load_ps(r_ptr_v4sf) ) );
679 }
680
681
682 r_ptr = (REAL32 *)r_ptr_v4sf;
683 d_ptr = (REAL32 *)d_ptr_v4sf;
684 for(unsigned int i=0; i < remainder; i++, r_ptr++, d_ptr++) {
685 *d_ptr += *r_ptr;
686 }
687 }
688 else {
689 // Unordered case
690 const int* tab = s.siteTable().slice();
691
692 // Loop through the sites
693 for(int j=0; j < s.numSiteTable(); j++) {
694 int i = tab[j];
695
696 // Do the copy in the dumb way -- this could become quite complex
697 // Depending on whether the individual matrices are aligned or not.
698 d.elem(i).elem() += l.elem(i).elem();
699
700 }
701
702 }
703}
704
705//-------------------------------------------------------------------
706// Specialization to optimize the case
707// LatticeColorMatrix -= LatticeColorMatrix
708template<>
709void evaluate(OLattice< TCol >& d,
710 const OpSubtractAssign& op,
711 const QDPExpr<
712 UnaryNode<OpIdentity, Reference< QDPType< TCol, OLattice< TCol > > > >,
713 OLattice< TCol > >& rhs,
714 const Subset& s)
715{
716 typedef OLattice<TCol> C;
717 const C& l = static_cast<const C&>(rhs.expression().child());
718 if (s.hasOrderedRep()) {
719 const int start = s.start();
720 const int end = s.end();
721
722 REAL32* d_ptr =&(d.elem(start).elem().elem(0,0).real());
723 const REAL32* r_ptr =&(l.elem(start).elem().elem(0,0).real());
724
725 const unsigned int total_reals = (end-start+1)*3*3*2;
726 const unsigned int total_v4sf = total_reals/4;
727 const unsigned int remainder = total_reals%4;
728
729 float* d_ptr_v4sf = (float *)d_ptr;
730 float* r_ptr_v4sf = (float *)r_ptr;
731
732
733 for(unsigned int i = 0 ; i < total_v4sf; i++, d_ptr_v4sf +=4, r_ptr_v4sf+=4 ) {
734 _mm_store_ps( d_ptr_v4sf, _mm_sub_ps( _mm_load_ps( d_ptr_v4sf),
735 _mm_load_ps( r_ptr_v4sf) ) );
736 }
737
738
739 r_ptr = (REAL32 *)r_ptr_v4sf;
740 d_ptr = (REAL32 *)d_ptr_v4sf;
741 for(unsigned int i=0; i < remainder; i++, r_ptr++, d_ptr++) {
742 *d_ptr -= *r_ptr;
743 }
744 }
745 else {
746 // Unordered case
747 const int* tab = s.siteTable().slice();
748
749 // Loop through the sites
750 for(int j=0; j < s.numSiteTable(); j++) {
751 int i = tab[j];
752
753 // Do the copy in the dumb way -- this could become quite complex
754 // Depending on whether the individual matrices are aligned or not.
755 d.elem(i).elem() -= l.elem(i).elem();
756
757 }
758
759 }
760}
761
762// Specialization to optimize the case
763// LatticeColorMatrix[Subset] -= LatticeColorMatrix * LatticeColorMatrix
764template<>
765void evaluate(OLattice< TCol >& d,
766 const OpSubtractAssign& op,
767 const QDPExpr<BinaryNode<OpMultiply,
768 Reference<QDPType< TCol, OLattice< TCol > > >,
769 Reference<QDPType< TCol, OLattice< TCol > > > >,
770 OLattice< TCol > >& rhs,
771 const Subset& s)
772{
773// cout << "call single site QDP_M_meq_M_times_M" << endl;
774
775 typedef OLattice< TCol > C;
776
777 const C& l = static_cast<const C&>(rhs.expression().left());
778 const C& r = static_cast<const C&>(rhs.expression().right());
779
780 PColorMatrix<RComplexFloat,3> tmp;
781
782
783 if( s.hasOrderedRep() ) {
784 for(int i=s.start(); i <= s.end(); i++) {
785 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
786 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
787 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
788
789 intrin_sse_mult_su3_nn(lm, rm, tmpm);
790
791 d.elem(i).elem().elem(0,0).real() -= tmp.elem(0,0).real();
792 d.elem(i).elem().elem(0,0).imag() -= tmp.elem(0,0).imag();
793 d.elem(i).elem().elem(0,1).real() -= tmp.elem(0,1).real();
794 d.elem(i).elem().elem(0,1).imag() -= tmp.elem(0,1).imag();
795 d.elem(i).elem().elem(0,2).real() -= tmp.elem(0,2).real();
796 d.elem(i).elem().elem(0,2).imag() -= tmp.elem(0,2).imag();
797
798 d.elem(i).elem().elem(1,0).real() -= tmp.elem(1,0).real();
799 d.elem(i).elem().elem(1,0).imag() -= tmp.elem(1,0).imag();
800 d.elem(i).elem().elem(1,1).real() -= tmp.elem(1,1).real();
801 d.elem(i).elem().elem(1,1).imag() -= tmp.elem(1,1).imag();
802 d.elem(i).elem().elem(1,2).real() -= tmp.elem(1,2).real();
803 d.elem(i).elem().elem(1,2).imag() -= tmp.elem(1,2).imag();
804
805 d.elem(i).elem().elem(2,0).real() -= tmp.elem(2,0).real();
806 d.elem(i).elem().elem(2,0).imag() -= tmp.elem(2,0).imag();
807 d.elem(i).elem().elem(2,1).real() -= tmp.elem(2,1).real();
808 d.elem(i).elem().elem(2,1).imag() -= tmp.elem(2,1).imag();
809 d.elem(i).elem().elem(2,2).real() -= tmp.elem(2,2).real();
810 d.elem(i).elem().elem(2,2).imag() -= tmp.elem(2,2).imag();
811 }
812 }
813 else {
814
815 const int *tab = s.siteTable().slice();
816 for(int j=0; j < s.numSiteTable(); j++) {
817 int i=tab[j];
818 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
819 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
820 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
821
822 intrin_sse_mult_su3_nn(lm, rm, tmpm);
823
824 d.elem(i).elem().elem(0,0).real() -= tmp.elem(0,0).real();
825 d.elem(i).elem().elem(0,0).imag() -= tmp.elem(0,0).imag();
826 d.elem(i).elem().elem(0,1).real() -= tmp.elem(0,1).real();
827 d.elem(i).elem().elem(0,1).imag() -= tmp.elem(0,1).imag();
828 d.elem(i).elem().elem(0,2).real() -= tmp.elem(0,2).real();
829 d.elem(i).elem().elem(0,2).imag() -= tmp.elem(0,2).imag();
830
831 d.elem(i).elem().elem(1,0).real() -= tmp.elem(1,0).real();
832 d.elem(i).elem().elem(1,0).imag() -= tmp.elem(1,0).imag();
833 d.elem(i).elem().elem(1,1).real() -= tmp.elem(1,1).real();
834 d.elem(i).elem().elem(1,1).imag() -= tmp.elem(1,1).imag();
835 d.elem(i).elem().elem(1,2).real() -= tmp.elem(1,2).real();
836 d.elem(i).elem().elem(1,2).imag() -= tmp.elem(1,2).imag();
837
838 d.elem(i).elem().elem(2,0).real() -= tmp.elem(2,0).real();
839 d.elem(i).elem().elem(2,0).imag() -= tmp.elem(2,0).imag();
840 d.elem(i).elem().elem(2,1).real() -= tmp.elem(2,1).real();
841 d.elem(i).elem().elem(2,1).imag() -= tmp.elem(2,1).imag();
842 d.elem(i).elem().elem(2,2).real() -= tmp.elem(2,2).real();
843 d.elem(i).elem().elem(2,2).imag() -= tmp.elem(2,2).imag();
844 }
845 }
846}
847
848
849// Specialization to optimize the case
850// LatticeColorMatrix[Subset] -= adj(LatticeColorMatrix) * LatticeColorMatrix
851template<>
852void evaluate(OLattice< TCol >& d,
853 const OpSubtractAssign& op,
854 const QDPExpr<BinaryNode<OpAdjMultiply,
855 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > >,
856 Reference<QDPType< TCol, OLattice< TCol > > > >,
857 OLattice< TCol > >& rhs,
858 const Subset& s)
859{
860// cout << "call single site QDP_M_meq_aM_times_M" << endl;
861
862 typedef OLattice< TCol > C;
863
864 const C& l = static_cast<const C&>(rhs.expression().left().child());
865 const C& r = static_cast<const C&>(rhs.expression().right());
866
867 PColorMatrix<RComplexFloat,3> tmp;
868 if( s.hasOrderedRep() ) {
869 for(int i=s.start(); i <= s.end(); i++) {
870 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
871 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
872 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
873
874 intrin_sse_mult_su3_an(lm, rm, tmpm);
875
876
877 d.elem(i).elem().elem(0,0).real() -= tmp.elem(0,0).real();
878 d.elem(i).elem().elem(0,0).imag() -= tmp.elem(0,0).imag();
879 d.elem(i).elem().elem(0,1).real() -= tmp.elem(0,1).real();
880 d.elem(i).elem().elem(0,1).imag() -= tmp.elem(0,1).imag();
881 d.elem(i).elem().elem(0,2).real() -= tmp.elem(0,2).real();
882 d.elem(i).elem().elem(0,2).imag() -= tmp.elem(0,2).imag();
883
884 d.elem(i).elem().elem(1,0).real() -= tmp.elem(1,0).real();
885 d.elem(i).elem().elem(1,0).imag() -= tmp.elem(1,0).imag();
886 d.elem(i).elem().elem(1,1).real() -= tmp.elem(1,1).real();
887 d.elem(i).elem().elem(1,1).imag() -= tmp.elem(1,1).imag();
888 d.elem(i).elem().elem(1,2).real() -= tmp.elem(1,2).real();
889 d.elem(i).elem().elem(1,2).imag() -= tmp.elem(1,2).imag();
890
891 d.elem(i).elem().elem(2,0).real() -= tmp.elem(2,0).real();
892 d.elem(i).elem().elem(2,0).imag() -= tmp.elem(2,0).imag();
893 d.elem(i).elem().elem(2,1).real() -= tmp.elem(2,1).real();
894 d.elem(i).elem().elem(2,1).imag() -= tmp.elem(2,1).imag();
895 d.elem(i).elem().elem(2,2).real() -= tmp.elem(2,2).real();
896 d.elem(i).elem().elem(2,2).imag() -= tmp.elem(2,2).imag();
897 }
898 }
899 else {
900
901 const int *tab = s.siteTable().slice();
902 for(int j=0; j < s.numSiteTable(); j++) {
903 int i=tab[j];
904 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
905 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
906 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
907
908 intrin_sse_mult_su3_an(lm, rm, tmpm);
909
910 d.elem(i).elem().elem(0,0).real() -= tmp.elem(0,0).real();
911 d.elem(i).elem().elem(0,0).imag() -= tmp.elem(0,0).imag();
912 d.elem(i).elem().elem(0,1).real() -= tmp.elem(0,1).real();
913 d.elem(i).elem().elem(0,1).imag() -= tmp.elem(0,1).imag();
914 d.elem(i).elem().elem(0,2).real() -= tmp.elem(0,2).real();
915 d.elem(i).elem().elem(0,2).imag() -= tmp.elem(0,2).imag();
916
917 d.elem(i).elem().elem(1,0).real() -= tmp.elem(1,0).real();
918 d.elem(i).elem().elem(1,0).imag() -= tmp.elem(1,0).imag();
919 d.elem(i).elem().elem(1,1).real() -= tmp.elem(1,1).real();
920 d.elem(i).elem().elem(1,1).imag() -= tmp.elem(1,1).imag();
921 d.elem(i).elem().elem(1,2).real() -= tmp.elem(1,2).real();
922 d.elem(i).elem().elem(1,2).imag() -= tmp.elem(1,2).imag();
923
924 d.elem(i).elem().elem(2,0).real() -= tmp.elem(2,0).real();
925 d.elem(i).elem().elem(2,0).imag() -= tmp.elem(2,0).imag();
926 d.elem(i).elem().elem(2,1).real() -= tmp.elem(2,1).real();
927 d.elem(i).elem().elem(2,1).imag() -= tmp.elem(2,1).imag();
928 d.elem(i).elem().elem(2,2).real() -= tmp.elem(2,2).real();
929 d.elem(i).elem().elem(2,2).imag() -= tmp.elem(2,2).imag();
930 }
931 }
932}
933
934
935// Specialization to optimize the case
936// LatticeColorMatrix[Subset] -= LatticeColorMatrix * adj(LatticeColorMatrix)
937template<>
938void evaluate(OLattice< TCol >& d,
939 const OpSubtractAssign& op,
940 const QDPExpr<BinaryNode<OpMultiplyAdj,
941 Reference<QDPType< TCol, OLattice< TCol > > >,
942 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > > >,
943 OLattice< TCol > >& rhs,
944 const Subset& s)
945{
946// cout << "call single site QDP_M_meq_M_times_aM" << endl;
947
948 typedef OLattice< TCol > C;
949
950 const C& l = static_cast<const C&>(rhs.expression().left());
951 const C& r = static_cast<const C&>(rhs.expression().right().child());
952
953 PColorMatrix<RComplexFloat,3> tmp;
954
955 if( s.hasOrderedRep() ) {
956 for(int i=s.start(); i <= s.end(); i++) {
957 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
958 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
959 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
960
961 intrin_sse_mult_su3_na(lm, rm, tmpm);
962
963 d.elem(i).elem().elem(0,0).real() -= tmp.elem(0,0).real();
964 d.elem(i).elem().elem(0,0).imag() -= tmp.elem(0,0).imag();
965 d.elem(i).elem().elem(0,1).real() -= tmp.elem(0,1).real();
966 d.elem(i).elem().elem(0,1).imag() -= tmp.elem(0,1).imag();
967 d.elem(i).elem().elem(0,2).real() -= tmp.elem(0,2).real();
968 d.elem(i).elem().elem(0,2).imag() -= tmp.elem(0,2).imag();
969
970 d.elem(i).elem().elem(1,0).real() -= tmp.elem(1,0).real();
971 d.elem(i).elem().elem(1,0).imag() -= tmp.elem(1,0).imag();
972 d.elem(i).elem().elem(1,1).real() -= tmp.elem(1,1).real();
973 d.elem(i).elem().elem(1,1).imag() -= tmp.elem(1,1).imag();
974 d.elem(i).elem().elem(1,2).real() -= tmp.elem(1,2).real();
975 d.elem(i).elem().elem(1,2).imag() -= tmp.elem(1,2).imag();
976
977 d.elem(i).elem().elem(2,0).real() -= tmp.elem(2,0).real();
978 d.elem(i).elem().elem(2,0).imag() -= tmp.elem(2,0).imag();
979 d.elem(i).elem().elem(2,1).real() -= tmp.elem(2,1).real();
980 d.elem(i).elem().elem(2,1).imag() -= tmp.elem(2,1).imag();
981 d.elem(i).elem().elem(2,2).real() -= tmp.elem(2,2).real();
982 d.elem(i).elem().elem(2,2).imag() -= tmp.elem(2,2).imag();
983 }
984 }
985 else {
986 const int *tab = s.siteTable().slice();
987 for(int j=0; j < s.numSiteTable(); j++) {
988 int i=tab[j];
989 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
990 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
991 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
992
993 intrin_sse_mult_su3_na(lm, rm, tmpm);
994
995 d.elem(i).elem().elem(0,0).real() -= tmp.elem(0,0).real();
996 d.elem(i).elem().elem(0,0).imag() -= tmp.elem(0,0).imag();
997 d.elem(i).elem().elem(0,1).real() -= tmp.elem(0,1).real();
998 d.elem(i).elem().elem(0,1).imag() -= tmp.elem(0,1).imag();
999 d.elem(i).elem().elem(0,2).real() -= tmp.elem(0,2).real();
1000 d.elem(i).elem().elem(0,2).imag() -= tmp.elem(0,2).imag();
1001
1002 d.elem(i).elem().elem(1,0).real() -= tmp.elem(1,0).real();
1003 d.elem(i).elem().elem(1,0).imag() -= tmp.elem(1,0).imag();
1004 d.elem(i).elem().elem(1,1).real() -= tmp.elem(1,1).real();
1005 d.elem(i).elem().elem(1,1).imag() -= tmp.elem(1,1).imag();
1006 d.elem(i).elem().elem(1,2).real() -= tmp.elem(1,2).real();
1007 d.elem(i).elem().elem(1,2).imag() -= tmp.elem(1,2).imag();
1008
1009 d.elem(i).elem().elem(2,0).real() -= tmp.elem(2,0).real();
1010 d.elem(i).elem().elem(2,0).imag() -= tmp.elem(2,0).imag();
1011 d.elem(i).elem().elem(2,1).real() -= tmp.elem(2,1).real();
1012 d.elem(i).elem().elem(2,1).imag() -= tmp.elem(2,1).imag();
1013 d.elem(i).elem().elem(2,2).real() -= tmp.elem(2,2).real();
1014 d.elem(i).elem().elem(2,2).imag() -= tmp.elem(2,2).imag();
1015 }
1016 }
1017}
1018
1019
1020// Specialization to optimize the case
1021// LatticeColorMatrix[Subset] -= adj(LatticeColorMatrix) * adj(LatticeColorMatrix)
1022template<>
1023void evaluate(OLattice< TCol >& d,
1024 const OpSubtractAssign& op,
1025 const QDPExpr<BinaryNode<OpAdjMultiplyAdj,
1026 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > >,
1027 UnaryNode<OpIdentity, Reference<QDPType< TCol, OLattice< TCol > > > > >,
1028 OLattice< TCol > >& rhs,
1029 const Subset& s)
1030{
1031// cout << "call single site QDP_M_meq_Ma_times_Ma" << endl;
1032
1033 typedef OLattice< TCol > C;
1034
1035 const C& l = static_cast<const C&>(rhs.expression().left().child());
1036 const C& r = static_cast<const C&>(rhs.expression().right().child());
1037
1038 PColorMatrix<RComplexFloat,3> tmp;
1039 if( s.hasOrderedRep() ) {
1040 for(int i=s.start(); i <= s.end(); i++) {
1041 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
1042 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
1043 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
1044
1045 intrin_sse_mult_su3_nn(rm, lm, tmpm);
1046
1047 // Take the adj(r*l) = adj(l)*adj(r)
1048 d.elem(i).elem().elem(0,0).real() -= tmp.elem(0,0).real();
1049 d.elem(i).elem().elem(0,0).imag() += tmp.elem(0,0).imag();
1050 d.elem(i).elem().elem(0,1).real() -= tmp.elem(1,0).real();
1051 d.elem(i).elem().elem(0,1).imag() += tmp.elem(1,0).imag();
1052 d.elem(i).elem().elem(0,2).real() -= tmp.elem(2,0).real();
1053 d.elem(i).elem().elem(0,2).imag() += tmp.elem(2,0).imag();
1054
1055 d.elem(i).elem().elem(1,0).real() -= tmp.elem(0,1).real();
1056 d.elem(i).elem().elem(1,0).imag() += tmp.elem(0,1).imag();
1057 d.elem(i).elem().elem(1,1).real() -= tmp.elem(1,1).real();
1058 d.elem(i).elem().elem(1,1).imag() += tmp.elem(1,1).imag();
1059 d.elem(i).elem().elem(1,2).real() -= tmp.elem(2,1).real();
1060 d.elem(i).elem().elem(1,2).imag() += tmp.elem(2,1).imag();
1061
1062 d.elem(i).elem().elem(2,0).real() -= tmp.elem(0,2).real();
1063 d.elem(i).elem().elem(2,0).imag() += tmp.elem(0,2).imag();
1064 d.elem(i).elem().elem(2,1).real() -= tmp.elem(1,2).real();
1065 d.elem(i).elem().elem(2,1).imag() += tmp.elem(1,2).imag();
1066 d.elem(i).elem().elem(2,2).real() -= tmp.elem(2,2).real();
1067 d.elem(i).elem().elem(2,2).imag() += tmp.elem(2,2).imag();
1068
1069 }
1070 }
1071 else {
1072 const int *tab = s.siteTable().slice();
1073 for(int j=0; j < s.numSiteTable(); j++) {
1074 int i=tab[j];
1075 su3_matrixf *lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
1076 su3_matrixf *rm = (su3_matrixf *)&(r.elem(i).elem().elem(0,0).real());
1077 su3_matrixf *tmpm = (su3_matrixf *)&(tmp.elem(0,0).real());
1078
1079 intrin_sse_mult_su3_nn(rm, lm, tmpm);
1080
1081
1082 // Take the adj(r*l) = adj(l)*adj(r)
1083 d.elem(i).elem().elem(0,0).real() -= tmp.elem(0,0).real();
1084 d.elem(i).elem().elem(0,0).imag() += tmp.elem(0,0).imag();
1085 d.elem(i).elem().elem(0,1).real() -= tmp.elem(1,0).real();
1086 d.elem(i).elem().elem(0,1).imag() += tmp.elem(1,0).imag();
1087 d.elem(i).elem().elem(0,2).real() -= tmp.elem(2,0).real();
1088 d.elem(i).elem().elem(0,2).imag() += tmp.elem(2,0).imag();
1089
1090 d.elem(i).elem().elem(1,0).real() -= tmp.elem(0,1).real();
1091 d.elem(i).elem().elem(1,0).imag() += tmp.elem(0,1).imag();
1092 d.elem(i).elem().elem(1,1).real() -= tmp.elem(1,1).real();
1093 d.elem(i).elem().elem(1,1).imag() += tmp.elem(1,1).imag();
1094 d.elem(i).elem().elem(1,2).real() -= tmp.elem(2,1).real();
1095 d.elem(i).elem().elem(1,2).imag() += tmp.elem(2,1).imag();
1096
1097 d.elem(i).elem().elem(2,0).real() -= tmp.elem(0,2).real();
1098 d.elem(i).elem().elem(2,0).imag() += tmp.elem(0,2).imag();
1099 d.elem(i).elem().elem(2,1).real() -= tmp.elem(1,2).real();
1100 d.elem(i).elem().elem(2,1).imag() += tmp.elem(1,2).imag();
1101 d.elem(i).elem().elem(2,2).real() -= tmp.elem(2,2).real();
1102 d.elem(i).elem().elem(2,2).imag() += tmp.elem(2,2).imag();
1103 }
1104 }
1105}
1106
1107
1108//-------------------------------------------------------------------
1109
1110// Specialization to optimize the case
1111// LatticeHalfFermion = LatticeColorMatrix * LatticeHalfFermion
1112// NOTE: let this be a subroutine to save space
1113template<>
1114void evaluate(OLattice< TVec2 >& d,
1115 const OpAssign& op,
1116 const QDPExpr<BinaryNode<OpMultiply,
1117 Reference<QDPType< TCol, OLattice< TCol > > >,
1118 Reference<QDPType< TVec2, OLattice< TVec2 > > > >,
1119 OLattice< TVec2 > >& rhs,
1120 const Subset& s)
1121{
1122#if defined(QDP_SCALARSITE_DEBUG)
1123 cout << "specialized QDP_H_M_times_H" << endl;
1124#endif
1125
1126 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
1127 typedef OLattice<PSpinVector<PColorVector<RComplexFloat, 3>, 2> > H;
1128
1129 const C& l = static_cast<const C&>(rhs.expression().left());
1130 const H& r = static_cast<const H&>(rhs.expression().right());
1131
1132
1133
1134 if( s.hasOrderedRep() ) {
1135 for(int i=s.start(); i <= s.end(); i++) {
1136
1137 su3_matrixf* lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
1138 half_wilson_vectorf* rh = (half_wilson_vectorf*)&(r.elem(i).elem(0).elem(0).real());
1139 half_wilson_vectorf* dh = (half_wilson_vectorf*)&(d.elem(i).elem(0).elem(0).real());
1140
1141 intrin_sse_mult_su3_mat_hwvec(lm,rh,dh);
1142
1143 }
1144 }
1145 else {
1146
1147 const int *tab = s.siteTable().slice();
1148 for(int j=0; j < s.numSiteTable(); j++) {
1149 int i=tab[j];
1150 su3_matrixf* lm = (su3_matrixf *)&(l.elem(i).elem().elem(0,0).real());
1151 half_wilson_vectorf* rh = (half_wilson_vectorf*)&(r.elem(i).elem(0).elem(0).real());
1152 half_wilson_vectorf* dh = (half_wilson_vectorf*)&(d.elem(i).elem(0).elem(0).real());
1153
1154 intrin_sse_mult_su3_mat_hwvec(lm,rh,dh);
1155 }
1156 }
1157}
1158#endif
1159
1160
1161//-------------------------------------------------------------------
1162// GNUC vector type
1163
1164//#define DEBUG_BLAS
1165
1166// AXPY and AXMY routines
1167void vaxpy3(REAL32 *Out,REAL32 *scalep,REAL32 *InScale, REAL32 *Add,int n_4vec)
1168{
1169#ifdef DEBUG_BLAS
1170 QDPIO::cout << "SSE_TEST: vaxpy3" << endl;
1171#endif
1172
1173// int n_loops = n_4vec >> 2; // only works on multiple of length 4 vectors
1174 int n_loops = n_4vec; // only works on multiple of length 24 vectors
1175
1176 __m128 vscalep = _mm_load_ss(scalep);
1177 vscalep = _mm_shuffle_ps(vscalep, vscalep, 0);
1178
1179
1180 for (; n_loops-- > 0; )
1181 {
1182 _mm_store_ps(Out+ 0, _mm_add_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+ 0)), _mm_load_ps(Add+ 0)));
1183 _mm_store_ps(Out+ 4, _mm_add_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+ 4)), _mm_load_ps(Add+ 4)));
1184 _mm_store_ps(Out+ 8, _mm_add_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+ 8)), _mm_load_ps(Add+ 8)));
1185 _mm_store_ps(Out+12, _mm_add_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+12)), _mm_load_ps(Add+12)));
1186 _mm_store_ps(Out+16, _mm_add_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+16)), _mm_load_ps(Add+16)));
1187 _mm_store_ps(Out+20, _mm_add_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+20)), _mm_load_ps(Add+20)));
1188
1189 Out += 24; InScale += 24; Add += 24;
1190 }
1191}
1192
1193
1194void vaxmy3(REAL32 *Out,REAL32 *scalep,REAL32 *InScale, REAL32 *Sub,int n_4vec)
1195{
1196#ifdef DEBUG_BLAS
1197 QDPIO::cout << "SSE_TEST: vaxmy3" << endl;
1198#endif
1199
1200// int n_loops = n_4vec >> 2; // only works on multiple of length 4 vectors
1201 int n_loops = n_4vec; // only works on multiple of length 24 vectors
1202
1203// v4sf va = load_v4sf((float *)&a);
1204 __m128 vscalep = _mm_load_ss(scalep);
1205 vscalep = _mm_shuffle_ps( vscalep, vscalep, 0);
1206
1207 for (; n_loops-- > 0; )
1208 {
1209 _mm_store_ps(Out+ 0, _mm_sub_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+ 0)), _mm_load_ps(Sub+ 0)));
1210 _mm_store_ps(Out+ 4, _mm_sub_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+ 4)), _mm_load_ps(Sub+ 4)));
1211 _mm_store_ps(Out+ 8, _mm_sub_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+ 8)), _mm_load_ps(Sub+ 8)));
1212 _mm_store_ps(Out+12, _mm_sub_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+12)), _mm_load_ps(Sub+12)));
1213 _mm_store_ps(Out+16, _mm_sub_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+16)), _mm_load_ps(Sub+16)));
1214 _mm_store_ps(Out+20, _mm_sub_ps(_mm_mul_ps(vscalep, _mm_load_ps(InScale+20)), _mm_load_ps(Sub+20)));
1215
1216 Out += 24; InScale += 24; Sub += 24;
1217 }
1218}
1219
1220
1221void vadd(REAL32 *Out, REAL32 *In1, REAL32 *In2, int n_4vec)
1222{
1223#ifdef DEBUG_BLAS
1224 QDPIO::cout << "SSE_TEST: vadd" << endl;
1225#endif
1226
1227// int n_loops = n_4vec >> 2; // only works on multiple of length 4 vectors
1228 int n_loops = n_4vec; // only works on multiple of length 24 vectors
1229
1230 for (; n_loops-- > 0; )
1231 {
1232 _mm_store_ps(Out+ 0, _mm_add_ps(_mm_load_ps(In1+ 0), _mm_load_ps(In2+ 0)));
1233 _mm_store_ps(Out+ 4, _mm_add_ps(_mm_load_ps(In1+ 4), _mm_load_ps(In2+ 4)));
1234 _mm_store_ps(Out+ 8, _mm_add_ps(_mm_load_ps(In1+ 8), _mm_load_ps(In2+ 8)));
1235 _mm_store_ps(Out+12, _mm_add_ps(_mm_load_ps(In1+12), _mm_load_ps(In2+12)));
1236 _mm_store_ps(Out+16, _mm_add_ps(_mm_load_ps(In1+16), _mm_load_ps(In2+16)));
1237 _mm_store_ps(Out+20, _mm_add_ps(_mm_load_ps(In1+20), _mm_load_ps(In2+20)));
1238
1239 Out += 24; In1 += 24; In2 += 24;
1240 }
1241}
1242
1243
1244void vsub(REAL32 *Out, REAL32 *In1, REAL32 *In2, int n_4vec)
1245{
1246#ifdef DEBUG_BLAS
1247 QDPIO::cout << "SSE_TEST: vsub" << endl;
1248#endif
1249
1250// int n_loops = n_4vec >> 2; // only works on multiple of length 4 vectors
1251 int n_loops = n_4vec; // only works on multiple of length 24 vectors
1252
1253 for (; n_loops-- > 0; )
1254 {
1255 _mm_store_ps(Out+ 0, _mm_sub_ps(_mm_load_ps(In1+ 0), _mm_load_ps(In2+ 0)));
1256 _mm_store_ps(Out+ 4, _mm_sub_ps(_mm_load_ps(In1+ 4), _mm_load_ps(In2+ 4)));
1257 _mm_store_ps(Out+ 8, _mm_sub_ps(_mm_load_ps(In1+ 8), _mm_load_ps(In2+ 8)));
1258 _mm_store_ps(Out+12, _mm_sub_ps(_mm_load_ps(In1+12), _mm_load_ps(In2+12)));
1259 _mm_store_ps(Out+16, _mm_sub_ps(_mm_load_ps(In1+16), _mm_load_ps(In2+16)));
1260 _mm_store_ps(Out+20, _mm_sub_ps(_mm_load_ps(In1+20), _mm_load_ps(In2+20)));
1261
1262 Out += 24; In1 += 24; In2 += 24;
1263 }
1264}
1265
1266void vscal(REAL32 *Out, REAL32 *scalep, REAL32 *In, int n_4vec)
1267{
1268#ifdef DEBUG_BLAS
1269 QDPIO::cout << "SSE_TEST: vadd" << endl;
1270#endif
1271
1272// int n_loops = n_4vec >> 2; // only works on multiple of length 4 vectors
1273 int n_loops = n_4vec; // only works on multiple of length 24 vectors
1274
1275// v4sf va = load_v4sf((float *)&a);
1276 __m128 vscalep = _mm_load_ss(scalep);
1277 vscalep = _mm_shuffle_ps(vscalep, vscalep, 0);
1278
1279 for (; n_loops-- > 0; )
1280 {
1281 _mm_store_ps(Out+ 0, _mm_mul_ps(vscalep, _mm_load_ps(In+ 0)));
1282 _mm_store_ps(Out+ 4, _mm_mul_ps(vscalep, _mm_load_ps(In+ 4)));
1283 _mm_store_ps(Out+ 8, _mm_mul_ps(vscalep, _mm_load_ps(In+ 8)));
1284 _mm_store_ps(Out+12, _mm_mul_ps(vscalep, _mm_load_ps(In+12)));
1285 _mm_store_ps(Out+16, _mm_mul_ps(vscalep, _mm_load_ps(In+16)));
1286 _mm_store_ps(Out+20, _mm_mul_ps(vscalep, _mm_load_ps(In+20)));
1287
1288 Out += 24; In += 24;
1289 }
1290}
1291
1292
1293void vaxpby3(REAL32 *Out, REAL32 *a, REAL32 *x, REAL32 *b, REAL32 *y, int n_4vec)
1294{
1295#ifdef DEBUG_BLAS
1296 QDPIO::cout << "SSE_TEST: vaxpby3: a*x+b*y" << endl;
1297#endif
1298
1299// int n_loops = n_4vec >> 2; // only works on multiple of length 4 vectors
1300 int n_loops = n_4vec; // only works on multiple of length 24 vectors
1301
1302// __m128 va = load_v4sf((float *)&a);
1303 __m128 va = _mm_load_ss(a);
1304 __m128 vb = _mm_load_ss(b);
1305 va = _mm_shuffle_ps(va, va, 0);
1306 vb = _mm_shuffle_ps(vb, vb, 0);
1307
1308
1309 for (; n_loops-- > 0; )
1310 {
1311 _mm_store_ps(Out+ 0, _mm_add_ps(_mm_mul_ps(va, _mm_load_ps(x+ 0)), _mm_mul_ps(vb, _mm_load_ps(y+ 0))));
1312 _mm_store_ps(Out+ 4, _mm_add_ps(_mm_mul_ps(va, _mm_load_ps(x+ 4)), _mm_mul_ps(vb, _mm_load_ps(y+ 4))));
1313 _mm_store_ps(Out+ 8, _mm_add_ps(_mm_mul_ps(va, _mm_load_ps(x+ 8)), _mm_mul_ps(vb, _mm_load_ps(y+ 8))));
1314 _mm_store_ps(Out+12, _mm_add_ps(_mm_mul_ps(va, _mm_load_ps(x+12)), _mm_mul_ps(vb, _mm_load_ps(y+12))));
1315 _mm_store_ps(Out+16, _mm_add_ps(_mm_mul_ps(va, _mm_load_ps(x+16)), _mm_mul_ps(vb, _mm_load_ps(y+16))));
1316 _mm_store_ps(Out+20, _mm_add_ps(_mm_mul_ps(va, _mm_load_ps(x+20)), _mm_mul_ps(vb, _mm_load_ps(y+20))));
1317
1318 Out += 24; x += 24; y += 24;
1319 }
1320}
1321
1322
1323void vaxmby3(REAL32 *Out, REAL32 *a, REAL32 *x, REAL32 *b, REAL32 *y, int n_4vec)
1324{
1325#ifdef DEBUG_BLAS
1326 QDPIO::cout << "SSE_TEST: vaxmby3: a*x-b*y" << endl;
1327#endif
1328
1329// int n_loops = n_4vec >> 2; // only works on multiple of length 4 vectors
1330 int n_loops = n_4vec; // only works on multiple of length 24 vectors
1331
1332// v4sf va = load_v4sf((float *)&a);
1333 __m128 va = _mm_load_ss(a);
1334 __m128 vb = _mm_load_ss(b);
1335 va = _mm_shuffle_ps( va, va, 0);
1336 vb = _mm_shuffle_ps( vb, vb, 0);
1337
1338 for (; n_loops-- > 0; )
1339 {
1340 _mm_store_ps(Out+ 0, _mm_sub_ps(_mm_mul_ps(va, _mm_load_ps(x+ 0)), _mm_mul_ps(vb, _mm_load_ps(y+ 0))));
1341 _mm_store_ps(Out+ 4, _mm_sub_ps(_mm_mul_ps(va, _mm_load_ps(x+ 4)), _mm_mul_ps(vb, _mm_load_ps(y+ 4))));
1342 _mm_store_ps(Out+ 8, _mm_sub_ps(_mm_mul_ps(va, _mm_load_ps(x+ 8)), _mm_mul_ps(vb, _mm_load_ps(y+ 8))));
1343 _mm_store_ps(Out+12, _mm_sub_ps(_mm_mul_ps(va, _mm_load_ps(x+12)), _mm_mul_ps(vb, _mm_load_ps(y+12))));
1344 _mm_store_ps(Out+16, _mm_sub_ps(_mm_mul_ps(va, _mm_load_ps(x+16)), _mm_mul_ps(vb, _mm_load_ps(y+16))));
1345 _mm_store_ps(Out+20, _mm_sub_ps(_mm_mul_ps(va, _mm_load_ps(x+20)), _mm_mul_ps(vb, _mm_load_ps(y+20))));
1346
1347 Out += 24; x += 24; y += 24;
1348 }
1349}
1350
1351
1352
1355
1356void local_sumsq_24_48(REAL64 *Out, REAL32 *In, int n_4vec)
1357{
1358
1359 __m128d vsum = _mm_setzero_pd();
1360 __m128d vsum2 = _mm_setzero_pd();
1361 __m128d lower_2, upper_2;
1362 __m128d dat_sq, dat_sq2;
1363 __m128 tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
1364
1365 REAL32* num=In;
1366 int loop_end;
1367
1368
1369
1370 loop_end = n_4vec-1;
1371 tmp1 = _mm_load_ps(num); num+=4;
1372
1373 for(int site=0; site < loop_end; site++) {
1374
1375 tmp3 = _mm_load_ps(num); num+=4; // Load 4
1376
1377 // First 4 numbers
1378 tmp2 = _mm_shuffle_ps(tmp1, tmp1, 0x0e); // flip numbers into tmp2
1379 lower_2 = _mm_cvtps_pd(tmp1); // convert to double
1380 upper_2 = _mm_cvtps_pd(tmp2);
1381
1382 // f*f
1383 dat_sq = _mm_mul_pd(lower_2, lower_2);
1384 vsum = _mm_add_pd(vsum, dat_sq);
1385 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1386 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1387 tmp4 = _mm_load_ps(num); num+=4; // Load 4
1388
1389 // 2nd 4 numbers
1390 tmp2 = _mm_shuffle_ps(tmp3, tmp3, 0x0e); // flip numbers into tmp2
1391 lower_2 = _mm_cvtps_pd(tmp3); // convert to double
1392 upper_2 = _mm_cvtps_pd(tmp2);
1393
1394 // f*f
1395 dat_sq = _mm_mul_pd(lower_2, lower_2);
1396 vsum = _mm_add_pd(vsum, dat_sq);
1397 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1398 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1399 tmp5 = _mm_load_ps(num); num+=4; // Load 4
1400
1401
1402 // 3rd 4 numbers
1403 tmp2 = _mm_shuffle_ps(tmp4, tmp4, 0x0e); // flip numbers into tmp2
1404 lower_2 = _mm_cvtps_pd(tmp4); // convert to double
1405 upper_2 = _mm_cvtps_pd(tmp2);
1406
1407 // f*f
1408 dat_sq = _mm_mul_pd(lower_2, lower_2);
1409 vsum = _mm_add_pd(vsum, dat_sq);
1410 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1411 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1412 tmp6 = _mm_load_ps(num); num+=4; // Load 4
1413
1414 // 4th 4 numbers
1415 tmp2 = _mm_shuffle_ps(tmp5, tmp5, 0x0e); // flip numbers into tmp2
1416 lower_2 = _mm_cvtps_pd(tmp5); // convert to double
1417 upper_2 = _mm_cvtps_pd(tmp2);
1418
1419 // f*f
1420 dat_sq = _mm_mul_pd(lower_2, lower_2);
1421 vsum = _mm_add_pd(vsum, dat_sq);
1422 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1423 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1424 tmp7 = _mm_load_ps(num); num+=4; // Load 4
1425
1426 // 5th 4 numbers
1427 tmp2 = _mm_shuffle_ps(tmp6, tmp6, 0x0e); // flip numbers into tmp2
1428 lower_2 = _mm_cvtps_pd(tmp6); // convert to double
1429 upper_2 = _mm_cvtps_pd(tmp2);
1430
1431 // f*f
1432 dat_sq = _mm_mul_pd(lower_2, lower_2);
1433 vsum = _mm_add_pd(vsum, dat_sq);
1434 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1435 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1436 tmp1 = _mm_load_ps(num); num+=4; // Load 4
1437
1438 // 6th 4 numbers
1439 tmp2 = _mm_shuffle_ps(tmp7, tmp7, 0x0e); // flip numbers into tmp2
1440 lower_2 = _mm_cvtps_pd(tmp7); // convert to double
1441 upper_2 = _mm_cvtps_pd(tmp2);
1442
1443 // f*f
1444 dat_sq = _mm_mul_pd(lower_2, lower_2);
1445 vsum = _mm_add_pd(vsum, dat_sq);
1446 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1447 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1448
1449
1450 }
1451 /* Last one */
1452 {
1453 tmp3 = _mm_load_ps(num); num+=4; // Load 4
1454
1455 // First 4 numbers
1456 tmp2 = _mm_shuffle_ps(tmp1, tmp1, 0x0e); // flip numbers into tmp2
1457 lower_2 = _mm_cvtps_pd(tmp1); // convert to double
1458 upper_2 = _mm_cvtps_pd(tmp2);
1459
1460 // f*f
1461 dat_sq = _mm_mul_pd(lower_2, lower_2);
1462 vsum = _mm_add_pd(vsum, dat_sq);
1463 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1464 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1465
1466 tmp4 = _mm_load_ps(num); num+=4; // Load 4
1467
1468 // 2nd 4 numbers
1469 tmp2 = _mm_shuffle_ps(tmp3, tmp3, 0x0e); // flip numbers into tmp2
1470 lower_2 = _mm_cvtps_pd(tmp3); // convert to double
1471 upper_2 = _mm_cvtps_pd(tmp2);
1472
1473 // f*f
1474 dat_sq = _mm_mul_pd(lower_2, lower_2);
1475 vsum = _mm_add_pd(vsum, dat_sq);
1476 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1477 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1478
1479 tmp5 = _mm_load_ps(num); num+=4; // Load 4
1480
1481 // 3rd 4 numbers
1482 tmp2 = _mm_shuffle_ps(tmp4, tmp4, 0x0e); // flip numbers into tmp2
1483 lower_2 = _mm_cvtps_pd(tmp4); // convert to double
1484 upper_2 = _mm_cvtps_pd(tmp2);
1485
1486 // f*f
1487 dat_sq = _mm_mul_pd(lower_2, lower_2);
1488 vsum = _mm_add_pd(vsum, dat_sq);
1489 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1490 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1491
1492 tmp6 = _mm_load_ps(num); num+=4; // Load 4
1493
1494 // 4th 4 numbers
1495 tmp2 = _mm_shuffle_ps(tmp5, tmp5, 0x0e); // flip numbers into tmp2
1496 lower_2 = _mm_cvtps_pd(tmp5); // convert to double
1497 upper_2 = _mm_cvtps_pd(tmp2);
1498
1499 // f*f
1500 dat_sq = _mm_mul_pd(lower_2, lower_2);
1501 vsum = _mm_add_pd(vsum, dat_sq);
1502 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1503 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1504
1505 tmp7 = _mm_load_ps(num); num+=4; // Load 4
1506
1507 // 5th 4 numbers
1508 tmp2 = _mm_shuffle_ps(tmp6, tmp6, 0x0e); // flip numbers into tmp2
1509 lower_2 = _mm_cvtps_pd(tmp6); // convert to double
1510 upper_2 = _mm_cvtps_pd(tmp2);
1511
1512 // f*f
1513 dat_sq = _mm_mul_pd(lower_2, lower_2);
1514 vsum = _mm_add_pd(vsum, dat_sq);
1515 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1516 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1517
1518 // 6th 4 numbers
1519 tmp2 = _mm_shuffle_ps(tmp7, tmp7, 0x0e); // flip numbers into tmp2
1520 lower_2 = _mm_cvtps_pd(tmp7); // convert to double
1521 upper_2 = _mm_cvtps_pd(tmp2);
1522
1523 // f*f
1524 dat_sq = _mm_mul_pd(lower_2, lower_2);
1525 vsum = _mm_add_pd(vsum, dat_sq);
1526 dat_sq2 = _mm_mul_pd(upper_2, upper_2);
1527 vsum2 = _mm_add_pd(vsum2, dat_sq2);
1528 }
1529
1530
1531 vsum = _mm_add_pd(vsum, vsum2);
1532 /* Now sum horizontally on a and b */
1533 /* Move high word of vsum to low of lower_2 */
1534 /* Cross operation: SRC = vsum = | sum_1 | sum_0 |
1535 DEST= vsum = | sum 1 | sum_0 |
1536
1537 Bit field 0 != 0 => DEST[63-0] = DEST[127-64] = sum_1
1538 Bit field 1 == 0 => DEST[127:64] ← SRC[63:0] = sum_0 ;
1539 So output register: sum_0 | sum_1
1540 */
1541 lower_2 = _mm_shuffle_pd(vsum, vsum, 0x01);
1542
1543 /* Now sum up: lower_2 = sum_0 | sum_1
1544 + vsum = sum_1 | sum_0
1545 ========================
1546 Sum | Sum
1547 */
1548
1549 upper_2 = _mm_add_pd(lower_2, vsum);
1550
1551 /* Sum a should now have 2 copies of the sum */
1552 /* I should store the result */
1553 _mm_storeh_pd((double *)Out, upper_2);
1554
1555 // *Out=result;
1556}
1557
1558
1559void local_vcdot(REAL64 *Out_re, REAL64 *Out_im, REAL32 *V1, REAL32 *V2, int n_3vec)
1560{
1561 double result_re;
1562 double result_im;
1563
1564
1565 double v1_0r;
1566 double v1_0i;
1567 double v1_1r;
1568 double v1_1i;
1569 double v1_2r;
1570 double v1_2i;
1571
1572 double v2_0r;
1573 double v2_0i;
1574 double v2_1r;
1575 double v2_1i;
1576 double v2_2r;
1577 double v2_2i;
1578
1579 int counter=0;
1580 unsigned long vecptr1=0;
1581 unsigned long vecptr2=0;
1582 result_re= 0;
1583 result_im= 0;
1584
1585
1586 if( n_3vec > 0 ) {
1587
1588 v1_0r = (REAL64)V1[vecptr1++];
1589 v2_0r = (REAL64)V2[vecptr2++];
1590
1591 v1_0i = (REAL64)V1[vecptr1++];
1592 v2_0i = (REAL64)V2[vecptr2++];
1593
1594 v1_1r = (REAL64)V1[vecptr1++];
1595 v2_1r = (REAL64)V2[vecptr2++];
1596
1597 for(counter=0; counter < n_3vec-1; counter++) {
1598
1599
1600
1601 result_re = result_re + v1_0r*v2_0r;
1602 v1_1i =(REAL64)V1[vecptr1++];
1603 result_im = result_im - v1_0i*v2_0r;
1604 v2_1i = (REAL64)V2[vecptr2++];
1605 result_im = result_im + v1_0r*v2_0i;
1606 v1_2r = (REAL64)V1[vecptr1++];
1607 result_re = result_re + v1_0i*v2_0i;
1608 v2_2r = (REAL64)V2[vecptr2++];
1609
1610 result_re = result_re + v1_1r*v2_1r;
1611 v1_2i = (REAL64)V1[vecptr1++];
1612 result_im = result_im - v1_1i*v2_1r;
1613 v2_2i = (REAL64)V2[vecptr2++];
1614 result_im = result_im + v1_1r*v2_1i;
1615 v1_0r = (REAL64)V1[vecptr1++];
1616 result_re = result_re + v1_1i*v2_1i;
1617 v2_0r = (REAL64)V2[vecptr2++];
1618
1619 result_re = result_re + v1_2r*v2_2r;
1620 v1_0i = (REAL64)V1[vecptr1++];
1621 result_im = result_im - v1_2i*v2_2r;
1622 v2_0i = (REAL64)V2[vecptr2++];
1623 result_im = result_im + v1_2r*v2_2i;
1624 v1_1r = (REAL64)V1[vecptr1++];
1625 result_re = result_re + v1_2i*v2_2i;
1626 v2_1r = (REAL64)V2[vecptr2++];
1627
1628 }
1629
1630 // Last one plus drain...
1631 result_re = result_re + v1_0r*v2_0r;
1632 v1_1i =(REAL64)V1[vecptr1++];
1633 result_im = result_im - v1_0i*v2_0r;
1634 v2_1i = (REAL64)V2[vecptr2++];
1635 result_im = result_im + v1_0r*v2_0i;
1636 v1_2r = (REAL64)V1[vecptr1++];
1637 result_re = result_re + v1_0i*v2_0i;
1638 v2_2r = (REAL64)V2[vecptr2++];
1639
1640 result_re = result_re + v1_1r*v2_1r;
1641 v1_2i = (REAL64)V1[vecptr1++];
1642 result_im = result_im - v1_1i*v2_1r;
1643 v2_2i = (REAL64)V2[vecptr2++];
1644
1645 result_im = result_im + v1_1r*v2_1i;
1646 result_re = result_re + v1_1i*v2_1i;
1647
1648
1649 result_re = result_re + v1_2r*v2_2r;
1650 result_im = result_im - v1_2i*v2_2r;
1651 result_im = result_im + v1_2r*v2_2i;
1652 result_re = result_re + v1_2i*v2_2i;
1653
1654 }
1655
1656 *Out_re=(REAL64)result_re;
1657 *Out_im=(REAL64)result_im;
1658}
1659
1660
1661void local_vcdot_real(REAL64 *Out, REAL32 *V1, REAL32 *V2, int n_3vec)
1662{
1663 REAL64 result;
1664
1665 REAL64 v1_0r;
1666 REAL64 v1_0i;
1667 REAL64 v1_1r;
1668 REAL64 v1_1i;
1669 REAL64 v1_2r;
1670 REAL64 v1_2i;
1671
1672 REAL64 v2_0r;
1673 REAL64 v2_0i;
1674 REAL64 v2_1r;
1675 REAL64 v2_1i;
1676 REAL64 v2_2r;
1677 REAL64 v2_2i;
1678
1679 int counter=0;
1680 unsigned long vecptr1=0;
1681 unsigned long vecptr2=0;
1682 result= 0;
1683
1684
1685 if( n_3vec > 0 ) {
1686
1687 // Prefetch
1688 v1_0r = (REAL64)V1[vecptr1++];
1689 v2_0r = (REAL64)V2[vecptr2++];
1690
1691 v1_0i = (REAL64)V1[vecptr1++];
1692 v2_0i = (REAL64)V2[vecptr2++];
1693
1694 v1_1r = (REAL64)V1[vecptr1++];
1695 v2_1r = (REAL64)V2[vecptr2++];
1696
1697 v1_1i =(REAL64)V1[vecptr1++];
1698 v2_1i = (REAL64)V2[vecptr2++];
1699
1700 v1_2r = (REAL64)V1[vecptr1++];
1701 v2_2r = (REAL64)V2[vecptr2++];
1702
1703 v1_2i = (REAL64)V1[vecptr1++];
1704 v2_2i = (REAL64)V2[vecptr2++];
1705
1706 for(counter=0; counter < n_3vec-1; counter++) {
1707 result = result + v1_0r*v2_0r;
1708 v1_0r = (REAL64)V1[vecptr1++];
1709 v2_0r = (REAL64)V2[vecptr2++];
1710
1711 result = result + v1_0i*v2_0i;
1712 v1_0i = (REAL64)V1[vecptr1++];
1713 v2_0i = (REAL64)V2[vecptr2++];
1714
1715 result = result + v1_1r*v2_1r;
1716 v1_1r = (REAL64)V1[vecptr1++];
1717 v2_1r = (REAL64)V2[vecptr2++];
1718
1719 result = result + v1_1i*v2_1i;
1720 v1_1i =(REAL64)V1[vecptr1++];
1721 v2_1i = (REAL64)V2[vecptr2++];
1722
1723 result = result + v1_2r*v2_2r;
1724 v1_2r = (REAL64)V1[vecptr1++];
1725 v2_2r = (REAL64)V2[vecptr2++];
1726
1727 result = result + v1_2i*v2_2i;
1728 v1_2i = (REAL64)V1[vecptr1++];
1729 v2_2i = (REAL64)V2[vecptr2++];
1730
1731
1732 }
1733
1734 // Last one plus drain...
1735 result = result + v1_0r*v2_0r;
1736 result = result + v1_0i*v2_0i;
1737 result = result + v1_1r*v2_1r;
1738 result = result + v1_1i*v2_1i;
1739 result = result + v1_2r*v2_2r;
1740 result = result + v1_2i*v2_2i;
1741
1742 }
1743
1744 *Out=(REAL64)result;
1745}
1746
1747
1748
1749
1750} // namespace QDP;
1751
1752#endif // defined(__GNUC__)
QDP::OLattice
Outer grid Lattice type.
Definition qdp_outer.h:264
QDP::PColorMatrix
Primitive color Matrix class.
Definition qdp_primcolormat.h:26
QDP::QDPExpr
Expression class for QDP.
Definition qdp_qdpexpr.h:16
QDP::QDPType
QDPType - major type class/container for all QDP objects.
Definition qdp_qdptype.h:29
QDP::Subset
Subsets - controls how lattices are looped.
Definition qdp_subset.h:39
QDP::REAL32
float REAL32
Definition qdp_scalarsite_defs.h:21
QDP::REAL64
double REAL64
Definition qdp_scalarsite_defs.h:22
QDP::H
OLattice< PSpinVector< PColorVector< RComplexFloat, 3 >, 2 > > H
Definition qdp_scalarsite_generic_linalg.h:13
QDP::evaluate
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)
Definition qdp_scalarsite_linalg_double.cc:78
QDP::C
OLattice< PScalar< PColorMatrix< RComplexFloat, 3 > > > C
Definition qdp_scalarsite_generic_linalg.h:12
QDP::QDPIO::cout
StandardOutputStream cout
Definition qdp_stdio.cc:21
QDP
Yet another random number generator.
Definition qdp_bgq_threadbind.cc:16
QDP::local_sumsq_24_48
void local_sumsq_24_48(REAL64 *Out, REAL32 *In, int n_3vec)
QDP::real
MakeReturn< UnaryNode< FnReal, typenameCreateLeaf< QDPExpr< T1, C1 > >::Leaf_t >, typenameUnaryReturn< C1, FnReal >::Type_t >::Expression_t real(const QDPExpr< T1, C1 > &l)
Definition qdp.h:4972
QDP::local_vcdot_real
void local_vcdot_real(REAL64 *Out_re, REAL32 *V1, REAL32 *V2, int n_3vec)
QDP::local_vcdot
void local_vcdot(REAL64 *Out_re, REAL64 *Out_im, REAL32 *V1, REAL32 *V2, int n_3vec)
QDP::vaxpby3
void vaxpby3(REAL *Out, REAL *ap, REAL *xp, REAL *bp, REAL *yp, int n_3vec)
Definition generic_blas_vaxpby3.h:15
QDP::vsub
void vsub(REAL *Out, REAL *In1, REAL *In2, int n_3vec)
Definition generic_blas_vsub.h:15
QDP::vaxpy3
void vaxpy3(REAL *Out, REAL *scalep, REAL *InScale, REAL *Add, int n_4vec)
Definition generic_blas_vaxpy3.h:14
QDP::vadd
void vadd(REAL *Out, REAL *In1, REAL *In2, int n_3vec)
Definition generic_blas_vadd.h:15
QDP::vaxmy3
void vaxmy3(REAL *Out, REAL *scalep, REAL *InScale, REAL *Sub, int n_3vec)
Definition generic_blas_vaxmy3.h:15
QDP::vaxmby3
void vaxmby3(REAL *Out, REAL *ap, REAL *xp, REAL *bp, REAL *yp, int n_3vec)
Definition generic_blas_vaxmby3.h:14
QDP::vscal
void vscal(REAL *Out, REAL *scalep, REAL *In, int n_3vec)
Definition generic_blas_vscal.h:15
qdp.h
Primary include file for QDP.
qdp_sse_intrin.h
Generated by doxygen 1.16.1