QDP++
qdp_scalarsite_sse_blas.h
Go to the documentation of this file.
1// $Id: qdp_scalarsite_sse_blas.h,v 1.30 2009-09-15 20:48:51 bjoo Exp $
7
8#ifndef QDP_SCALARSITE_SSE_BLAS_H
9#define QDP_SCALARSITE_SSE_BLAS_H
10
11#include "qdp_config.h"
12
13namespace QDP {
14
15 namespace ThreadReductions{
17 };
19 // Forward declarations of BLAS routines
20 void vaxpy3(REAL32 *Out, REAL32 *scalep,REAL32 *InScale, REAL32 *Add,int n_3vec);
21 void vaxmy3(REAL32 *Out, REAL32 *scalep,REAL32 *InScale, REAL32 *Sub,int n_3vec);
22 void vadd(REAL32 *Out, REAL32 *In1, REAL32 *In2, int n_3vec);
23 void vsub(REAL32 *Out, REAL32 *In1, REAL32 *In2, int n_3vec);
24 void vscal(REAL32 *Out, REAL32 *scalep, REAL32 *In, int n_3vec);
25 void vaxpby3(REAL32* Out, REAL32* a, REAL32* x, REAL32* b, REAL32* y, int n_3vec);
26 void vaxmby3(REAL32* Out, REAL32* a, REAL32* x, REAL32* b, REAL32* y, int n_3vec);
27
28 void local_sumsq_24_48(REAL64 *Out, REAL32 *In, int n_3vec);
29
30 void local_vcdot(REAL64 *Out_re, REAL64 *Out_im, REAL32 *V1, REAL32 *V2, int n_3vec);
31 void local_vcdot_real(REAL64 *Out_re, REAL32 *V1, REAL32 *V2, int n_3vec);
32
33
36
37
39 // Threading evaluates
40 //
41 // by Xu Guo, EPCC, 6 October, 2008
43 // the wrappers for the functions to be threaded
45
48 /* #define DEBUG_BLAS_VAXMBY */
49 /* #define DEBUG_BLAS_VAXPBY */
50
51#define QDP_SCALARSITE_USE_EVALUATE
52
53
54// TVec is the LatticeFermion from qdp_dwdefs.h with the OLattice<> stripped
55// from around it
57// TScalar is the usual Real, with the OScalar<> stripped from it
58//
59// THis is simply to make the code more readable, and reduces < < s and > >s
60// in the template arguments
61
62
63#if defined(QDP_SCALARSITE_USE_EVALUATE)
64
65// d += Scalar*Vec
66template<>
67inline
69 const OpAddAssign& op,
71 Reference< QDPType< TScal, OScalar < TScal > > >,
73 OLattice< TVec > > &rhs,
74 const Subset& s)
75{
76
77#ifdef DEBUG_BLAS
78 QDPIO::cout << "SSE: y += a*x" << endl;
79#endif
80
81 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().right());
82 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().left());
83
84 REAL32 ar = a.elem().elem().elem().elem();
85 REAL32* aptr = &ar;
86 if( s.hasOrderedRep() ) {
87 REAL32* xptr = (REAL32 *)&(x.elem(s.start()).elem(0).elem(0).real());
88 REAL32* yptr = &(d.elem(s.start()).elem(0).elem(0).real());
89 // cout << "Specialised axpy a ="<< ar << endl;
90
91 int total_n_3vec = (s.end()-s.start()+1);
92
93 ordered_sse_vaxOpy3_user_arg arg = {yptr, aptr, xptr, yptr, vaxpy3};
94
98 // Original code
100 //int n_3vec = (s.end()-s.start()+1)*24;
101 //vaxpy3(yptr, aptr, xptr, yptr, n_3vec);
102 }
103 else {
104
105 const int* tab = s.siteTable().slice();
107 int totalSize = s.numSiteTable();
112
113
114 // Original code
116 /*
117 for(int j=0; j < s.numSiteTable(); j++) {
118 int i = tab[j];
119
120 REAL32* xptr = (REAL32 *)&(x.elem(i).elem(0).elem(0).real());
121 REAL32* yptr = &(d.elem(i).elem(0).elem(0).real());
122
123 vaxpy3(yptr, aptr, xptr, yptr, 24);
124 }
125 */
126 }
127
128
129 }
130
131// d -= Scalar*Vec
132template<>
133inline
135 const OpSubtractAssign& op,
137 Reference< QDPType< TScal, OScalar < TScal > > >,
139 OLattice< TVec > > &rhs,
140 const Subset& s)
141 {
142
143#ifdef DEBUG_BLAS
144 QDPIO::cout << "SSE: y -= a*x" << endl;
145#endif
147 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().right());
148 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().left());
149
150 // - sign as y -= ax <=> y = y-ax = -ax + y = axpy with -a
151 REAL32 ar = -( a.elem().elem().elem().elem());
152 REAL32* aptr = &ar;
153 if( s.hasOrderedRep() ) {
154 REAL32* xptr = (REAL32 *)&(x.elem(s.start()).elem(0).elem(0).real());
155 REAL32* yptr = &(d.elem(s.start()).elem(0).elem(0).real());
156
157 int total_n_3vec = (s.end()-s.start()+1);
158
159 ordered_sse_vaxOpy3_user_arg arg = {yptr, aptr, xptr, yptr, vaxpy3};
160
162
164 // Original code
166 //int n_3vec = (s.end()-s.start()+1)*24;
167 //vaxpy3(yptr, aptr, xptr, yptr, n_3vec);
168 }
169 else {
170 const int* tab = s.siteTable().slice();
172 int totalSize = s.numSiteTable();
173
175
177
179 // Original code
181 /*
182 for(int j=0; j < s.numSiteTable(); j++) {
183 int i = tab[j];
184
185 REAL32* xptr = (REAL32 *)&(x.elem(i).elem(0).elem(0).real());
186 REAL32* yptr = &(d.elem(i).elem(0).elem(0).real());
187 vaxpy3(yptr, aptr, xptr, yptr, 24);
188
189 }*/
190 }
191
192 }
193
194// z = ax + y
195template<>
196inline
198 const OpAssign &op,
199 const QDPExpr<
206 const Subset& s)
207{
208
209#ifdef DEBUG_BLAS
210 QDPIO::cout << "SSE: z = a*x + y" << endl;
211#endif
213 // Peel the stuff out of the expression
214 // y is the right side of rhs
215 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right());
217 // ax is the left side of rhs and is in a binary node
218 typedef BinaryNode<OpMultiply,
222 // get the binary node
223 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
224
225 // get a and x out of the bynary node
226 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
227 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
228 // Set pointers
229 REAL32 ar = a.elem().elem().elem().elem();
230 REAL32 *aptr = (REAL32 *)&ar;
231 if( s.hasOrderedRep() ) {
232 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
233 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
234 REAL32* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
235
236 int total_n_3vec = (s.end()-s.start()+1);
237
238 ordered_sse_vaxOpy3_user_arg arg = {zptr, aptr, xptr, yptr, vaxpy3};
239
241
243 // Original code
245 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
246 //int n_3vec = (s.end()-s.start()+1)*24;
247 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
248 }
249 else {
250 const int* tab = s.siteTable().slice();
252 int totalSize = s.numSiteTable();
253
254 unordered_sse_vaxOpy3_z_user_arg arg(x, y, d, aptr, tab, vaxpy3);
255
259 // Original code
261 /*
262 for(int j=0; j < s.numSiteTable(); j++) {
263 int i = tab[j];
264 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
265 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
266 REAL32* zptr = &(d.elem(i).elem(0).elem(0).real());
267
268 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
269 vaxpy3(zptr, aptr, xptr, yptr, 24);
270
271 }*/
273
274 }
275
276// Vec = Vec + Scal*Vec
277template<>
278inline
280 const OpAssign &op,
281 const QDPExpr<
287 OLattice< TVec > > &rhs,
288 const Subset& s)
289{
290#ifdef DEBUG_BLAS
291 QDPIO::cout << "SSE: z = y + a*x" << endl;
292#endif
293
294
295 // Peel the stuff out of the expression
296
297 // y is the left side of rhs
298 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().left());
299
300 // ax is the right side of rhs and is in a binary node
301 typedef BinaryNode<OpMultiply,
304
305 // get the binary node
306 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
307
308 // get a and x out of the bynary node
309 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
310 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
311 // Set pointers
312 REAL32 ar = a.elem().elem().elem().elem();
313 REAL32 *aptr = (REAL32 *)&ar;
314 if( s.hasOrderedRep() ) {
315 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
316 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
317 REAL32* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
318
319 int total_n_3vec = (s.end()-s.start()+1);
320
321 ordered_sse_vaxOpy3_user_arg arg = {zptr, aptr, xptr, yptr, vaxpy3};
322
324
326 // Original code
328 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
329 //int n_3vec = (s.end()-s.start()+1)*24;
330 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
331 }
332 else {
333 const int* tab = s.siteTable().slice();
334
335 int totalSize = s.numSiteTable();
336
337 unordered_sse_vaxOpy3_z_user_arg arg(x, y, d, aptr, tab, vaxpy3);
338
342 // Original code
344 /*
345 for(int j=0; j < s.numSiteTable(); j++) {
346 int i = tab[j];
347 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
348 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
349 REAL32* zptr = &(d.elem(i).elem(0).elem(0).real());
350
351 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
352 vaxpy3(zptr, aptr, xptr, yptr, 24);
355 }*/
358
359// Vec = Scalar*Vec - Vec
360template<>
361inline
363 const OpAssign &op,
364 const QDPExpr<
370 OLattice< TVec > > &rhs,
371 const Subset& s)
372{
373#ifdef DEBUG_BLAS
374 QDPIO::cout << "SSE: z = a*x - y" << endl;
375#endif
376
377
378 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right());
379
380 // ax is the left side of rhs and is in a binary node
381 typedef BinaryNode<OpMultiply,
384
385 // get the binary node
386 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
388 // get a and x out of the bynary node
389 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
390 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
391 // Set pointers
392 REAL32 ar = a.elem().elem().elem().elem();
393 REAL32 *aptr = (REAL32 *)&ar;
394 if( s.hasOrderedRep() ) {
395
396 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
397 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
398 REAL32* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
399
400 int total_n_3vec = (s.end()-s.start()+1);
401
402 ordered_sse_vaxOpy3_user_arg arg = {zptr, aptr, xptr, yptr, vaxmy3};
403
405
407 // Original code
409 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
410 //int n_3vec = (s.end()-s.start()+1)*24;
411 //vaxmy3(zptr, aptr, xptr, yptr, n_3vec);
412 }
413 else {
414 const int* tab = s.siteTable().slice();
415
416 int totalSize = s.numSiteTable();
417
418 unordered_sse_vaxOpy3_z_user_arg arg(x, y, d, aptr, tab, vaxmy3);
419
421
423 // Original code
425 /*
426 for(int j=0; j < s.numSiteTable(); j++) {
427 int i = tab[j];
428 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
429 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
430 REAL32* zptr = &(d.elem(i).elem(0).elem(0).real());
431
432
433 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
434
435 vaxmy3(zptr, aptr, xptr, yptr,24);
436 }*/
437 }
438 }
439
440template<>
441inline
443 const OpAssign &op,
444 const QDPExpr<
450 OLattice< TVec > > &rhs,
451 const Subset& s)
452{
453#ifdef DEBUG_BLAS
454 QDPIO::cout << "SSE: z = y - a*x" << endl;
455#endif
456
457 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().left());
458
459 // ax is the right side of rhs and is in a binary node
460 typedef BinaryNode<OpMultiply,
463
464 // get the binary node
465 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
466
467 // get a and x out of the bynary node
468 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
469 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
470 // Set pointers etc.
471
472 // -ve sign as y - ax = -ax + y = axpy with -a.
473 REAL32 ar = -a.elem().elem().elem().elem();
474 REAL32 *aptr = (REAL32 *)&ar;
475 if( s.hasOrderedRep()) {
476 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
477 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
478 REAL32* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
479
480 int total_n_3vec = (s.end()-s.start()+1);
481
482 ordered_sse_vaxOpy3_user_arg arg = {zptr, aptr, xptr, yptr, vaxpy3};
483
485
487 // Original code
489 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
490 //int n_3vec = (s.end()-s.start()+1)*24;
491 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
492 }
493 else {
494 const int* tab = s.siteTable().slice();
495
496 int totalSize = s.numSiteTable();
497
498 unordered_sse_vaxOpy3_z_user_arg arg(x, y, d, aptr, tab, vaxpy3);
499
501
503 // Original code
505 /*
506 for(int j=0; j < s.numSiteTable(); j++) {
507 int i = tab[j];
508 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
509 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
510 REAL32* zptr = &(d.elem(i).elem(0).elem(0).real());
511 vaxpy3(zptr, aptr, xptr, yptr, 24);
512 }*/
513 }
514 }
515
516// Vec += Vec * Scalar (AXPY)
517template<>
518inline
520 const OpAddAssign& op,
523 Reference< QDPType< TScal, OScalar < TScal > > > >,
524 OLattice< TVec > > &rhs,
525 const Subset& s)
526{
527
528#ifdef DEBUG_BLAS
529 QDPIO::cout << "SSE: y += x*a" << endl;
530#endif
531
532 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().left());
533 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().right());
534
535 REAL32 ar = a.elem().elem().elem().elem();
536 REAL32* aptr = &ar;
537 if( s.hasOrderedRep() ) {
538
539 REAL32* xptr = (REAL32 *)&(x.elem(s.start()).elem(0).elem(0).real());
540 REAL32* yptr = &(d.elem(s.start()).elem(0).elem(0).real());
541
542 int total_n_3vec = (s.end()-s.start()+1);
543
544 ordered_sse_vaxOpy3_user_arg arg = {yptr, aptr, xptr, yptr, vaxpy3};
545
547
549 // Original code
551 //int n_3vec = (s.end()-s.start()+1)*24;
552 //vaxpy3(yptr, aptr, xptr, yptr, n_3vec);
553 }
554 else {
555 const int* tab = s.siteTable().slice();
556
557 int totalSize = s.numSiteTable();
558
559 unordered_sse_vaxOpy3_y_user_arg arg(x, d, aptr, tab, 1, vaxpy3);
560
562
564 // Original code
566 /*
567 for(int j=0; j < s.numSiteTable(); j++) {
568 int i = tab[j];
569 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
570 REAL32 *yptr = (REAL32 *) &(d.elem(i).elem(0).elem(0).real());
571 vaxpy3(yptr, aptr, xptr, yptr,24);
572 }*/
573 }
574
575 }
576
577
578// Vec -= Vec *Scalar
579template<>
580inline
582 const OpSubtractAssign& op,
585 Reference< QDPType< TScal, OScalar < TScal > > > >,
586 OLattice< TVec > > &rhs,
587 const Subset& s)
588{
589
590#ifdef DEBUG_BLAS
591 QDPIO::cout << "SSE: y -= x*a" << endl;
592#endif
593
594 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().left());
595 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().right());
596
597 // - sign as y -= ax <=> y = y-ax = -ax + y = axpy with -a
598 REAL32 ar = -( a.elem().elem().elem().elem());
599 REAL32* aptr = &ar;
600 if( s.hasOrderedRep() ) {
601
602 REAL32* xptr = (REAL32 *)&(x.elem(s.start()).elem(0).elem(0).real());
603 REAL32* yptr = &(d.elem(s.start()).elem(0).elem(0).real());
604
605 int total_n_3vec = (s.end()-s.start()+1);
606
607 ordered_sse_vaxOpy3_user_arg arg = {yptr, aptr, xptr, yptr, vaxpy3};
608
610
612 // Original code
614 //int n_3vec = (s.end()-s.start()+1)*24;
615 //vaxpy3(yptr, aptr, xptr, yptr, n_3vec);
616 }
617 else {
618 const int* tab = s.siteTable().slice();
619
620 int totalSize = s.numSiteTable();
621
622 unordered_sse_vaxOpy3_y_user_arg arg(x, d, aptr, tab, 1, vaxpy3);
623
625
627 // Original code
629 /*
630 for(int j=0; j < s.numSiteTable(); j++) {
631 int i = tab[j];
632 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
633 REAL32 *yptr = (REAL32 *) &(d.elem(i).elem(0).elem(0).real());
634
635 vaxpy3(yptr, aptr, xptr, yptr, 24);
636 }*/
637 }
638
639 }
640
641
642// Vec = Vec *Scalar + Vec (AXPY)
643template<>
644inline
646 const OpAssign &op,
647 const QDPExpr<
653 OLattice< TVec > > &rhs,
654 const Subset& s)
655{
656
657#ifdef DEBUG_BLAS
658 QDPIO::cout << "SSE: z = x*a + y" << endl;
659#endif
660
661 // Peel the stuff out of the expression
662 // y is the right side of rhs
663 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right());
664
665 // ax is the right side of rhs and is in a binary node
666 typedef BinaryNode<OpMultiply,
669
670 // get the binary node
671 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
672
673 // get a and x out of the bynary node
674 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.right());
675 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.left());
676 // Set pointers
677 REAL32 ar = a.elem().elem().elem().elem();
678 REAL32 *aptr = (REAL32 *)&ar;
679 if( s.hasOrderedRep() ) {
680 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
681 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
682 REAL32* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
683
684 int total_n_3vec = (s.end()-s.start()+1);
685
686 ordered_sse_vaxOpy3_user_arg arg = {zptr, aptr, xptr, yptr, vaxpy3};
687
689
691 // Original code
693 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
694 //int n_3vec = (s.end()-s.start()+1)*24;
695 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
696 }
697 else {
698 const int* tab = s.siteTable().slice();
699
700 int totalSize = s.numSiteTable();
701
702 unordered_sse_vaxOpy3_z_user_arg arg(x, y, d, aptr, tab, vaxpy3);
703
705
707 // Original code
709 /*
710 for(int j=0; j < s.numSiteTable(); j++) {
711 int i = tab[j];
712 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
713 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
714 REAL32* zptr = &(d.elem(i).elem(0).elem(0).real());
715 vaxpy3(zptr, aptr, xptr, yptr, 24);
716 }*/
717 }
718 }
719
720
721// Vec = Vec + Vec * Scalar (AXPY)
722template<>
723inline
725 const OpAssign &op,
726 const QDPExpr<
732 OLattice< TVec > > &rhs,
733 const Subset& s)
734{
735#ifdef DEBUG_BLAS
736 QDPIO::cout << "SSE: z = y + x*a" << endl;
737#endif
738
739
740 // Peel the stuff out of the expression
741
742 // y is the left side of rhs
743 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().left());
744
745 // ax is the right side of rhs and is in a binary node
746 typedef BinaryNode<OpMultiply,
749
750 // get the binary node
751 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
752
753 // get a and x out of the bynary node
754 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.right());
755 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.left());
756 // Set pointers
757 REAL32 ar = a.elem().elem().elem().elem();
758 REAL32 *aptr = (REAL32 *)&ar;
759 if( s.hasOrderedRep() ) {
760 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
761 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
762 REAL32* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
763
764 int total_n_3vec = (s.end()-s.start()+1);
765
766 ordered_sse_vaxOpy3_user_arg arg = {zptr, aptr, xptr, yptr, vaxpy3};
767
769
771 // Original code
773 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
774 //int n_3vec = (s.end()-s.start()+1)*24;
775 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
776 }
777 else {
778 const int* tab = s.siteTable().slice();
779
780 int totalSize = s.numSiteTable();
781
782 unordered_sse_vaxOpy3_z_user_arg arg(x, y, d, aptr, tab, vaxpy3);
783
785
787 // Original code
789 /*
790 for(int j=0; j < s.numSiteTable(); j++) {
791 int i = tab[j];
792 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
793 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
794 REAL32* zptr = &(d.elem(i).elem(0).elem(0).real());
795
796 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
797 vaxpy3(zptr, aptr, xptr, yptr, 24);
798 }*/
799 }
800}
801
802
803// Vec = Vec*Scalar - Vec (AXMY)
804template<>
805inline
807 const OpAssign &op,
808 const QDPExpr<
814 OLattice< TVec > > &rhs,
815 const Subset& s)
816{
817#ifdef DEBUG_BLAS
818 QDPIO::cout << "SSE: z = x*a - y" << endl;
819#endif
820
821 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right());
822
823
824 typedef BinaryNode<OpMultiply,
827
828 // get the binary node
829 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
830
831 // get a and x out of the bynary node
832 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.right());
833 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.left());
834 // Set pointers
835 REAL32 ar = a.elem().elem().elem().elem();
836 REAL32 *aptr = (REAL32 *)&ar;
837 if ( s.hasOrderedRep() ) {
838 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
839 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
840 REAL32* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
841
842 int total_n_3vec = (s.end()-s.start()+1);
843
844 ordered_sse_vaxOpy3_user_arg arg = {zptr, aptr, xptr, yptr, vaxmy3};
845
847
849 // Original code
851 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
852 //int n_3vec = (s.end()-s.start()+1)*24;
853 //vaxmy3(zptr, aptr, xptr, yptr, n_3vec);
854 }
855 else {
856 const int* tab = s.siteTable().slice();
857
858 int totalSize = s.numSiteTable();
859
860 unordered_sse_vaxOpy3_z_user_arg arg(x, y, d, aptr, tab, vaxmy3);
861
863
865 // Original code
867 /*
868 for(int j=0; j < s.numSiteTable(); j++) {
869 int i = tab[j];
870 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
871 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
872 REAL32* zptr = &(d.elem(i).elem(0).elem(0).real());
873
874
875 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
876 vaxmy3(zptr, aptr, xptr, yptr, 24);
877 }*/
878 }
879}
880
881
882// Vec = Vec - Vec*Scalar (AXPY with -Scalar)
883template<>
884inline
886 const OpAssign &op,
887 const QDPExpr<
893 OLattice< TVec > > &rhs,
894 const Subset& s)
895{
896#ifdef DEBUG_BLAS
897 QDPIO::cout << "SSE: z = y - x*a" << endl;
898#endif
899
900 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().left());
901
902 typedef BinaryNode<OpMultiply,
905
906 // get the binary node
907 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
908
909 // get a and x out of the bynary node
910 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.right());
911 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.left());
912 // Set pointers etc.
913
914 // -ve sign as y - ax = -ax + y = axpy with -a.
915 REAL32 ar = -a.elem().elem().elem().elem();
916 REAL32 *aptr = (REAL32 *)&ar;
917 if( s.hasOrderedRep() ) {
918 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
919 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
920 REAL32* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
921
922 int total_n_3vec = (s.end()-s.start()+1);
923
924 ordered_sse_vaxOpy3_user_arg arg = {zptr, aptr, xptr, yptr, vaxpy3};
925
927
929 // Original code
931 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
932 //int n_3vec = (s.end()-s.start()+1)*24;
933 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
934 }
935 else {
936 const int* tab = s.siteTable().slice();
937
938 int totalSize = s.numSiteTable();
939
940 unordered_sse_vaxOpy3_z_user_arg arg(x, y, d, aptr, tab, vaxpy3);
941
943
945 // Original code
947 /*
948 for(int j=0; j < s.numSiteTable(); j++) {
949 int i = tab[j];
950 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
951 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
952 REAL32* zptr = &(d.elem(i).elem(0).elem(0).real());
953
954 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
955 vaxpy3(zptr, aptr, xptr, yptr, 24);
956 }*/
957 }
958}
959
960
961template<>
962inline
964 const OpAssign &op,
965 const QDPExpr<
969 OLattice< TVec > > &rhs,
970 const Subset& s)
971{
972#ifdef DEBUG_BLAS
973 cout << "SSE: v+v " << endl;
974#endif
975
976 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(rhs.expression().left());
977 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(rhs.expression().right());
978
979 if( s.hasOrderedRep() ) {
980 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
981 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
982 REAL32* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
983
984 int total_n_3vec = (s.end()-s.start()+1);
985
986 ordered_sse_vOp_user_arg arg = {zptr, xptr, yptr, vadd};
987
989
991 // Original code
993 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
994 //int n_3vec = (s.end()-s.start()+1)*24;
995 //vadd(zptr, xptr, yptr, n_3vec);
996 }
997 else {
998 const int* tab = s.siteTable().slice();
999
1000 int totalSize = s.numSiteTable();
1001
1002 unordered_sse_vOp_z_user_arg arg(x, y, d, tab, vadd);
1003
1005
1007 // Original code
1009 /*
1010 for(int j=0; j < s.numSiteTable(); j++) {
1011 int i = tab[j];
1012 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
1013 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
1014 REAL32* zptr = &(d.elem(i).elem(0).elem(0).real());
1015
1016
1017 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1018 vadd(zptr, xptr, yptr, 24);
1019 }*/
1020 }
1021}
1022
1023template<>
1024inline
1026 const OpAssign &op,
1027 const QDPExpr<
1031 OLattice< TVec > > &rhs,
1032 const Subset& s)
1033{
1034#ifdef DEBUG_BLAS
1035 cout << "SSE: v-v " << endl;
1036#endif
1037
1038 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(rhs.expression().left());
1039 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(rhs.expression().right());
1040
1041 if( s.hasOrderedRep() ) {
1042 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
1043 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
1044 REAL32* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1045
1046 int total_n_3vec = (s.end()-s.start()+1);
1047
1048 ordered_sse_vOp_user_arg arg = {zptr, xptr, yptr, vsub};
1049
1051
1053 // Original code
1055 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1056 //int n_3vec = (s.end()-s.start()+1)*24;
1057 //vsub(zptr, xptr, yptr, n_3vec);
1058 }
1059 else {
1060 const int* tab = s.siteTable().slice();
1061
1062 int totalSize = s.numSiteTable();
1063
1064 unordered_sse_vOp_z_user_arg arg(x, y, d, tab, vsub);
1065
1067
1069 // Original code
1071 /*
1072 for(int j=0; j < s.numSiteTable(); j++) {
1073 int i = tab[j];
1074 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
1075 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
1076 REAL32* zptr = &(d.elem(i).elem(0).elem(0).real());
1077
1078
1079 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1080 vsub(zptr, xptr, yptr, 24);
1081 }*/
1082 }
1083}
1084
1085// Vec = Scal * Vec
1086template<>
1087inline
1089 const OpAssign &op,
1090 const QDPExpr<
1094 OLattice< TVec > > &rhs,
1095 const Subset& s)
1096{
1097#ifdef DEBUG_BLAS
1098 cout << "SSE: v = a*v " << endl;
1099#endif
1100 const OLattice< TVec > &x = static_cast<const OLattice< TVec >&>(rhs.expression().right());
1101 const OScalar< TScal > &a = static_cast<const OScalar< TScal >&>(rhs.expression().left());
1102
1103 REAL32 ar = a.elem().elem().elem().elem();
1104 REAL32 *aptr = &ar;
1105 if( s.hasOrderedRep() ) {
1106 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
1107 REAL32 *zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1108
1109 int total_n_3vec = (s.end()-s.start()+1);
1110
1111 ordered_sse_vscal_user_arg arg = {zptr, aptr, xptr};
1112
1114
1116 // Original code
1118 //int n_3vec = (s.end()-s.start()+1)*24;
1119
1120 //vscal(zptr, aptr, xptr, n_3vec);
1121 }
1122 else {
1123 const int* tab = s.siteTable().slice();
1124
1125 int totalSize = s.numSiteTable();
1126
1127 unordered_sse_vscal_user_arg arg(x, d, aptr, tab);
1128
1130
1132 // Original code
1134 /*
1135 for(int j=0; j < s.numSiteTable(); j++) {
1136 int i = tab[j];
1137 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
1138 REAL32 *zptr = &(d.elem(i).elem(0).elem(0).real());
1139
1140 vscal(zptr, aptr, xptr, 24);
1141 }*/
1142 }
1143}
1144
1145template<>
1146inline
1148 const OpAssign &op,
1149 const QDPExpr<
1153 OLattice< TVec > > &rhs,
1154 const Subset& s)
1155{
1156#ifdef DEBUG_BLAS
1157 cout << "SSE: v = v*a " << endl;
1158#endif
1159
1160 const OLattice< TVec > &x = static_cast<const OLattice< TVec >&>(rhs.expression().left());
1161 const OScalar< TScal > &a = static_cast<const OScalar< TScal >&>(rhs.expression().right());
1162
1163 REAL32 ar = a.elem().elem().elem().elem();
1164 REAL32 *aptr = &ar;
1165 if( s.hasOrderedRep() ) {
1166 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
1167 REAL32 *zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1168
1169 int total_n_3vec = (s.end()-s.start()+1);
1170
1171 ordered_sse_vscal_user_arg arg = {zptr, aptr, xptr};
1172
1174
1176 // Original code
1178 //int n_3vec = (s.end()-s.start()+1)*24;
1179
1180 //vscal(zptr, aptr, xptr, n_3vec);
1181 }
1182 else {
1183 const int* tab = s.siteTable().slice();
1184
1185 int totalSize = s.numSiteTable();
1186
1187 unordered_sse_vscal_user_arg arg(x, d, aptr, tab);
1188
1190
1192 // Original code
1194 /*
1195 for(int j=0; j < s.numSiteTable(); j++) {
1196 int i = tab[j];
1197 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
1198 REAL32 *zptr = &(d.elem(i).elem(0).elem(0).real());
1199
1200 vscal(zptr, aptr, xptr, 24);
1201 }*/
1202 }
1203}
1204
1205
1206//-----------------------------------------------------------------------------
1207// v *= a
1208template<>
1209inline
1211 const OpMultiplyAssign &op,
1212 const QDPExpr<
1215 OScalar< TScal > > &rhs,
1216 const Subset& s)
1217{
1218 const OScalar< TScal >& a = static_cast< const OScalar<TScal >&>(rhs.expression().child());
1219
1220
1221#ifdef DEBUG_BLAS
1222 QDPIO::cout << "BJ: v *= a, a = " << a << endl;
1223#endif
1224
1225 REAL32 ar = a.elem().elem().elem().elem();
1226 if( s.hasOrderedRep() ) {
1227 REAL32 * xptr = &(d.elem(s.start()).elem(0).elem(0).real());
1228 REAL32 * zptr = xptr;
1229
1230 int total_n_3vec = (s.end()-s.start()+1);
1231
1232 ordered_sse_vscal_user_arg arg = {zptr, &ar, xptr};
1233
1235
1237 // Original code
1239 //int n_3vec = (s.end()-s.start()+1)*24;
1240 //vscal(zptr,&ar, xptr, n_3vec);
1241 }
1242 else {
1243 const int* tab = s.siteTable().slice();
1244
1245 int totalSize = s.numSiteTable();
1246
1247 unordered_sse_vscal_user_arg arg(d, d, &ar, tab);
1248
1250
1252 // Original code
1254 /*
1255 for(int j=0; j < s.numSiteTable(); j++) {
1256 int i = tab[j];
1257 REAL32 * xptr = &(d.elem(i).elem(0).elem(0).real());
1258 REAL32 * zptr = xptr;
1259
1260 vscal(zptr, &ar, xptr, 24);
1261 }*/
1262 }
1263}
1264
1265// v /= a
1266template<>
1267inline
1269 const OpDivideAssign &op,
1270 const QDPExpr<
1273 OScalar< TScal > > &rhs,
1274 const Subset& s)
1275{
1276 const OScalar< TScal >& a = static_cast< const OScalar<TScal >&>(rhs.expression().child());
1277
1278
1279#ifdef DEBUG_BLAS
1280 QDPIO::cout << "BJ: v /= a, a = " << a << endl;
1281#endif
1282
1283 REAL32 ar = (REAL)1/a.elem().elem().elem().elem();
1284 if( s.hasOrderedRep() ) {
1285 REAL32 * xptr = &(d.elem(s.start()).elem(0).elem(0).real());
1286 REAL32 * zptr = xptr;
1287
1288 int total_n_3vec = (s.end()-s.start()+1);
1289
1290 ordered_sse_vscal_user_arg arg = {zptr, &ar, xptr};
1291
1293
1295 // Original code
1297 //int n_3vec = (s.end()-s.start()+1)*24;
1298 //vscal(zptr,&ar, xptr, n_3vec);
1299 }
1300 else {
1301 const int* tab = s.siteTable().slice();
1302
1303 int totalSize = s.numSiteTable();
1304
1305 unordered_sse_vscal_user_arg arg(d, d, &ar, tab);
1306
1308
1310 // Original code
1312 /*
1313 for(int j=0; j < s.numSiteTable(); j++) {
1314 int i = tab[j];
1315 REAL32 * xptr = &(d.elem(i).elem(0).elem(0).real());
1316 REAL32 * zptr = xptr;
1317
1318 vscal(zptr,&ar, xptr, 24);
1319 }*/
1320 }
1321}
1322
1323// v += v
1324template<>
1325inline
1327 const OpAddAssign &op,
1328 const QDPExpr<
1331 OLattice< TVec > > &rhs,
1332 const Subset& s)
1333{
1334 const OLattice< TVec >& x = static_cast< const OLattice<TVec >&>(rhs.expression().child());
1335
1336
1337
1338#ifdef DEBUG_BLAS
1339 QDPIO::cout << "BJ: v += v" << endl;
1340#endif
1341 if(s.hasOrderedRep() ) {
1342 REAL32 *xptr = (REAL32 *)(&x.elem(s.start()).elem(0).elem(0).real());
1343 REAL32 *yptr = (REAL32 *)(&d.elem(s.start()).elem(0).elem(0).real());
1344
1345 int total_n_3vec = (s.end()-s.start()+1);
1346
1347 ordered_sse_vOp_user_arg arg = {yptr, yptr, xptr, vadd};
1348
1350
1352 // Original code
1354 /*int n_3vec = (s.end() - s.start()+1)*24;
1355 REAL32 *xptr = (REAL32 *)(&x.elem(s.start()).elem(0).elem(0).real());
1356 REAL32 *yptr = (REAL32 *)(&d.elem(s.start()).elem(0).elem(0).real());
1357 REAL one = 1;
1358 vadd(yptr, yptr, xptr,n_3vec);*/
1359 }
1360 else {
1361 const int* tab = s.siteTable().slice();
1362
1363 int totalSize = s.numSiteTable();
1364
1365 unordered_sse_vOp_y_user_arg arg(x, d, tab, vadd);
1366
1368
1370 // Original code
1372 /*
1373 for(int j=0; j < s.numSiteTable(); j++) {
1374 int i = tab[j];
1375 REAL32 *xptr = (REAL32 *)(&x.elem(i).elem(0).elem(0).real());
1376 REAL32 *yptr = (REAL32 *)(&d.elem(i).elem(0).elem(0).real());
1377 vadd(yptr, yptr, xptr,24);
1378 }*/
1379 }
1380}
1381
1382// v -= v
1383template<>
1384inline
1386 const OpSubtractAssign &op,
1387 const QDPExpr<
1390 OLattice< TVec > > &rhs,
1391 const Subset& s)
1392{
1393 const OLattice< TVec >& x = static_cast< const OLattice<TVec >&>(rhs.expression().child());
1394
1395
1396
1397#ifdef DEBUG_BLAS
1398 QDPIO::cout << "BJ: v -= v" << endl;
1399#endif
1400 if( s.hasOrderedRep() ) {
1401 //int n_3vec = (s.end() - s.start()+1)*24;
1402 REAL32 *xptr = (REAL32 *)(&x.elem(s.start()).elem(0).elem(0).real());
1403 REAL32 *yptr = (REAL32 *)(&d.elem(s.start()).elem(0).elem(0).real());
1404
1405 int total_n_3vec = (s.end()-s.start()+1);
1406
1407 ordered_sse_vOp_user_arg arg = {yptr, yptr, xptr, vsub};
1408
1410
1412 // Original code
1414 //vsub(yptr, yptr, xptr, n_3vec);
1415 }
1416 else {
1417 const int* tab = s.siteTable().slice();
1418
1419 int totalSize = s.numSiteTable();
1420
1421 unordered_sse_vOp_y_user_arg arg(x, d, tab, vsub);
1422
1424
1426 // Original code
1428 /*
1429 for(int j=0; j < s.numSiteTable(); j++) {
1430 int i = tab[j];
1431 REAL32 *xptr = (REAL32 *)(&x.elem(i).elem(0).elem(0).real());
1432 REAL32 *yptr = (REAL32 *)(&d.elem(i).elem(0).elem(0).real());
1433
1434 vsub(yptr, yptr, xptr, 24);
1435 }*/
1436 }
1437}
1438
1439
1440// z = ax + by
1441template<>
1442inline
1444 const OpAssign &op,
1445 const QDPExpr<
1453 OLattice< TVec > > &rhs,
1454 const Subset& s)
1455{
1456
1457#ifdef DEBUG_BLAS
1458 QDPIO::cout << "z = a*x + b*y" << endl;
1459#endif
1460
1461 // Peel the stuff out of the expression
1462 // y is the right side of rhs
1463
1464 // ax is the left side of rhs and is in a binary node
1465 typedef BinaryNode<OpMultiply,
1468
1469 // get the binary node
1470 const BN &mulNode1 = static_cast<const BN&> (rhs.expression().left());
1471 const BN &mulNode2 = static_cast<const BN&> (rhs.expression().right());
1472
1473 // get a and x out of the binary node
1474 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
1475 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
1476
1477 // get b and y out of the binary node
1478 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1479 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right());
1480
1481
1482 // Set pointers
1483 REAL32 *aptr = (REAL32 *)&(a.elem().elem().elem().elem());
1484 REAL32 *bptr = (REAL32 *)&(b.elem().elem().elem().elem());
1485 if( s.hasOrderedRep() ) {
1486 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
1487 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
1488 REAL32 * zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1489
1490 int total_n_3vec = (s.end()-s.start()+1);
1491
1492 ordered_sse_vaxOpby3_user_arg arg = {zptr, aptr, xptr, bptr, yptr, vaxpby3};
1493
1495
1497 // Original code
1499 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1500 //int n_3vec = (s.end()-s.start()+1)*24;
1501 //vaxpby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1502 }
1503 else {
1504 const int* tab = s.siteTable().slice();
1505
1506 int totalSize = s.numSiteTable();
1507
1508 unordered_sse_vaxOpby3_user_arg arg(aptr, x, bptr, y, d, tab, vaxpby3);
1509
1511
1513 // Original code
1515 /*
1516 for(int j=0; j < s.numSiteTable(); j++) {
1517 int i = tab[j];
1518 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
1519 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
1520 REAL32 * zptr = &(d.elem(i).elem(0).elem(0).real());
1521
1522 vaxpby3(zptr, aptr, xptr, bptr, yptr, 24);
1523 }*/
1524 }
1525}
1526
1527
1528// z = xa + by
1529template<>
1530inline
1532 const OpAssign &op,
1533 const QDPExpr<
1541 OLattice< TVec > > &rhs,
1542 const Subset& s)
1543{
1544
1545#ifdef DEBUG_BLAS
1546 QDPIO::cout << "z = x*a + b*y" << endl;
1547#endif
1548
1549 // Peel the stuff out of the expression
1550 // y is the right side of rhs
1551
1552 // ax is the left side of rhs and is in a binary node
1553 typedef BinaryNode<OpMultiply,
1556
1557 typedef BinaryNode<OpMultiply,
1560
1561
1562
1563 // get the binary node
1564 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
1565 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
1566
1567 // get a and x out of the binary node
1568 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.left());
1569
1570 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.right());
1571
1572 // get b and y out of the binary node
1573 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1574 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right());
1575
1576
1577 // Set pointers
1578 REAL32 *aptr = (REAL32 *)&(a.elem().elem().elem().elem());
1579 REAL32 *bptr = (REAL32 *)&(b.elem().elem().elem().elem());
1580 if( s.hasOrderedRep() ) {
1581 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
1582 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
1583 REAL32 * zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1584
1585 int total_n_3vec = (s.end()-s.start()+1);
1586
1587 ordered_sse_vaxOpby3_user_arg arg = {zptr, aptr, xptr, bptr, yptr, vaxpby3};
1588
1590
1592 // Original code
1594 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1595 //int n_3vec = (s.end()-s.start()+1)*24;
1596 //vaxpby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1597 }
1598 else {
1599 const int* tab = s.siteTable().slice();
1600
1601 int totalSize = s.numSiteTable();
1602
1603 unordered_sse_vaxOpby3_user_arg arg(aptr, x, bptr, y, d, tab, vaxpby3);
1604
1606
1608 // Original code
1610 /*
1611 for(int j=0; j < s.numSiteTable(); j++) {
1612 int i = tab[j];
1613 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
1614 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
1615 REAL32 * zptr = &(d.elem(i).elem(0).elem(0).real());
1616
1617
1618 vaxpby3(zptr, aptr, xptr, bptr, yptr, 24);
1619 }*/
1620 }
1621}
1622
1623// z = ax + yb
1624template<>
1625inline
1627 const OpAssign &op,
1628 const QDPExpr<
1636 OLattice< TVec > > &rhs,
1637 const Subset& s)
1638{
1639
1640#ifdef DEBUG_BLAS
1641 QDPIO::cout << "z = a*x + y*b" << endl;
1642#endif
1643
1644 // Peel the stuff out of the expression
1645 // y is the right side of rhs
1646
1647 // type of a*x
1648 typedef BinaryNode<OpMultiply,
1651
1652 // type of y*b
1653 typedef BinaryNode<OpMultiply,
1656
1657
1658
1659 // get the binary nodes
1660 // a*x node
1661 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
1662
1663 // y*b node
1664 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
1665
1666 // get a and x out of the binary node
1667 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
1668
1669 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
1670
1671
1672 // get b and y out of the binary node
1673 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.left());
1674
1675 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.right());
1676
1677
1678 // Set pointers
1679 REAL32 *aptr = (REAL32 *)&(a.elem().elem().elem().elem());
1680 REAL32 *bptr = (REAL32 *)&(b.elem().elem().elem().elem());
1681 if( s.hasOrderedRep() ) {
1682
1683 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
1684 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
1685 REAL32 * zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1686
1687 int total_n_3vec = (s.end()-s.start()+1);
1688
1689 ordered_sse_vaxOpby3_user_arg arg = {zptr, aptr, xptr, bptr, yptr, vaxpby3};
1690
1692
1694 // Original code
1696 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1697 //int n_3vec = (s.end()-s.start()+1)*24;
1698 //vaxpby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1699 }
1700 else {
1701 const int* tab = s.siteTable().slice();
1702
1703 int totalSize = s.numSiteTable();
1704
1705 unordered_sse_vaxOpby3_user_arg arg(aptr, x, bptr, y, d, tab, vaxpby3);
1706
1708
1710 // Original code
1712 /*
1713 for(int j=0; j < s.numSiteTable(); j++) {
1714 int i = tab[j];
1715 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
1716 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
1717 REAL32 * zptr = &(d.elem(i).elem(0).elem(0).real());
1718
1719
1720 vaxpby3(zptr, aptr, xptr, bptr, yptr, 24);
1721 }*/
1722 }
1723}
1724
1725// z = xa + yb
1726template<>
1727inline
1729 const OpAssign &op,
1730 const QDPExpr<
1738 OLattice< TVec > > &rhs,
1739 const Subset& s)
1740{
1741
1742#ifdef DEBUG_BLAS
1743 QDPIO::cout << "z = x*a + y*b" << endl;
1744#endif
1745
1746 // Peel the stuff out of the expression
1747 // y is the right side of rhs
1748
1749 // ax is the left side of rhs and is in a binary node
1750 typedef BinaryNode<OpMultiply,
1753
1754 // get the binary node
1755 const BN &mulNode1 = static_cast<const BN&> (rhs.expression().left());
1756 const BN &mulNode2 = static_cast<const BN&> (rhs.expression().right());
1757
1758 // get a and x out of the binary node
1759 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.left());
1760 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.right());
1761
1762 // get b and y out of the binary node
1763 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.left());
1764
1765 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.right());
1766
1767 // Set pointers
1768 REAL32 *aptr = (REAL32 *)&(a.elem().elem().elem().elem());
1769 REAL32 *bptr = (REAL32 *)&(b.elem().elem().elem().elem());
1770 if( s.hasOrderedRep() ) {
1771 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
1772 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
1773 REAL32 * zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1774
1775 int total_n_3vec = (s.end()-s.start()+1);
1776
1777 ordered_sse_vaxOpby3_user_arg arg = {zptr, aptr, xptr, bptr, yptr, vaxpby3};
1778
1780
1782 // Original code
1784 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1785 //int n_3vec = (s.end()-s.start()+1)*24;
1786 //vaxpby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1787 }
1788 else {
1789 const int* tab = s.siteTable().slice();
1790
1791 int totalSize = s.numSiteTable();
1792
1793 unordered_sse_vaxOpby3_user_arg arg(aptr, x, bptr, y, d, tab, vaxpby3);
1794
1796
1798 // Original code
1800 /*
1801 for(int j=0; j < s.numSiteTable(); j++) {
1802 int i = tab[j];
1803 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
1804 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
1805 REAL32 * zptr = &(d.elem(i).elem(0).elem(0).real());
1806
1807
1808 // Get the no of 3vecs. i and s.end() are inclusive so add +1
1809 vaxpby3(zptr, aptr, xptr, bptr, yptr, 24);
1810 }*/
1811 }
1812}
1813
1814// z = ax - by
1815template<>
1816inline
1818 const OpAssign &op,
1819 const QDPExpr<
1827 OLattice< TVec > > &rhs,
1828 const Subset& s)
1829{
1830
1831#ifdef DEBUG_BLAS
1832 QDPIO::cout << "z = a*x - b*y" << endl;
1833#endif
1834
1835 // Peel the stuff out of the expression
1836 // y is the right side of rhs
1837
1838 // ax is the left side of rhs and is in a binary node
1839 typedef BinaryNode<OpMultiply,
1842
1843 // get the binary node
1844 const BN &mulNode1 = static_cast<const BN&> (rhs.expression().left());
1845 const BN &mulNode2 = static_cast<const BN&> (rhs.expression().right());
1846
1847 // get a and x out of the binary node
1848 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
1849 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
1850
1851 // get b and y out of the binary node
1852 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1853 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right());
1854
1855
1856 // Set pointers
1857 REAL32 *aptr = (REAL32 *)&(a.elem().elem().elem().elem());
1858 REAL32 *bptr = (REAL32 *)&(b.elem().elem().elem().elem());
1859 if( s.hasOrderedRep() ) {
1860 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
1861 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
1862 REAL32 * zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1863
1864 int total_n_3vec = (s.end()-s.start()+1);
1865
1866 ordered_sse_vaxOpby3_user_arg arg = {zptr, aptr, xptr, bptr, yptr, vaxmby3};
1867
1869
1871 // Original code
1873 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1874 //int n_3vec = (s.end()-s.start()+1)*24;
1875 //vaxmby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1876 }
1877 else {
1878 const int* tab = s.siteTable().slice();
1879
1880 int totalSize = s.numSiteTable();
1881
1882 unordered_sse_vaxOpby3_user_arg arg(aptr, x, bptr, y, d, tab, vaxmby3);
1883
1885
1887 // Original code
1889 /*
1890 for(int j=0; j < s.numSiteTable(); j++) {
1891 int i = tab[j];
1892 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
1893 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
1894 REAL32 * zptr = &(d.elem(i).elem(0).elem(0).real());
1895
1896
1897 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1898 vaxmby3(zptr, aptr, xptr, bptr, yptr, 24);
1899
1900 }*/
1901 }
1902}
1903
1904
1905// z = xa - by
1906template<>
1907inline
1909 const OpAssign &op,
1910 const QDPExpr<
1918 OLattice< TVec > > &rhs,
1919 const Subset& s)
1920{
1921
1922#ifdef DEBUG_BLAS
1923 QDPIO::cout << "z = x*a - b*y" << endl;
1924#endif
1925
1926 // Peel the stuff out of the expression
1927 // y is the right side of rhs
1928
1929 // ax is the left side of rhs and is in a binary node
1930 typedef BinaryNode<OpMultiply,
1933
1934 typedef BinaryNode<OpMultiply,
1937
1938
1939
1940 // get the binary node
1941 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
1942 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
1943
1944 // get a and x out of the binary node
1945 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.left());
1946
1947 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.right());
1948
1949 // get b and y out of the binary node
1950 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1951 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right());
1952
1953
1954 // Set pointers
1955 REAL32 *aptr = (REAL32 *)&(a.elem().elem().elem().elem());
1956 REAL32 *bptr = (REAL32 *)&(b.elem().elem().elem().elem());
1957 if( s.hasOrderedRep() ) {
1958 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
1959 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
1960 REAL32 * zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1961
1962 int total_n_3vec = (s.end()-s.start()+1);
1963
1964 ordered_sse_vaxOpby3_user_arg arg = {zptr, aptr, xptr, bptr, yptr, vaxmby3};
1965
1967
1969 // Original code
1971 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1972 //int n_3vec = (s.end()-s.start()+1)*24;
1973 //vaxmby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1974 }
1975 else {
1976 const int* tab = s.siteTable().slice();
1977
1978 int totalSize = s.numSiteTable();
1979
1980 unordered_sse_vaxOpby3_user_arg arg(aptr, x, bptr, y, d, tab, vaxmby3);
1981
1983
1985 // Original code
1987 /*
1988 for(int j=0; j < s.numSiteTable(); j++) {
1989 int i = tab[j];
1990 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
1991 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
1992 REAL32 * zptr = &(d.elem(i).elem(0).elem(0).real());
1993
1994
1995 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1996 vaxmby3(zptr, aptr, xptr, bptr, yptr, 24);
1997 }*/
1998 }
1999}
2000
2001// z = ax - yb
2002template<>
2003inline
2005 const OpAssign &op,
2006 const QDPExpr<
2014 OLattice< TVec > > &rhs,
2015 const Subset& s)
2016{
2017
2018#ifdef DEBUG_BLAS
2019 QDPIO::cout << "z = a*x - y*b" << endl;
2020#endif
2021
2022 // Peel the stuff out of the expression
2023 // y is the right side of rhs
2024
2025 // type of a*x
2026 typedef BinaryNode<OpMultiply,
2029
2030 // type of y*b
2031 typedef BinaryNode<OpMultiply,
2034
2035
2036
2037 // get the binary nodes
2038 // a*x node
2039 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
2040
2041 // y*b node
2042 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
2043
2044 // get a and x out of the binary node
2045 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
2046
2047 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
2048
2049
2050 // get b and y out of the binary node
2051 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.left());
2052
2053 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.right());
2054
2055
2056 // Set pointers
2057 REAL32 *aptr = (REAL32 *)&(a.elem().elem().elem().elem());
2058 REAL32 *bptr = (REAL32 *)&(b.elem().elem().elem().elem());
2059 if( s.hasOrderedRep() ) {
2060 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
2061 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
2062 REAL32 * zptr = &(d.elem(s.start()).elem(0).elem(0).real());
2063
2064 int total_n_3vec = (s.end()-s.start()+1);
2065
2066 ordered_sse_vaxOpby3_user_arg arg = {zptr, aptr, xptr, bptr, yptr, vaxmby3};
2067
2069
2071 // Original code
2073 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2074 //int n_3vec = (s.end()-s.start()+1)*24;
2075 //vaxmby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
2076 }
2077 else {
2078 const int* tab = s.siteTable().slice();
2079
2080 int totalSize = s.numSiteTable();
2081
2082 unordered_sse_vaxOpby3_user_arg arg(aptr, x, bptr, y, d, tab, vaxmby3);
2083
2085
2087 // Original code
2089 /*
2090 for(int j=0; j < s.numSiteTable(); j++) {
2091 int i = tab[j];
2092 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
2093 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
2094 REAL32 * zptr = &(d.elem(i).elem(0).elem(0).real());
2095
2096
2097 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2098 vaxmby3(zptr, aptr, xptr, bptr, yptr, 24);
2099 }*/
2100 }
2101}
2102
2103// z = xa - yb
2104template<>
2105inline
2107 const OpAssign &op,
2108 const QDPExpr<
2116 OLattice< TVec > > &rhs,
2117 const Subset& s)
2118{
2119
2120#ifdef DEBUG_BLAS
2121 QDPIO::cout << "z = x*a - y*b" << endl;
2122#endif
2123
2124 // Peel the stuff out of the expression
2125 // y is the right side of rhs
2126
2127 // ax is the left side of rhs and is in a binary node
2128 typedef BinaryNode<OpMultiply,
2131
2132 // get the binary node
2133 const BN &mulNode1 = static_cast<const BN&> (rhs.expression().left());
2134 const BN &mulNode2 = static_cast<const BN&> (rhs.expression().right());
2135
2136 // get a and x out of the binary node
2137 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.left());
2138 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.right());
2139
2140 // get b and y out of the binary node
2141 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.left());
2142
2143 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.right());
2144
2145 // Set pointers
2146 REAL32 *aptr = (REAL32 *)&(a.elem().elem().elem().elem());
2147 REAL32 *bptr = (REAL32 *)&(b.elem().elem().elem().elem());
2148 if( s.hasOrderedRep() ) {
2149 REAL32 *xptr = (REAL32 *) &(x.elem(s.start()).elem(0).elem(0).real());
2150 REAL32 *yptr = (REAL32 *) &(y.elem(s.start()).elem(0).elem(0).real());
2151 REAL32 *zptr = &(d.elem(s.start()).elem(0).elem(0).real());
2152
2153 int total_n_3vec = (s.end()-s.start()+1);
2154
2155 ordered_sse_vaxOpby3_user_arg arg = {zptr, aptr, xptr, bptr, yptr, vaxmby3};
2156
2158
2160 // Original code
2162 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2163 //int n_3vec = (s.end()-s.start()+1)*24;
2164 //vaxmby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
2165 }
2166 else {
2167 const int* tab = s.siteTable().slice();
2168
2169 int totalSize = s.numSiteTable();
2170
2171 unordered_sse_vaxOpby3_user_arg arg(aptr, x, bptr, y, d, tab, vaxmby3);
2172
2174
2176 // Original code
2178 /*
2179 for(int j=0; j < s.numSiteTable(); j++) {
2180 int i = tab[j];
2181 REAL32 *xptr = (REAL32 *) &(x.elem(i).elem(0).elem(0).real());
2182 REAL32 *yptr = (REAL32 *) &(y.elem(i).elem(0).elem(0).real());
2183 REAL32 *zptr = &(d.elem(i).elem(0).elem(0).real());
2184
2185
2186 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2187 vaxmby3(zptr, aptr, xptr, bptr, yptr, 24);
2188 }*/
2189 }
2190}
2191
2192
2193
2194//-----------------------------------------------------------------------------
2195
2196#endif // if defined(QDP_SCALARSITE_USE_EVALUATE)
2197
2198
2199
2200#if 1
2201// Global norm squared of a vector...
2202template<>
2205{
2206#ifdef DEBUG_BLAS
2207 QDPIO::cout << "Using SSE sumsq" << endl;
2208#endif
2209
2210 if ( s.hasOrderedRep() ) {
2211
2212#ifdef DEBUG_BLAS
2213 QDPIO::cout << "BJ sumsq " << endl;
2214#endif
2215 int n_real = (s.end() - s.start() + 1);
2217 arg.vptr = (REAL32*) &(s1.elem(s.start()).elem(0).elem(0).real());
2219 arg.func = local_sumsq_24_48;
2221
2222 REAL64 ltmp=arg.results[0];
2223 for(int i=1; i < qdpNumThreads(); i++) {
2224 ltmp += arg.results[i];
2225 }
2226 // Use specialized sum for REAL64
2229 return lsum;
2230 }
2231 else {
2232 REAL64 ltmp1=0;
2233
2234
2235 const int* tab=s.siteTable().slice();
2236 for(int j=0; j < s.numSiteTable(); j++) {
2237 int i=tab[j];
2238 REAL64 ltmp2=0;
2239 REAL32* s1ptr = (REAL32 *)&(s1.elem(i).elem(0).elem(0).real());
2240 local_sumsq_24_48(&ltmp2, s1ptr, 1);
2241 ltmp1 += ltmp2;
2242 }
2243
2244 // Use specialized sum for REAL64
2247
2248 return lsum;
2249
2250 }
2251}
2252
2253template<>
2256{
2257#ifdef DEBUG_BLAS
2258 QDPIO::cout << "Using SSE sumsq all" << endl;
2259#endif
2260
2261 int n_real = (all.end() - all.start() + 1);
2263 arg.vptr = (REAL32*) &(s1.elem(all.start()).elem(0).elem(0).real());
2265 arg.func = local_sumsq_24_48;
2267
2268
2269 // I am relying on this being a Double here
2270 REAL64 ltmp=arg.results[0];
2271 for(int i=1; i < qdpNumThreads(); i++) {
2272 ltmp += arg.results[i];
2273 }
2274 // Do the sum with specialized QMP_sum_double
2277 return lsum;
2278}
2279#endif
2280
2281
2282template<>
2283inline BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t
2285 const QDPType< TVec, OLattice<TVec> > &v2)
2286{
2287#ifdef DEBUG_BLAS
2288 QDPIO::cout << "BJ: innerProduct all" << endl;
2289#endif
2290
2291 // This BinaryReturn has Type_t
2292 // OScalar<OScalar<OScalar<RComplex<PScalar<REAL> > > > >
2294 // Inner product is accumulated internally in REAL64
2295 REAL64 ip[2];
2296 ip[0]=0;
2297 ip[1]=0;
2298
2299 // Length of subset
2300 unsigned long n_3vec = (all.end() - all.start() + 1)*Ns;
2301
2302 // Call My CDOT
2303 local_vcdot((REAL64 *)&(ip[0]), (REAL64 *)&(ip[1]),
2304 (REAL32 *)&(v1.elem(all.start()).elem(0).elem(0).real()),
2305 (REAL32 *)&(v2.elem(all.start()).elem(0).elem(0).real()),
2306 (int)n_3vec);
2307
2308
2309 // Global sum -- still on a vector of doubles
2311
2312 // Downcast (and possibly lose precision) here
2313 lprod.elem().elem().elem().real() = ip[0];
2314 lprod.elem().elem().elem().imag() = ip[1];
2315
2316 // Return
2317 return lprod;
2318}
2319
2320template<>
2321inline BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t
2323 const QDPType< TVec, OLattice<TVec> > &v2,
2324 const Subset& s)
2325{
2326
2328
2329 REAL64 ip[2];
2330 if( s.hasOrderedRep() ) {
2331#ifdef DEBUG_BLAS
2332 QDPIO::cout << "BJ: innerProduct s" << endl;
2333#endif
2334
2335 // This BinaryReturn has Type_t
2336 // OScalar<OScalar<OScalar<RComplex<PScalar<REAL> > > > >
2337
2338
2339 ip[0] = 0;
2340 ip[1] = 0;
2341
2342 unsigned long n_3vec = (s.end() - s.start() + 1)*Ns;
2343 local_vcdot((REAL64 *)&(ip[0]), (REAL64 *)&(ip[1]),
2344 (REAL32 *)&(v1.elem(s.start()).elem(0).elem(0).real()),
2345 (REAL32 *)&(v2.elem(s.start()).elem(0).elem(0).real()),
2346 (int)n_3vec);
2347
2348
2349 }
2350 else {
2351 // This BinaryReturn has Type_t
2352 // OScalar<OScalar<OScalar<RComplex<PScalar<REAL> > > > >
2353
2354 REAL64 ip_tmp[2];
2355 ip[0] = 0;
2356 ip[1] = 0;
2357
2358 const int* tab = s.siteTable().slice();
2359 for(int j=0; j < s.numSiteTable(); j++) {
2360 int i=tab[j];
2361
2362 local_vcdot((REAL64 *)&(ip_tmp[0]), (REAL64 *)&(ip_tmp[1]),
2363 (REAL32 *)&(v1.elem(i).elem(0).elem(0).real()),
2364 (REAL32 *)&(v2.elem(i).elem(0).elem(0).real()),
2365 (int)Ns);
2366
2367 ip[0] += ip_tmp[0];
2368 ip[1] += ip_tmp[1];
2369 }
2370 }
2371
2373
2374 lprod.elem().elem().elem().real() = ip[0];
2375 lprod.elem().elem().elem().imag() = ip[1];
2376
2377 return lprod;
2378}
2379
2380
2381
2382// Inner Product Real
2383template<>
2384inline
2385BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t
2387 const QDPType< TVec, OLattice<TVec> > &v2)
2388{
2389#ifdef DEBUG_BLAS
2390 QDPIO::cout << "BJ: innerProductReal all" << endl;
2391#endif
2392
2393 // This BinaryReturn has Type_t
2394 // OScalar<OScalar<OScalar<RScalar<PScalar<REAL> > > > >
2396 // Inner product is accumulated internally in REAL64
2397 REAL64 ip_re=0;
2398
2399 // Length of subset
2400 unsigned long n_3vec = (all.end() - all.start() + 1)*Ns;
2401
2402 // Call My CDOT
2403 local_vcdot_real((DOUBLE*)&ip_re,
2404 (REAL32 *)&(v1.elem(all.start()).elem(0).elem(0).real()),
2405 (REAL32 *)&(v2.elem(all.start()).elem(0).elem(0).real()),
2406 (int)n_3vec);
2407
2408 // Global sum
2410
2411 // Whether CDOT did anything or not ip_re and ip_im should
2412 // now be right. Assign them to the ReturnType
2413 lprod.elem().elem().elem().elem() = ip_re;
2414
2415
2416 // Return
2417 return lprod;
2418}
2419
2420
2421template<>
2422inline
2423BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t
2425 const QDPType< TVec, OLattice<TVec> > &v2,
2426 const Subset& s)
2427{
2428
2430
2431 REAL64 ip_re=0;
2432
2433 if( s.hasOrderedRep() ) {
2434#ifdef DEBUG_BLAS
2435 QDPIO::cout << "BJ: innerProductReal s" << endl;
2436#endif
2437
2438 // This BinaryReturn has Type_t
2439 // OScalar<OScalar<OScalar<RScalar<PScalar<REAL> > > > >
2440 unsigned long n_3vec = (s.end() - s.start() + 1)*Ns;
2441 local_vcdot_real((REAL64 *)&ip_re,
2442 (REAL32 *)&(v1.elem(s.start()).elem(0).elem(0).real()),
2443 (REAL32 *)&(v2.elem(s.start()).elem(0).elem(0).real()),
2444 (int)n_3vec);
2445
2446
2447
2448
2449 }
2450 else {
2451 REAL64 ip_re_tmp=0;
2452 const int* tab = s.siteTable().slice();
2453 for(int j=0; j < s.numSiteTable(); j++) {
2454 int i=tab[j];
2455
2456 local_vcdot_real((REAL64 *)&ip_re_tmp,
2457 (REAL32 *)&(v1.elem(i).elem(0).elem(0).real()),
2458 (REAL32 *)&(v2.elem(i).elem(0).elem(0).real()),
2459 (int)Ns);
2460
2461 ip_re += ip_re_tmp;
2462
2463 }
2464
2465
2466 }
2467
2469 lprod.elem().elem().elem().elem() = ip_re;
2470 return lprod;
2471}
2472
2473#if 1
2474template<>
2477{
2478#ifdef DEBUG_BLAS
2479 QDPIO::cout << "Using SSE multi1d sumsq all" << endl;
2480#endif
2481
2482 int n_real = (all.end() - all.start() + 1);
2483 REAL64 ltmp = 0;
2484 for(int n=0; n < s1.size(); ++n)
2485 {
2486 const REAL32 *s1ptr = &(s1[n].elem(all.start()).elem(0).elem(0).real());
2487
2488 // I am relying on this being a Double here
2489 REAL64 lltmp=0;
2490 local_sumsq_24_48(&lltmp, (REAL32 *)s1ptr, n_real);
2491
2492 ltmp += lltmp;
2493 }
2494
2497
2498 return lsum;
2499}
2500
2501
2502template<>
2504 norm2(const multi1d< OLattice< TVec > >& s1, const Subset& s)
2505{
2506#ifdef DEBUG_BLAS
2507 QDPIO::cout << "Using SSE multi1d sumsq all" << endl;
2508#endif
2509
2510 REAL64 ltmp = 0;
2511
2512 if( s.hasOrderedRep() ) {
2513 int n_real = (s.end() - s.start() + 1);
2514
2515 for(int n=0; n < s1.size(); ++n) {
2516
2517 const REAL32 *s1ptr = &(s1[n].elem(s.start()).elem(0).elem(0).real());
2518
2519 // I am relying on this being a Double here
2520 REAL64 lltmp=0;
2521 local_sumsq_24_48(&lltmp, (REAL32 *)s1ptr, n_real);
2522
2523 ltmp += lltmp;
2524 }
2525 }
2526 else {
2527
2528 const int* tab = s.siteTable().slice();
2529 for(int n=0; n < s1.size(); ++n) {
2530 for(int j=0; j < s.numSiteTable(); j++) {
2531 int i=tab[j];
2532 REAL64 lltmp=0;
2533 const REAL32 *s1ptr = &(s1[n].elem(i).elem(0).elem(0).real());
2534 local_sumsq_24_48(&lltmp,(REAL32 *)s1ptr,1);
2535 ltmp += lltmp;
2536 }
2537 }
2538
2539 }
2542
2543 return lsum;
2544}
2545#endif
2546
2547template<>
2548inline BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t
2550 const multi1d< OLattice<TVec> > &v2)
2551{
2552#ifdef DEBUG_BLAS
2553 QDPIO::cout << "BJ: multi1d innerProduct all" << endl;
2554#endif
2555
2556 // This BinaryReturn has Type_t
2557 // OScalar<OScalar<OScalar<RComplex<PScalar<REAL> > > > >
2559 // Inner product is accumulated internally in REAL64
2560 REAL64 ip[2];
2561 ip[0]=0;
2562 ip[1]=0;
2563
2564 // Length of subset
2565 unsigned long n_3vec = (all.end() - all.start() + 1)*Ns;
2566
2567 for(int n=0; n < v1.size(); ++n)
2568 {
2569 REAL64 iip[2];
2570 iip[0]=0;
2571 iip[1]=0;
2572
2573 // Call My CDOT
2574 local_vcdot((REAL64 *)&(iip[0]), (REAL64 *)&(iip[1]),
2575 (REAL32 *)&(v1[n].elem(all.start()).elem(0).elem(0).real()),
2576 (REAL32 *)&(v2[n].elem(all.start()).elem(0).elem(0).real()),
2577 (int)n_3vec);
2578
2579 ip[0] += iip[0];
2580 ip[1] += iip[1];
2581 }
2582
2583 // Global sum -- still on a vector of doubles
2585
2586 // Downcast (and possibly lose precision) here
2587 lprod.elem().elem().elem().real() = ip[0];
2588 lprod.elem().elem().elem().imag() = ip[1];
2589
2590 // Return
2591 return lprod;
2592}
2593
2594template<>
2595inline BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t
2597 const multi1d< OLattice<TVec> > &v2,
2598 const Subset& s)
2599{
2600#ifdef DEBUG_BLAS
2601 QDPIO::cout << "BJ: multi1d innerProduct subset" << endl;
2602#endif
2603
2604
2605 // This BinaryReturn has Type_t
2606 // OScalar<OScalar<OScalar<RComplex<PScalar<REAL> > > > >
2608 // Inner product is accumulated internally in REAL64
2609 REAL64 ip[2];
2610 ip[0]=0;
2611 ip[1]=0;
2612
2613 if( s.hasOrderedRep() ) {
2614 // Length of subset
2615 unsigned long n_3vec = (s.end() - s.start() + 1)*Ns;
2616
2617 for(int n=0; n < v1.size(); ++n) {
2618
2619 REAL64 iip[2];
2620 iip[0]=0;
2621 iip[1]=0;
2622
2623 // Call My CDOT
2624 local_vcdot((REAL64 *)&(iip[0]),
2625 (REAL64 *)&(iip[1]),
2626 (REAL32 *)&(v1[n].elem(s.start()).elem(0).elem(0).real()),
2627 (REAL32 *)&(v2[n].elem(s.start()).elem(0).elem(0).real()),
2628 (int)n_3vec);
2629
2630 ip[0] += iip[0];
2631 ip[1] += iip[1];
2632 }
2633 }
2634 else {
2635 // Length of an atom
2636 unsigned long n_3vec = Ns;
2637
2638 // Site table
2639 const int* tab = s.siteTable().slice();
2640
2641 // Loop through N5
2642 for(int n=0; n < v1.size(); ++n) {
2643
2644 // Loop through site table
2645 for(int j=0; j < s.numSiteTable(); j++) {
2646 int i=tab[j];
2647 REAL64 iip[2];
2648 iip[0]=0;
2649 iip[1]=0;
2650
2651 // Call My CDOT
2652 local_vcdot((REAL64 *)&(iip[0]), (REAL64 *)&(iip[1]),
2653 (REAL32 *)&(v1[n].elem(i).elem(0).elem(0).real()),
2654 (REAL32 *)&(v2[n].elem(i).elem(0).elem(0).real()),
2655 (int)n_3vec);
2656
2657 ip[0] += iip[0];
2658 ip[1] += iip[1];
2659 }
2660 }
2661 }
2662 // Global sum -- still on a vector of doubles
2664
2665 // Downcast (and possibly lose precision) here
2666 lprod.elem().elem().elem().real() = ip[0];
2667 lprod.elem().elem().elem().imag() = ip[1];
2668
2669 // Return
2670 return lprod;
2671}
2672
2673// Inner Product Real
2674template<>
2675inline
2676BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t
2678 const multi1d< OLattice<TVec> > &v2)
2679{
2680#ifdef DEBUG_BLAS
2681 QDPIO::cout << "BJ: innerProductReal(multi1d) all" << endl;
2682#endif
2683
2684 // This BinaryReturn hasType_t
2685 // OScalar<OScalar<OScalar<RScalar<PScalar<REAL> > > > >
2687 // Inner product is accumulated internally in REAL64
2688 REAL64 ip_re=0;
2689
2690 // Length of subset
2691 unsigned long n_3vec = (all.end() - all.start() + 1)*Ns;
2692
2693 for(int n=0; n < v1.size(); ++n)
2694 {
2695 REAL64 iip_re=0;
2696
2697 // Call My CDOT
2698 local_vcdot_real((REAL64 *)&iip_re,
2699 (REAL32 *)&(v1[n].elem(all.start()).elem(0).elem(0).real()),
2700 (REAL32 *)&(v2[n].elem(all.start()).elem(0).elem(0).real()),
2701 (int)n_3vec);
2702
2703 ip_re += iip_re;
2704 }
2705
2706 // Global sum
2708
2709 // Whether CDOT did anything or not ip_re and ip_im should
2710 // now be right. Assign them to the ReturnType
2711 lprod.elem().elem().elem().elem() = ip_re;
2712
2713
2714 // Return
2715 return lprod;
2716}
2717
2718// Inner Product Real
2719template<>
2720inline
2721BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t
2723 const multi1d< OLattice<TVec> > &v2,
2724 const Subset& s)
2725{
2726#ifdef DEBUG_BLAS
2727 QDPIO::cout << "BJ: innerProductReal(multi1d) all" << endl;
2728#endif
2729
2730 // This BinaryReturn hasType_t
2731 // OScalar<OScalar<OScalar<RScalar<PScalar<REAL> > > > >
2733 // Inner product is accumulated internally in REAL64
2734 REAL64 ip_re=0;
2735
2736 // Length of subset
2737 if ( s.hasOrderedRep() ) {
2738 unsigned long n_3vec = (s.end() - s.start() + 1)*Ns;
2739
2740 for(int n=0; n < v1.size(); ++n) {
2741
2742 REAL64 iip_re=0;
2743
2744 // Call My CDOT
2745 local_vcdot_real((REAL64 *)&iip_re,
2746 (REAL32 *)&(v1[n].elem(s.start()).elem(0).elem(0).real()),
2747 (REAL32 *)&(v2[n].elem(s.start()).elem(0).elem(0).real()),
2748 (int)n_3vec);
2749
2750 ip_re += iip_re;
2751 }
2752 }
2753 else {
2754 unsigned long n_3vec = Ns;
2755 const int* tab = s.siteTable().slice();
2756
2757 for(int n=0; n < v1.size(); ++n) {
2758 for(int j=0; j< s.numSiteTable(); j++) {
2759 const int i=tab[j];
2760 REAL64 iip_re=0;
2761
2762 // Call My CDOT
2763 local_vcdot_real((REAL64 *)&iip_re,
2764 (REAL32 *)&(v1[n].elem(i).elem(0).elem(0).real()),
2765 (REAL32 *)&(v2[n].elem(i).elem(0).elem(0).real()),
2766 (int)n_3vec);
2767
2768 ip_re += iip_re;
2769 }
2770 }
2771 }
2772 // Global sum
2774
2775 // Whether CDOT did anything or not ip_re and ip_im should
2776 // now be right. Assign them to the ReturnType
2777 lprod.elem().elem().elem().elem() = ip_re;
2778
2779
2780 // Return
2781 return lprod;
2782}
2783
2784
2785
2786
2787#if defined(QDP_SCALARSITE_DEBUG)
2788#undef QDP_SCALARSITE_DEBUG
2789#endif
2790
2791#if defined(QDP_SCALARSITE_USE_EVALUATE)
2792#undef QDP_SCALARSITE_USE_EVALUATE
2793#endif
2794
2795
2796} // namespace QDP;
2797
2798#endif // guard
2799
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
Primitive Scalar.
Primitive spin Vector class.
Expression class for QDP.
Definition qdp_qdpexpr.h:16
QDPType - major type class/container for all QDP objects.
Definition qdp_qdptype.h:29
Subsets - controls how lattices are looped.
Definition qdp_subset.h:39
int end() const
Definition qdp_subset.h:81
const multi1d< int > & siteTable() const
Definition qdp_subset.h:83
int numSiteTable() const
Definition qdp_subset.h:84
bool hasOrderedRep() const
Definition qdp_subset.h:79
int start() const
Definition qdp_subset.h:80
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
REAL32 REAL
double REAL64
REAL64 DOUBLE
BinaryReturn< C1, C2, FnInnerProductReal >::Type_t innerProductReal(const QDPType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
OScalar = innerProductReal(adj(source1)*source2).
UnaryReturn< C, FnNorm2 >::Type_t norm2(const QDPType< T, C > &s1)
OScalar = norm2(trace(adj(source)*source)).
BinaryReturn< C1, C2, FnInnerProduct >::Type_t innerProduct(const QDPType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
OScalar = innerProduct(adj(source1)*source2).
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)
const int Ns
Definition qdp_params.h:26
Subset all
Default all subset.
Definition qdp_subset.cc:16
StandardOutputStream cout
Definition qdp_stdio.cc:21
void globalSum(T &dest)
Sum across all nodes.
void globalSumArray(unsigned int *dest, int len)
Wrapper to get a functional unsigned global sum.
Yet another random number generator.
void local_sumsq_24_48(REAL64 *Out, REAL32 *In, int n_3vec)
PScalar< PScalar< RScalar< REAL > > > TScal
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
void dispatch_to_threads(int numSiteTable, Arg a, void(*func)(int, int, int, Arg *))
void unordered_sse_vOp_y_evaluate_function(int lo, int hi, int myId, unordered_sse_vOp_y_user_arg *a)
void local_vcdot_real(REAL64 *Out_re, REAL32 *V1, REAL32 *V2, int n_3vec)
void unordered_sse_vaxOpy3_y_evaluate_function(int lo, int hi, int myId, unordered_sse_vaxOpy3_y_user_arg *a)
void local_vcdot(REAL64 *Out_re, REAL64 *Out_im, REAL32 *V1, REAL32 *V2, int n_3vec)
void unordered_sse_vscal_evaluate_function(int lo, int hi, int myId, unordered_sse_vscal_user_arg *a)
void vaxpby3(REAL *Out, REAL *ap, REAL *xp, REAL *bp, REAL *yp, int n_3vec)
void ordered_sse_vaxOpy3_evaluate_function(int lo, int hi, int myId, ordered_sse_vaxOpy3_user_arg *a)
void vsub(REAL *Out, REAL *In1, REAL *In2, int n_3vec)
void ordered_norm_single_func(int lo, int hi, int myId, ordered_sse_norm_single_user_arg *a)
void vaxpy3(REAL *Out, REAL *scalep, REAL *InScale, REAL *Add, int n_4vec)
void vadd(REAL *Out, REAL *In1, REAL *In2, int n_3vec)
void vaxmy3(REAL *Out, REAL *scalep, REAL *InScale, REAL *Sub, int n_3vec)
void vaxmby3(REAL *Out, REAL *ap, REAL *xp, REAL *bp, REAL *yp, int n_3vec)
void unordered_vOp_z_evaluate_function(int lo, int hi, int myId, unordered_sse_vOp_z_user_arg *a)
void ordered_sse_vOp_evaluate_function(int lo, int hi, int myId, ordered_sse_vOp_user_arg *a)
void vscal(REAL *Out, REAL *scalep, REAL *In, int n_3vec)
PSpinVector< PColorVector< RComplex< REAL >, 3 >, Ns > TVec
void unordered_sse_vaxOpby3_evaluate_function(int lo, int hi, int myId, unordered_sse_vaxOpby3_user_arg *arg)
void ordered_sse_vscal_evaluate_function(int lo, int hi, int myId, ordered_sse_vscal_user_arg *a)
void unordered_sse_vaxOpy3_z_evaluate_function(int lo, int hi, int myId, unordered_sse_vaxOpy3_z_user_arg *a)
int qdpNumThreads()
void ordered_sse_vaxOpby3_evaluate_function(int lo, int hi, int myId, ordered_sse_vaxOpby3_user_arg *arg)
void unordered_sse_vaxOpy3_z_evaluate_function(int lo, int hi, int myId, unordered_sse_vaxOpy3_z_user_arg *a)
void ordered_sse_vaxOpy3_evaluate_function(int lo, int hi, int myId, ordered_sse_vaxOpy3_user_arg *a)
void(* func)(REAL64 *, REAL32 *, int)