QDP++
qdp_scalarsite_bagel_qdp_blas_g5.h
Go to the documentation of this file.
1// $Id: qdp_scalarsite_bagel_qdp_blas_g5.h,v 1.8 2007-06-10 14:32:10 edwards Exp $
2
8
9
10#ifndef QDP_SCALARSITE_BAGEL_QDP_BLAS_G5_H
11#define QDP_SCALARSITE_BAGEL_QDP_BLAS_G5_H
12
13#include "bagel_qdp.h"
14
15
16namespace QDP {
17
18// Types needed for the expression templates.
19// TVec has outer Ns template so it ought to work for staggered as well
20typedef PSpinVector<PColorVector<RComplex<REAL>, 3>, Ns> TVec;
21typedef PScalar<PScalar<RScalar<REAL> > > TScal;
22
23// #define DEBUG_BLAS_G6
24// TVec is the LatticeFermion from qdp_dwdefs.h with the OLattice<> stripped
25// from around it
26
27// TScalar is the usual Real, with the OScalar<> stripped from it
28//
29// THis is simply to make the code more readable, and reduces < < s and > >s
30// in the template arguments
31
32// d += Scalar*ChiralProjPlus(Vec);
33template<>
34inline
35void evaluate(OLattice< TVec >& d,
36 const OpAddAssign& op,
37 const QDPExpr<BinaryNode<OpMultiply,
38 Reference< QDPType< TScal, OScalar < TScal > > >,
39 UnaryNode< FnChiralProjectPlus, Reference< QDPType<TVec,OLattice<TVec> > > >
40>,
41 OLattice< TVec > > &rhs,
42 const Subset& s)
43{
44
45#ifdef DEBUG_BLAS_G5
46 QDPIO::cout << "y += a*P{+}x" << endl;
47#endif
48
49 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().right().child());
50 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().left());
51
52
53 REAL ar = a.elem().elem().elem().elem();
54 REAL* aptr = &ar;
55
56 if( s.hasOrderedRep() ) {
57
58 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
59 REAL* yptr = (REAL *)&(d.elem(s.start()).elem(0).elem(0).real());
60
61
62 int n_4vec = (s.end()-s.start()+1);
63 qdp_vxpay_chp(yptr, aptr,yptr, xptr, n_4vec);
64 }
65 else {
66 const int* tab = s.siteTable().slice();
67 for(int j = 0; j < s.numSiteTable() ; j++) {
68 int i = tab[j];
69
70 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
71 REAL* yptr = (REAL *)&(d.elem(i).elem(0).elem(0).real());
72 qdp_vxpay_chp(yptr, aptr,yptr, xptr, 1);
73 }
74 }
75
76}
77
78// d += Scalar*ChiralProjMinus(Vec);
79template<>
80inline
81void evaluate(OLattice< TVec >& d,
82 const OpAddAssign& op,
83 const QDPExpr<BinaryNode<OpMultiply,
84 Reference< QDPType< TScal, OScalar < TScal > > >,
85 UnaryNode< FnChiralProjectMinus, Reference< QDPType<TVec,OLattice<TVec> > > >
86>,
87 OLattice< TVec > > &rhs,
88 const Subset& s)
89{
90
91#ifdef DEBUG_BLAS_G5
92 QDPIO::cout << "y += a*P{-}x" << endl;
93#endif
94
95 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().right().child());
96 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().left());
97
98
99 REAL ar = a.elem().elem().elem().elem();
100 REAL* aptr = &ar;
101
102 if( s.hasOrderedRep() ) {
103 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
104 REAL* yptr = (REAL *)&(d.elem(s.start()).elem(0).elem(0).real());
105
106
107 int n_4vec = (s.end()-s.start()+1);
108 qdp_vxpay_chm(yptr, aptr,yptr, xptr, n_4vec);
109 }
110 else {
111 const int* tab = s.siteTable().slice();
112 for(int j = 0; j < s.numSiteTable() ; j++) {
113 int i = tab[j];
114
115 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
116 REAL* yptr = (REAL *)&(d.elem(i).elem(0).elem(0).real());
117
118 qdp_vxpay_chm(yptr, aptr,yptr, xptr, 1);
119 }
120 }
121
122
123}
124
125
126// d -= Scalar*ChiralProjPlus(Vec);
127template<>
128inline
129void evaluate(OLattice< TVec >& d,
130 const OpSubtractAssign& op,
131 const QDPExpr<BinaryNode<OpMultiply,
132 Reference< QDPType< TScal, OScalar < TScal > > >,
133 UnaryNode< FnChiralProjectPlus, Reference< QDPType<TVec,OLattice<TVec> > > >
134>,
135 OLattice< TVec > > &rhs,
136 const Subset& s)
137{
138
139#ifdef DEBUG_BLAS_G5
140 QDPIO::cout << "y -= a*P{+}x" << endl;
141#endif
142
143 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().right().child());
144 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().left());
145
146
147 REAL ar = a.elem().elem().elem().elem();
148 REAL* aptr = &ar;
149 if( s.hasOrderedRep() ) {
150 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
151 REAL* yptr = (REAL *)&(d.elem(s.start()).elem(0).elem(0).real());
152
153 int n_4vec = (s.end()-s.start()+1);
154 qdp_vxmay_chp(yptr, aptr,yptr, xptr, n_4vec);
155 }
156 else {
157 const int* tab = s.siteTable().slice();
158 for(int j = 0; j < s.numSiteTable() ; j++) {
159 int i = tab[j];
160
161 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
162 REAL* yptr = (REAL *)&(d.elem(i).elem(0).elem(0).real());
163
164 qdp_vxmay_chp(yptr, aptr,yptr, xptr, 1);
165 }
166 }
167
168}
169
170// d -= Scalar*ChiralProjMinus(Vec);
171template<>
172inline
173void evaluate(OLattice< TVec >& d,
174 const OpSubtractAssign& op,
175 const QDPExpr<BinaryNode<OpMultiply,
176 Reference< QDPType< TScal, OScalar < TScal > > >,
177 UnaryNode< FnChiralProjectMinus, Reference< QDPType<TVec,OLattice<TVec> > > >
178>,
179 OLattice< TVec > > &rhs,
180 const Subset& s)
181{
182
183#ifdef DEBUG_BLAS_G5
184 QDPIO::cout << "y -= a*P{-}x" << endl;
185#endif
186
187 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().right().child());
188 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().left());
189
190
191 REAL ar = a.elem().elem().elem().elem();
192 REAL* aptr = &ar;
193
194 if( s.hasOrderedRep() ) {
195 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
196 REAL* yptr = (REAL *)&(d.elem(s.start()).elem(0).elem(0).real());
197
198
199 int n_4vec = (s.end()-s.start()+1);
200 qdp_vxmay_chm(yptr, aptr,yptr, xptr, n_4vec);
201 }
202 else {
203 const int* tab = s.siteTable().slice();
204 for(int j = 0; j < s.numSiteTable() ; j++) {
205 int i = tab[j];
206
207 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
208 REAL* yptr = (REAL *)&(d.elem(i).elem(0).elem(0).real());
209
210 qdp_vxmay_chm(yptr, aptr,yptr, xptr, 1);
211 }
212 }
213
214
215}
216
217
218// d += ChiralProjPlus(Vec);
219template<>
220inline
221void evaluate(OLattice< TVec >& d,
222 const OpAddAssign& op,
223 const QDPExpr<
224 UnaryNode< FnChiralProjectPlus, Reference< QDPType<TVec,OLattice<TVec> > > >,
225 OLattice< TVec > > &rhs,
226 const Subset& s)
227{
228
229#ifdef DEBUG_BLAS_G5
230 QDPIO::cout << "y += P{+}x" << endl;
231#endif
232
233 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().child());
234
235 if ( s.hasOrderedRep() ) {
236 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
237 REAL* yptr = (REAL *)&(d.elem(s.start()).elem(0).elem(0).real());
238
239
240 int n_4vec = (s.end()-s.start()+1);
241 qdp_vadd_chp(yptr, yptr, xptr, n_4vec);
242 }
243 else {
244 const int* tab = s.siteTable().slice();
245 for(int j = 0; j < s.numSiteTable() ; j++) {
246 int i = tab[j];
247
248 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
249 REAL* yptr = (REAL *)&(d.elem(i).elem(0).elem(0).real());
250
251 qdp_vadd_chp(yptr, yptr, xptr, 1);
252 }
253 }
254
255
256}
257
258
259// d += ChiralProjMinus(Vec);
260template<>
261inline
262void evaluate(OLattice< TVec >& d,
263 const OpAddAssign& op,
264 const QDPExpr<
265 UnaryNode< FnChiralProjectMinus, Reference< QDPType<TVec,OLattice<TVec> > > >,
266 OLattice< TVec > > &rhs,
267 const Subset& s)
268{
269
270#ifdef DEBUG_BLAS_G5
271 QDPIO::cout << "y += P{-}x" << endl;
272#endif
273
274 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().child());
275
276 if( s.hasOrderedRep() ) {
277 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
278 REAL* yptr = (REAL *)&(d.elem(s.start()).elem(0).elem(0).real());
279
280
281 int n_4vec = (s.end()-s.start()+1);
282 qdp_vadd_chm(yptr, yptr, xptr, n_4vec);
283 }
284 else {
285 const int* tab = s.siteTable().slice();
286 for(int j = 0; j < s.numSiteTable() ; j++) {
287 int i = tab[j];
288
289 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
290 REAL* yptr = (REAL *)&(d.elem(i).elem(0).elem(0).real());
291 qdp_vadd_chm(yptr, yptr, xptr, 1);
292
293 }
294 }
295
296}
297
298
299// d -= ChiralProjPlus(Vec);
300template<>
301inline
302void evaluate(OLattice< TVec >& d,
303 const OpSubtractAssign& op,
304 const QDPExpr<
305 UnaryNode< FnChiralProjectPlus, Reference< QDPType<TVec,OLattice<TVec> > > >,
306 OLattice< TVec > > &rhs,
307 const Subset& s)
308{
309
310#ifdef DEBUG_BLAS_G5
311 QDPIO::cout << "y -= P{+}x" << endl;
312#endif
313
314 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().child());
315
316
317 if( s.hasOrderedRep() ) {
318 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
319 REAL* yptr = (REAL *)&(d.elem(s.start()).elem(0).elem(0).real());
320
321
322 int n_4vec = (s.end()-s.start()+1);
323 qdp_vsub_chp(yptr, yptr, xptr, n_4vec);
324 }
325 else {
326
327 const int* tab=s.siteTable().slice();
328 for(int j=0; j < s.numSiteTable(); j++) {
329 int i = tab[j];
330
331 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
332 REAL* yptr = (REAL *)&(d.elem(i).elem(0).elem(0).real());
333
334 qdp_vsub_chp(yptr, yptr, xptr, 1);
335 }
336
337 }
338
339}
340
341
342// d += ChiralProjMinus(Vec);
343template<>
344inline
345void evaluate(OLattice< TVec >& d,
346 const OpSubtractAssign& op,
347 const QDPExpr<
348 UnaryNode< FnChiralProjectMinus, Reference< QDPType<TVec,OLattice<TVec> > > >,
349 OLattice< TVec > > &rhs,
350 const Subset& s)
351{
352
353#ifdef DEBUG_BLAS_G5
354 QDPIO::cout << "y -= P{-}x" << endl;
355#endif
356
357 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().child());
358
359 if( s.hasOrderedRep() ) {
360 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
361 REAL* yptr = (REAL *)&(d.elem(s.start()).elem(0).elem(0).real());
362
363
364 int n_4vec = (s.end()-s.start()+1);
365 qdp_vsub_chm(yptr, yptr, xptr, n_4vec);
366 }
367 else {
368 const int* tab = s.siteTable().slice();
369 for(int j = 0; j < s.numSiteTable() ; j++) {
370 int i = tab[j];
371
372 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
373 REAL* yptr = (REAL *)&(d.elem(i).elem(0).elem(0).real());
374 qdp_vsub_chm(yptr, yptr, xptr, 1);
375
376 }
377 }
378
379
380}
381
382// d = x + a P_{+} y
383template<>
384inline
385void evaluate( OLattice< TVec > &d,
386 const OpAssign &op,
387 const QDPExpr<
388 BinaryNode<OpAdd,
389 Reference< QDPType< TVec, OLattice< TVec > > >,
390 BinaryNode<OpMultiply,
391 Reference< QDPType< TScal, OScalar< TScal > > >,
392 UnaryNode< FnChiralProjectPlus, Reference< QDPType< TVec, OLattice<TVec> > > >
393 >
394 >,
395 OLattice< TVec >
396 > &rhs,
397 const Subset& s)
398{
399#ifdef DEBUG_BLAS_G5
400 QDPIO::cout << "z = x + a*P{+} y" << endl;
401#endif
402
403
404 // Peel the stuff out of the expression
405
406 // y is the left side of rhs
407 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&> (rhs.expression().left());
408
409 // ax is the right side of rhs and is in a binary node
410 typedef BinaryNode<OpMultiply,
411 Reference< QDPType< TScal, OScalar< TScal > > >,
412 UnaryNode< FnChiralProjectPlus, Reference< QDPType< TVec, OLattice< TVec > > > > > BN;
413
414 // get the binary node
415 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
416
417 // get a and x out of the bynary node
418 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
419 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode.right().child());
420 // Set pointers
421 REAL ar = a.elem().elem().elem().elem();
422 REAL *aptr = (REAL *)&ar;
423 if( s.hasOrderedRep() ) {
424 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
425 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
426 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
427
428 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
429 int n_4vec = (s.end()-s.start()+1);
430 qdp_vxpay_chp(zptr, aptr, xptr, yptr, n_4vec);
431 }
432 else {
433 const int* tab = s.siteTable().slice();
434 for(int j = 0; j < s.numSiteTable() ; j++) {
435 int i = tab[j];
436
437 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
438 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
439 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
440
441 qdp_vxpay_chp(zptr, aptr, xptr, yptr, 1);
442
443
444 }
445 }
446
447}
448
449// d = x + a P_{-} y
450template<>
451inline
452void evaluate( OLattice< TVec > &d,
453 const OpAssign &op,
454 const QDPExpr<
455 BinaryNode<OpAdd,
456 Reference< QDPType< TVec, OLattice< TVec > > >,
457 BinaryNode<OpMultiply,
458 Reference< QDPType< TScal, OScalar< TScal > > >,
459 UnaryNode< FnChiralProjectMinus, Reference< QDPType< TVec, OLattice<TVec> > > >
460 >
461 >,
462 OLattice< TVec >
463 > &rhs,
464 const Subset& s)
465{
466#ifdef DEBUG_BLAS_G5
467 QDPIO::cout << "z = x + a*P{-} y" << endl;
468#endif
469
470
471 // Peel the stuff out of the expression
472
473 // y is the left side of rhs
474 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&> (rhs.expression().left());
475
476 // ax is the right side of rhs and is in a binary node
477 typedef BinaryNode<OpMultiply,
478 Reference< QDPType< TScal, OScalar< TScal > > >,
479 UnaryNode< FnChiralProjectMinus, Reference< QDPType< TVec, OLattice< TVec > > > > > BN;
480
481 // get the binary node
482 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
483
484 // get a and x out of the bynary node
485 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
486 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode.right().child());
487 // Set pointers
488 REAL ar = a.elem().elem().elem().elem();
489 REAL *aptr = (REAL *)&ar;
490
491 if( s.hasOrderedRep() ) {
492 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
493 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
494 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
495
496 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
497 int n_4vec = (s.end()-s.start()+1);
498 qdp_vxpay_chm(zptr, aptr, xptr, yptr, n_4vec);
499 }
500 else {
501 const int* tab = s.siteTable().slice();
502 for(int j = 0; j < s.numSiteTable() ; j++) {
503 int i = tab[j];
504
505 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
506 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
507 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
508
509 qdp_vxpay_chm(zptr, aptr, xptr, yptr, 1);
510 }
511 }
512
513}
514
515// d = x - a P_{+} y
516template<>
517inline
518void evaluate( OLattice< TVec > &d,
519 const OpAssign &op,
520 const QDPExpr<
521 BinaryNode<OpSubtract,
522 Reference< QDPType< TVec, OLattice< TVec > > >,
523 BinaryNode<OpMultiply,
524 Reference< QDPType< TScal, OScalar< TScal > > >,
525 UnaryNode< FnChiralProjectPlus, Reference< QDPType< TVec, OLattice<TVec> > > >
526 >
527 >,
528 OLattice< TVec >
529 > &rhs,
530 const Subset& s)
531{
532#ifdef DEBUG_BLAS_G5
533 QDPIO::cout << "z = x - a*P{+} y" << endl;
534#endif
535
536
537 // Peel the stuff out of the expression
538
539 // y is the left side of rhs
540 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&> (rhs.expression().left());
541
542 // ax is the right side of rhs and is in a binary node
543 typedef BinaryNode<OpMultiply,
544 Reference< QDPType< TScal, OScalar< TScal > > >,
545 UnaryNode< FnChiralProjectPlus, Reference< QDPType< TVec, OLattice< TVec > > > > > BN;
546
547 // get the binary node
548 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
549
550 // get a and x out of the bynary node
551 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
552 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode.right().child());
553 // Set pointers
554 REAL ar = a.elem().elem().elem().elem();
555 REAL *aptr = (REAL *)&ar;
556
557 if( s.hasOrderedRep() ) {
558 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
559 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
560 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
561
562 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
563 int n_4vec = (s.end()-s.start()+1);
564 qdp_vxmay_chp(zptr, aptr, xptr, yptr, n_4vec);
565 }
566 else {
567 const int* tab = s.siteTable().slice();
568 for(int j = 0; j < s.numSiteTable() ; j++) {
569 int i = tab[j];
570
571 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
572 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
573 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
574
575 qdp_vxmay_chp(zptr, aptr, xptr, yptr, 1);
576 }
577 }
578
579}
580
581// d = x - a P_{-} y
582template<>
583inline
584void evaluate( OLattice< TVec > &d,
585 const OpAssign &op,
586 const QDPExpr<
587 BinaryNode<OpSubtract,
588 Reference< QDPType< TVec, OLattice< TVec > > >,
589 BinaryNode<OpMultiply,
590 Reference< QDPType< TScal, OScalar< TScal > > >,
591 UnaryNode< FnChiralProjectMinus, Reference< QDPType< TVec, OLattice<TVec> > > >
592 >
593 >,
594 OLattice< TVec >
595 > &rhs,
596 const Subset& s)
597{
598#ifdef DEBUG_BLAS_G5
599 QDPIO::cout << "z = x - a*P{-} y" << endl;
600#endif
601
602
603 // Peel the stuff out of the expression
604
605 // y is the left side of rhs
606 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&> (rhs.expression().left());
607
608 // ax is the right side of rhs and is in a binary node
609 typedef BinaryNode<OpMultiply,
610 Reference< QDPType< TScal, OScalar< TScal > > >,
611 UnaryNode< FnChiralProjectMinus, Reference< QDPType< TVec, OLattice< TVec > > > > > BN;
612
613 // get the binary node
614 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
615
616 // get a and x out of the bynary node
617 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
618 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode.right().child());
619 // Set pointers
620 REAL ar = a.elem().elem().elem().elem();
621 REAL *aptr = (REAL *)&ar;
622
623 if( s.hasOrderedRep() ) {
624 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
625 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
626 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
627
628 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
629 int n_4vec = (s.end()-s.start()+1);
630 qdp_vxmay_chm(zptr, aptr, xptr, yptr, n_4vec);
631 }
632 else {
633 const int* tab = s.siteTable().slice();
634 for(int j = 0; j < s.numSiteTable() ; j++) {
635 int i = tab[j];
636
637 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
638 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
639 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
640
641 qdp_vxmay_chm(zptr, aptr, xptr, yptr, 1);
642 }
643 }
644}
645
646// d = ax + P+ y
647template<>
648inline
649void evaluate( OLattice< TVec > &d,
650 const OpAssign &op,
651 const QDPExpr<
652 BinaryNode<OpAdd,
653 BinaryNode<OpMultiply,
654 Reference< QDPType< TScal, OScalar< TScal > > >,
655 Reference< QDPType< TVec, OLattice< TVec > > >
656 >,
657 UnaryNode<FnChiralProjectPlus, Reference< QDPType<TVec, OLattice<TVec> > > >
658 >,
659 OLattice< TVec > > &rhs,
660 const Subset& s)
661{
662#ifdef DEBUG_BLAS_G5
663 QDPIO::cout << "z = a*x + P{+} y" << endl;
664#endif
665
666
667 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right().child());
668
669 // ax is the left side of rhs and is in a binary node
670 typedef BinaryNode<OpMultiply,
671 Reference< QDPType< TScal, OScalar< TScal > > >,
672 Reference< QDPType< TVec, OLattice< TVec > > > > BN;
673
674 // get the binary node
675 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
676
677 // get a and x out of the bynary node
678 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
679 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
680 // Set pointers
681 REAL ar = a.elem().elem().elem().elem();
682 REAL *aptr = (REAL *)&ar;
683
684 if( s.hasOrderedRep() ) {
685 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
686 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
687 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
688
689
690 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
691 int n_4vec = (s.end()-s.start()+1);
692 qdp_vaxpy_chp(zptr, aptr, xptr, yptr, n_4vec);
693 }
694 else {
695 const int* tab = s.siteTable().slice();
696 for(int j = 0; j < s.numSiteTable() ; j++) {
697 int i = tab[j];
698
699 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
700 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
701 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
702
703 qdp_vaxpy_chp(zptr, aptr, xptr, yptr, 1);
704 }
705 }
706
707}
708
709// d = ax + P- y
710template<>
711inline
712void evaluate( OLattice< TVec > &d,
713 const OpAssign &op,
714 const QDPExpr<
715 BinaryNode<OpAdd,
716 BinaryNode<OpMultiply,
717 Reference< QDPType< TScal, OScalar< TScal > > >,
718 Reference< QDPType< TVec, OLattice< TVec > > >
719 >,
720 UnaryNode<FnChiralProjectMinus, Reference< QDPType<TVec, OLattice<TVec> > > >
721 >,
722 OLattice< TVec > > &rhs,
723 const Subset& s)
724{
725#ifdef DEBUG_BLAS_G5
726 QDPIO::cout << "z = a*x + P{-} y" << endl;
727#endif
728
729
730 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right().child());
731
732 // ax is the left side of rhs and is in a binary node
733 typedef BinaryNode<OpMultiply,
734 Reference< QDPType< TScal, OScalar< TScal > > >,
735 Reference< QDPType< TVec, OLattice< TVec > > > > BN;
736
737 // get the binary node
738 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
739
740 // get a and x out of the bynary node
741 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
742 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
743 // Set pointers
744 REAL ar = a.elem().elem().elem().elem();
745 REAL *aptr = (REAL *)&ar;
746
747 if( s.hasOrderedRep() ) {
748 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
749 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
750 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
751
752
753 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
754 int n_4vec = (s.end()-s.start()+1);
755 qdp_vaxpy_chm(zptr, aptr, xptr, yptr, n_4vec);
756 }
757 else {
758 const int* tab = s.siteTable().slice();
759 for(int j = 0; j < s.numSiteTable() ; j++) {
760 int i = tab[j];
761
762 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
763 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
764 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
765
766 qdp_vaxpy_chm(zptr, aptr, xptr, yptr, 1);
767 }
768 }
769
770}
771
772
773// d = ax - P+ y
774template<>
775inline
776void evaluate( OLattice< TVec > &d,
777 const OpAssign &op,
778 const QDPExpr<
779 BinaryNode<OpSubtract,
780 BinaryNode<OpMultiply,
781 Reference< QDPType< TScal, OScalar< TScal > > >,
782 Reference< QDPType< TVec, OLattice< TVec > > >
783 >,
784 UnaryNode<FnChiralProjectPlus, Reference< QDPType<TVec, OLattice<TVec> > > >
785 >,
786 OLattice< TVec > > &rhs,
787 const Subset& s)
788{
789#ifdef DEBUG_BLAS_G5
790 QDPIO::cout << "z = a*x + P{+} y" << endl;
791#endif
792
793
794 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right().child());
795
796 // ax is the left side of rhs and is in a binary node
797 typedef BinaryNode<OpMultiply,
798 Reference< QDPType< TScal, OScalar< TScal > > >,
799 Reference< QDPType< TVec, OLattice< TVec > > > > BN;
800
801 // get the binary node
802 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
803
804 // get a and x out of the bynary node
805 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
806 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
807 // Set pointers
808 REAL ar = a.elem().elem().elem().elem();
809 REAL *aptr = (REAL *)&ar;
810
811 if( s.hasOrderedRep()) {
812 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
813 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
814 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
815
816
817 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
818 int n_4vec = (s.end()-s.start()+1);
819 qdp_vaxmy_chp(zptr, aptr, xptr, yptr, n_4vec);
820 }
821 else {
822 const int* tab = s.siteTable().slice();
823 for(int j = 0; j < s.numSiteTable() ; j++) {
824 int i = tab[j];
825
826 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
827 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
828 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
829
830 qdp_vaxmy_chp(zptr, aptr, xptr, yptr, 1);
831 }
832 }
833
834}
835
836// d = ax - P- y
837template<>
838inline
839void evaluate( OLattice< TVec > &d,
840 const OpAssign &op,
841 const QDPExpr<
842 BinaryNode<OpSubtract,
843 BinaryNode<OpMultiply,
844 Reference< QDPType< TScal, OScalar< TScal > > >,
845 Reference< QDPType< TVec, OLattice< TVec > > >
846 >,
847 UnaryNode<FnChiralProjectMinus, Reference< QDPType<TVec, OLattice<TVec> > > >
848 >,
849 OLattice< TVec > > &rhs,
850 const Subset& s)
851{
852#ifdef DEBUG_BLAS_G5
853 QDPIO::cout << "z = a*x + P{-} y" << endl;
854#endif
855
856
857 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right().child());
858
859 // ax is the left side of rhs and is in a binary node
860 typedef BinaryNode<OpMultiply,
861 Reference< QDPType< TScal, OScalar< TScal > > >,
862 Reference< QDPType< TVec, OLattice< TVec > > > > BN;
863
864 // get the binary node
865 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
866
867 // get a and x out of the bynary node
868 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
869 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
870 // Set pointers
871 REAL ar = a.elem().elem().elem().elem();
872 REAL *aptr = (REAL *)&ar;
873
874 if( s.hasOrderedRep() ) {
875 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
876 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
877 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
878
879
880 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
881 int n_4vec = (s.end()-s.start()+1);
882 qdp_vaxmy_chm(zptr, aptr, xptr, yptr, n_4vec);
883 }
884 else {
885 const int* tab = s.siteTable().slice();
886 for(int j = 0; j < s.numSiteTable() ; j++) {
887 int i = tab[j];
888
889 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
890 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
891 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
892
893 qdp_vaxmy_chm(zptr, aptr, xptr, yptr, 1);
894 }
895 }
896
897}
898
899// Vec = Scal * P_{+} Vec
900template<>
901inline
902void evaluate( OLattice< TVec > &d,
903 const OpAssign &op,
904 const QDPExpr<
905 BinaryNode<OpMultiply,
906 Reference< QDPType< TScal, OScalar< TScal > > >,
907 UnaryNode< FnChiralProjectPlus, Reference< QDPType< TVec, OLattice< TVec > > > >
908 >,
909 OLattice< TVec > > &rhs,
910 const Subset& s)
911{
912
913#ifdef DEBUG_BLAS_G5
914 cout << "BJ: v = a*P{+}v " << endl;
915#endif
916
917 const OLattice< TVec > &x = static_cast<const OLattice< TVec >&>(rhs.expression().right().child());
918 const OScalar< TScal > &a = static_cast<const OScalar< TScal >&>(rhs.expression().left());
919
920 REAL ar = a.elem().elem().elem().elem();
921 REAL *aptr = &ar;
922
923 if( s.hasOrderedRep() ) {
924 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
925 REAL *zptr = &(d.elem(s.start()).elem(0).elem(0).real());
926 int n_4vec = (s.end()-s.start()+1);
927
928 qdp_vscal_chp(zptr, aptr, xptr, n_4vec);
929 }
930 else {
931 const int* tab = s.siteTable().slice();
932 for(int j = 0; j < s.numSiteTable() ; j++) {
933 int i = tab[j];
934
935 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
936 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
937
938 qdp_vscal_chp(zptr, aptr, xptr, 1);
939 }
940 }
941
942}
943
944// Vec = Scal * P_{-} Vec
945template<>
946inline
947void evaluate( OLattice< TVec > &d,
948 const OpAssign &op,
949 const QDPExpr<
950 BinaryNode<OpMultiply,
951 Reference< QDPType< TScal, OScalar< TScal > > >,
952 UnaryNode< FnChiralProjectMinus, Reference< QDPType< TVec, OLattice< TVec > > > >
953 >,
954 OLattice< TVec > > &rhs,
955 const Subset& s)
956{
957
958#ifdef DEBUG_BLAS_G5
959 cout << "BJ: v = a*P{-}v " << endl;
960#endif
961
962 const OLattice< TVec > &x = static_cast<const OLattice< TVec >&>(rhs.expression().right().child());
963 const OScalar< TScal > &a = static_cast<const OScalar< TScal >&>(rhs.expression().left());
964
965 REAL ar = a.elem().elem().elem().elem();
966 REAL *aptr = &ar;
967
968 if( s.hasOrderedRep() ) {
969 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
970 REAL *zptr = &(d.elem(s.start()).elem(0).elem(0).real());
971 int n_4vec = (s.end()-s.start()+1);
972
973 qdp_vscal_chm(zptr, aptr, xptr, n_4vec);
974 }
975 else {
976 const int* tab = s.siteTable().slice();
977 for(int j = 0; j < s.numSiteTable() ; j++) {
978 int i = tab[j];
979
980 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
981 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
982
983 qdp_vscal_chm(zptr, aptr, xptr, 1);
984 }
985 }
986
987}
988
989// z = ax + bP+ y
990template<>
991inline
992void evaluate( OLattice< TVec > &d,
993 const OpAssign &op,
994 const QDPExpr<
995 BinaryNode<OpAdd,
996 BinaryNode<OpMultiply,
997 Reference< QDPType< TScal, OScalar< TScal > > >,
998 Reference< QDPType< TVec, OLattice< TVec > > > >,
999 BinaryNode<OpMultiply,
1000 Reference< QDPType< TScal, OScalar< TScal > > >,
1001 UnaryNode< FnChiralProjectPlus, Reference< QDPType< TVec, OLattice< TVec > > > >
1002 >
1003 >,
1004 OLattice< TVec > > &rhs,
1005 const Subset& s)
1006{
1007
1008#ifdef DEBUG_BLAS_G5
1009 QDPIO::cout << "z = a*x + b*P+y" << endl;
1010#endif
1011
1012 // Peel the stuff out of the expression
1013 // y is the right side of rhs
1014
1015 // ax is the left side of rhs and is in a binary node
1016 typedef BinaryNode<OpMultiply,
1017 Reference< QDPType< TScal, OScalar< TScal > > >,
1018 Reference< QDPType< TVec, OLattice< TVec > > > > BN1;
1019
1020 typedef BinaryNode<OpMultiply,
1021 Reference< QDPType< TScal, OScalar< TScal > > >,
1022 UnaryNode< FnChiralProjectPlus, Reference< QDPType< TVec, OLattice< TVec > > > > > BN2;
1023
1024 // get the binary node
1025 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
1026 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
1027
1028 // get a and x out of the binary node
1029 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
1030 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
1031
1032 // get b and y out of the binary node
1033 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1034 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right().child());
1035
1036
1037 // Set pointers
1038 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1039 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1040
1041 if( s.hasOrderedRep() ) {
1042 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1043 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1044 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1045
1046
1047 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1048 int n_4vec = (s.end()-s.start()+1);
1049 qdp_vaxpby_chp(zptr, aptr, xptr, bptr, yptr, n_4vec);
1050 }
1051 else {
1052 const int* tab = s.siteTable().slice();
1053 for(int j = 0; j < s.numSiteTable() ; j++) {
1054 int i = tab[j];
1055 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1056 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1057 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1058
1059 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1060 qdp_vaxpby_chp(zptr, aptr, xptr, bptr, yptr, 1);
1061
1062 }
1063 }
1064
1065}
1066
1067// z = ax + bP- y
1068template<>
1069inline
1070void evaluate( OLattice< TVec > &d,
1071 const OpAssign &op,
1072 const QDPExpr<
1073 BinaryNode<OpAdd,
1074 BinaryNode<OpMultiply,
1075 Reference< QDPType< TScal, OScalar< TScal > > >,
1076 Reference< QDPType< TVec, OLattice< TVec > > > >,
1077 BinaryNode<OpMultiply,
1078 Reference< QDPType< TScal, OScalar< TScal > > >,
1079 UnaryNode< FnChiralProjectMinus, Reference< QDPType< TVec, OLattice< TVec > > > >
1080 >
1081 >,
1082 OLattice< TVec > > &rhs,
1083 const Subset& s)
1084{
1085
1086#ifdef DEBUG_BLAS_G5
1087 QDPIO::cout << "z = a*x + b*P-y" << endl;
1088#endif
1089
1090 // Peel the stuff out of the expression
1091 // y is the right side of rhs
1092
1093 // ax is the left side of rhs and is in a binary node
1094 typedef BinaryNode<OpMultiply,
1095 Reference< QDPType< TScal, OScalar< TScal > > >,
1096 Reference< QDPType< TVec, OLattice< TVec > > > > BN1;
1097
1098 typedef BinaryNode<OpMultiply,
1099 Reference< QDPType< TScal, OScalar< TScal > > >,
1100 UnaryNode< FnChiralProjectMinus, Reference< QDPType< TVec, OLattice< TVec > > > > > BN2;
1101
1102 // get the binary node
1103 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
1104 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
1105
1106 // get a and x out of the binary node
1107 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
1108 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
1109
1110 // get b and y out of the binary node
1111 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1112 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right().child());
1113
1114
1115 // Set pointers
1116 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1117 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1118
1119 if( s.hasOrderedRep() ) {
1120 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1121 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1122 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1123
1124
1125 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1126 int n_4vec = (s.end()-s.start()+1);
1127 qdp_vaxpby_chm(zptr, aptr, xptr, bptr, yptr, n_4vec);
1128 }
1129 else {
1130 const int* tab = s.siteTable().slice();
1131 for(int j = 0; j < s.numSiteTable() ; j++) {
1132 int i = tab[j];
1133 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1134 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1135 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1136
1137 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1138 qdp_vaxpby_chm(zptr, aptr, xptr, bptr, yptr, 1);
1139
1140 }
1141 }
1142
1143}
1144
1145// z = ax - bP+ y
1146template<>
1147inline
1148void evaluate( OLattice< TVec > &d,
1149 const OpAssign &op,
1150 const QDPExpr<
1151 BinaryNode<OpSubtract,
1152 BinaryNode<OpMultiply,
1153 Reference< QDPType< TScal, OScalar< TScal > > >,
1154 Reference< QDPType< TVec, OLattice< TVec > > > >,
1155 BinaryNode<OpMultiply,
1156 Reference< QDPType< TScal, OScalar< TScal > > >,
1157 UnaryNode< FnChiralProjectPlus, Reference< QDPType< TVec, OLattice< TVec > > > >
1158 >
1159 >,
1160 OLattice< TVec > > &rhs,
1161 const Subset& s)
1162{
1163
1164#ifdef DEBUG_BLAS_G5
1165 QDPIO::cout << "z = a*x - b*P+y" << endl;
1166#endif
1167
1168 // Peel the stuff out of the expression
1169 // y is the right side of rhs
1170
1171 // ax is the left side of rhs and is in a binary node
1172 typedef BinaryNode<OpMultiply,
1173 Reference< QDPType< TScal, OScalar< TScal > > >,
1174 Reference< QDPType< TVec, OLattice< TVec > > > > BN1;
1175
1176 typedef BinaryNode<OpMultiply,
1177 Reference< QDPType< TScal, OScalar< TScal > > >,
1178 UnaryNode< FnChiralProjectPlus, Reference< QDPType< TVec, OLattice< TVec > > > > > BN2;
1179
1180 // get the binary node
1181 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
1182 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
1183
1184 // get a and x out of the binary node
1185 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
1186 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
1187
1188 // get b and y out of the binary node
1189 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1190 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right().child());
1191
1192
1193 // Set pointers
1194 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1195 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1196
1197 if( s.hasOrderedRep() ) {
1198 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1199 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1200 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1201
1202
1203 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1204 int n_4vec = (s.end()-s.start()+1);
1205 qdp_vaxmby_chp(zptr, aptr, xptr, bptr, yptr, n_4vec);
1206 }
1207 else {
1208 const int* tab = s.siteTable().slice();
1209 for(int j = 0; j < s.numSiteTable() ; j++) {
1210 int i = tab[j];
1211
1212 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1213 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1214 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1215
1216 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1217 qdp_vaxmby_chp(zptr, aptr, xptr, bptr, yptr, 1);
1218
1219 }
1220 }
1221
1222}
1223
1224// z = ax - bP- y
1225template<>
1226inline
1227void evaluate( OLattice< TVec > &d,
1228 const OpAssign &op,
1229 const QDPExpr<
1230 BinaryNode<OpSubtract,
1231 BinaryNode<OpMultiply,
1232 Reference< QDPType< TScal, OScalar< TScal > > >,
1233 Reference< QDPType< TVec, OLattice< TVec > > > >,
1234 BinaryNode<OpMultiply,
1235 Reference< QDPType< TScal, OScalar< TScal > > >,
1236 UnaryNode< FnChiralProjectMinus, Reference< QDPType< TVec, OLattice< TVec > > > >
1237 >
1238 >,
1239 OLattice< TVec > > &rhs,
1240 const Subset& s)
1241{
1242
1243#ifdef DEBUG_BLAS_G5
1244 QDPIO::cout << "z = a*x - b*P-y" << endl;
1245#endif
1246
1247 // Peel the stuff out of the expression
1248 // y is the right side of rhs
1249
1250 // ax is the left side of rhs and is in a binary node
1251 typedef BinaryNode<OpMultiply,
1252 Reference< QDPType< TScal, OScalar< TScal > > >,
1253 Reference< QDPType< TVec, OLattice< TVec > > > > BN1;
1254
1255 typedef BinaryNode<OpMultiply,
1256 Reference< QDPType< TScal, OScalar< TScal > > >,
1257 UnaryNode< FnChiralProjectMinus, Reference< QDPType< TVec, OLattice< TVec > > > > > BN2;
1258
1259 // get the binary node
1260 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
1261 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
1262
1263 // get a and x out of the binary node
1264 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
1265 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
1266
1267 // get b and y out of the binary node
1268 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1269 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right().child());
1270
1271
1272 // Set pointers
1273 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1274 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1275
1276 if( s.hasOrderedRep() ) {
1277 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1278 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1279 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1280
1281
1282 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1283 int n_4vec = (s.end()-s.start()+1);
1284 qdp_vaxmby_chm(zptr, aptr, xptr, bptr, yptr, n_4vec);
1285
1286 }
1287 else {
1288 const int* tab = s.siteTable().slice();
1289 for(int j = 0; j < s.numSiteTable() ; j++) {
1290 int i = tab[j];
1291
1292 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1293 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1294 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1295
1296 qdp_vaxmby_chm(zptr, aptr, xptr, bptr, yptr, 1);
1297 }
1298 }
1299
1300}
1301
1302// Vec = Scal * GammaConst<Ns,Ns-1>* Vec
1303template<>
1304inline
1305void evaluate( OLattice< TVec > &d,
1306 const OpAssign &op,
1307 const QDPExpr<
1308 BinaryNode<OpMultiply,
1309 Reference< QDPType< TScal, OScalar< TScal > > >,
1310 BinaryNode<
1311 OpGammaConstMultiply,
1312 GammaConst<Ns,Ns*Ns-1>,
1313 Reference< QDPType< TVec, OLattice< TVec > > >
1314 >
1315 >,
1316 OLattice< TVec >
1317 > &rhs,
1318 const Subset& s)
1319{
1320
1321#ifdef DEBUG_BLAS_G5
1322 QDPIO::cout << "z = a*(GammaConst<Ns,Ns*Ns-1>()*x)" << endl;
1323#endif
1324
1325
1326 typedef BinaryNode<
1327 OpGammaConstMultiply,
1328 GammaConst<Ns,Ns*Ns-1>,
1329 Reference< QDPType< TVec, OLattice< TVec > > >
1330 > BN1;
1331
1332 const OScalar< TScal > &a = static_cast<const OScalar< TScal >&>(rhs.expression().left());
1333 const BN1 &node = static_cast<const BN1&>(rhs.expression().right());
1334
1335 const OLattice< TVec > &x = static_cast<const OLattice< TVec >&>(node.right());
1336
1337 REAL ar = a.elem().elem().elem().elem();
1338 REAL *aptr = &ar;
1339
1340 if( s.hasOrderedRep() ) {
1341 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1342 REAL *zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1343 int n_4vec = (s.end()-s.start()+1);
1344
1345 qdp_vscal_g5(zptr, aptr, xptr, n_4vec);
1346 }
1347 else {
1348 const int* tab = s.siteTable().slice();
1349 for(int j = 0; j < s.numSiteTable() ; j++) {
1350 int i = tab[j];
1351
1352 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1353 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1354
1355 qdp_vscal_g5(zptr, aptr, xptr, 1);
1356 }
1357 }
1358
1359}
1360
1361
1362// Vec = Vec - a*Gamma5*Vec
1363//
1364template<>
1365inline
1366void evaluate( OLattice< TVec > &d,
1367 const OpAssign &op,
1368 const QDPExpr<
1369 BinaryNode<OpSubtract,
1370 Reference< QDPType< TVec, OLattice< TVec > > >,
1371 BinaryNode<OpMultiply,
1372 Reference< QDPType< TScal, OScalar< TScal > > >,
1373 BinaryNode<
1374 OpGammaConstMultiply,
1375 GammaConst<Ns,Ns*Ns-1>,
1376 Reference< QDPType< TVec, OLattice< TVec > > >
1377 >
1378 >
1379 >,
1380 OLattice< TVec > > &rhs,
1381 const Subset& s)
1382{
1383
1384#ifdef DEBUG_BLAS_G5
1385 QDPIO::cout << "z = x - a GammaConst<Ns,Ns*Ns-1>()*y" << endl;
1386#endif
1387 const OLattice<TVec>& x = static_cast<const OLattice<TVec>&>(rhs.expression().left());
1388
1389 typedef BinaryNode<OpMultiply,
1390 Reference< QDPType< TScal, OScalar< TScal > > >,
1391 BinaryNode<
1392 OpGammaConstMultiply,
1393 GammaConst<Ns,Ns*Ns-1>,
1394 Reference< QDPType< TVec, OLattice< TVec > > >
1395 >
1396 > MN;
1397 const MN& mul_node = static_cast<const MN&>(rhs.expression().right());
1398
1399 const OScalar<TScal>& a = static_cast<const OScalar<TScal>& >(mul_node.left());
1400
1401 typedef BinaryNode<
1402 OpGammaConstMultiply,
1403 GammaConst<Ns,Ns*Ns-1>,
1404 Reference< QDPType< TVec, OLattice< TVec > > >
1405 > GN;
1406
1407 const GN& gamma_node = static_cast<const GN&>(mul_node.right());
1408 const OLattice<TVec>& y = static_cast<const OLattice<TVec>&>(gamma_node.right());
1409
1410 REAL ar = a.elem().elem().elem().elem();
1411 REAL *aptr = (REAL *)&ar;
1412
1413 if( s.hasOrderedRep() ) {
1414 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1415 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1416 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1417
1418
1419 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1420 int n_4vec = (s.end()-s.start()+1);
1421 qdp_vxmag5y(zptr, aptr, xptr, yptr, n_4vec);
1422 }
1423 else {
1424 const int* tab = s.siteTable().slice();
1425 for(int j = 0; j < s.numSiteTable() ; j++) {
1426 int i = tab[j];
1427
1428 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1429 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1430 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1431 qdp_vxmag5y(zptr, aptr, xptr, yptr, 1);
1432 }
1433 }
1434
1435}
1436
1437
1438// Vec = a*Vec + b*Gamma5*Vec
1439//
1440template<>
1441inline
1442void evaluate( OLattice< TVec > &d,
1443 const OpAssign &op,
1444 const QDPExpr<
1445 BinaryNode<OpAdd,
1446 BinaryNode<OpMultiply,
1447 Reference< QDPType<TScal, OScalar< TScal > > >,
1448 Reference< QDPType<TVec, OLattice< TVec > > >
1449 >,
1450 BinaryNode<OpMultiply,
1451 Reference< QDPType< TScal, OScalar< TScal > > >,
1452 BinaryNode<
1453 OpGammaConstMultiply,
1454 GammaConst<Ns,Ns*Ns-1>,
1455 Reference< QDPType< TVec, OLattice< TVec > > >
1456 >
1457 >
1458 >,
1459 OLattice< TVec > > &rhs,
1460 const Subset& s)
1461{
1462
1463#ifdef DEBUG_BLAS_G5
1464 QDPIO::cout << "z = a*x + b*GammaConst<Ns,Ns*Ns-1>()*y" << endl;
1465#endif
1466
1467 typedef BinaryNode<OpMultiply,
1468 Reference< QDPType<TScal, OScalar< TScal > > >,
1469 Reference< QDPType<TVec, OLattice< TVec > > >
1470 > MN1;
1471
1472 typedef BinaryNode<OpMultiply,
1473 Reference< QDPType< TScal, OScalar< TScal > > >,
1474 BinaryNode<
1475 OpGammaConstMultiply,
1476 GammaConst<Ns,Ns*Ns-1>,
1477 Reference< QDPType< TVec, OLattice< TVec > > >
1478 >
1479 > MN2;
1480
1481 typedef BinaryNode<
1482 OpGammaConstMultiply,
1483 GammaConst<Ns,Ns*Ns-1>,
1484 Reference< QDPType< TVec, OLattice< TVec > > >
1485 > GN;
1486
1487 const MN1& mulNode1 = static_cast< const MN1& >(rhs.expression().left());
1488 const MN2& mulNode2 = static_cast< const MN2& >(rhs.expression().right());
1489 const GN& gammaNode = static_cast< const GN& >(mulNode2.right());
1490
1491 const OScalar<TScal>& a = static_cast<const OScalar<TScal>&>(mulNode1.left());
1492 const OLattice<TVec>& x = static_cast<const OLattice<TVec>&>(mulNode1.right());
1493
1494 const OScalar<TScal>& b = static_cast<const OScalar<TScal>&>(mulNode2.left());
1495 const OLattice<TVec>& y = static_cast<const OLattice<TVec>&>(gammaNode.right());
1496
1497
1498 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1499 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1500
1501 if( s.hasOrderedRep() ) {
1502 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1503 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1504 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1505
1506
1507 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1508 int n_4vec = (s.end()-s.start()+1);
1509 qdp_vaxpbg5y(zptr, aptr, xptr, bptr, yptr, n_4vec);
1510 }
1511 else {
1512 const int* tab = s.siteTable().slice();
1513 for(int j = 0; j < s.numSiteTable() ; j++) {
1514 int i = tab[j];
1515
1516 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1517 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1518 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1519 qdp_vaxpbg5y(zptr, aptr, xptr, bptr, yptr, 1);
1520
1521 }
1522 }
1523
1524}
1525
1526// Vec = Gamma_5 *( a*Vec - b*Vec )
1527//
1528template<>
1529inline
1530void evaluate( OLattice< TVec > &d,
1531 const OpAssign &op,
1532 const QDPExpr<
1533 BinaryNode<OpGammaConstMultiply,
1534 GammaConst<Ns,Ns*Ns-1>,
1535 BinaryNode<OpSubtract,
1536 BinaryNode<OpMultiply,
1537 Reference< QDPType<TScal, OScalar< TScal > > >,
1538 Reference< QDPType<TVec, OLattice< TVec > > >
1539 >,
1540 BinaryNode<OpMultiply,
1541 Reference< QDPType< TScal, OScalar< TScal > > >,
1542 Reference< QDPType< TVec, OLattice< TVec > > >
1543 >
1544 >
1545 >,
1546 OLattice< TVec > > &rhs,
1547 const Subset& s)
1548{
1549
1550#ifdef DEBUG_BLAS_G5
1551 QDPIO::cout << "z = GammaConst<Ns,Ns*Ns-1>()*(ax - by)" << endl;
1552#endif
1553
1554 typedef BinaryNode<OpSubtract,
1555 BinaryNode<OpMultiply,
1556 Reference< QDPType<TScal, OScalar< TScal > > >,
1557 Reference< QDPType<TVec, OLattice< TVec > > >
1558 >,
1559 BinaryNode<OpMultiply,
1560 Reference< QDPType< TScal, OScalar< TScal > > >,
1561 Reference< QDPType< TVec, OLattice< TVec > > >
1562 >
1563 > AXMBY;
1564
1565 const AXMBY& axmby_node = static_cast<const AXMBY&>(rhs.expression().right());
1566
1567
1568 typedef BinaryNode<OpMultiply,
1569 Reference< QDPType<TScal, OScalar< TScal > > >,
1570 Reference< QDPType<TVec, OLattice< TVec > > >
1571 > MN;
1572
1573 const MN& mulNode1 = static_cast<const MN& >( axmby_node.left());
1574 const MN& mulNode2 = static_cast<const MN& >( axmby_node.right());
1575
1576 const OScalar<TScal>& a = static_cast<const OScalar<TScal>&>(mulNode1.left());
1577 const OLattice<TVec>& x = static_cast<const OLattice<TVec>&>(mulNode1.right());
1578
1579 const OScalar<TScal>& b = static_cast<const OScalar<TScal>&>(mulNode2.left());
1580 const OLattice<TVec>& y = static_cast<const OLattice<TVec>&>(mulNode2.right());
1581
1582
1583 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1584 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1585
1586 if( s.hasOrderedRep() ) {
1587 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1588 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1589 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1590
1591
1592 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1593 int n_4vec = (s.end()-s.start()+1);
1594 qdp_vg5axmby(zptr, aptr, xptr, bptr, yptr, n_4vec);
1595 }
1596 else {
1597 const int* tab = s.siteTable().slice();
1598 for(int j = 0; j < s.numSiteTable() ; j++) {
1599 int i = tab[j];
1600
1601 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1602 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1603 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1604
1605 qdp_vg5axmby(zptr, aptr, xptr, bptr, yptr, 1);
1606 }
1607 }
1608
1609}
1610
1611// --------------- These here come from generics ----------------------
1612// Vec = a*Vec + b*Gamma5*i*Vec
1613//
1614template<>
1615inline
1616void evaluate( OLattice< TVec > &d,
1617 const OpAssign &op,
1618 const QDPExpr<
1619 BinaryNode<OpAdd,
1620 BinaryNode<OpMultiply,
1621 Reference< QDPType<TScal, OScalar< TScal > > >,
1622 Reference< QDPType<TVec, OLattice< TVec > > >
1623 >,
1624 BinaryNode<OpMultiply,
1625 Reference< QDPType< TScal, OScalar< TScal > > >,
1626 BinaryNode<
1627 OpGammaConstMultiply,
1628 GammaConst<Ns,Ns*Ns-1>,
1629 UnaryNode<FnTimesI,
1630 Reference< QDPType<TVec, OLattice< TVec > > >
1631 >
1632 >
1633 >
1634 >,
1635 OLattice< TVec > > &rhs,
1636 const Subset& s)
1637{
1638
1639#ifdef DEBUG_BLAS_G5
1640 QDPIO::cout << "z = a*x + b*GammaConst<Ns,Ns*Ns-1>()*timesI(y)" << endl;
1641#endif
1642
1643 typedef BinaryNode<OpMultiply,
1644 Reference< QDPType<TScal, OScalar< TScal > > >,
1645 Reference< QDPType<TVec, OLattice< TVec > > >
1646 > MN1;
1647
1648 typedef BinaryNode<OpMultiply,
1649 Reference< QDPType< TScal, OScalar< TScal > > >,
1650 BinaryNode<
1651 OpGammaConstMultiply,
1652 GammaConst<Ns,Ns*Ns-1>,
1653 UnaryNode< FnTimesI,
1654 Reference< QDPType< TVec, OLattice< TVec > > >
1655 >
1656 >
1657 > MN2;
1658
1659 typedef BinaryNode<
1660 OpGammaConstMultiply,
1661 GammaConst<Ns,Ns*Ns-1>,
1662 UnaryNode< FnTimesI,
1663 Reference< QDPType< TVec, OLattice< TVec > > >
1664 >
1665 > GN;
1666
1667 typedef UnaryNode< FnTimesI,
1668 Reference< QDPType< TVec, OLattice< TVec > > >
1669 > IN;
1670
1671 const MN1& mulNode1 = static_cast< const MN1& >(rhs.expression().left());
1672 const MN2& mulNode2 = static_cast< const MN2& >(rhs.expression().right());
1673 const GN& gammaNode = static_cast< const GN& >(mulNode2.right());
1674 const IN& mulINode = static_cast< const IN& >(gammaNode.right());
1675
1676 const OScalar<TScal>& a = static_cast<const OScalar<TScal>&>(mulNode1.left());
1677 const OLattice<TVec>& x = static_cast<const OLattice<TVec>&>(mulNode1.right());
1678
1679 const OScalar<TScal>& b = static_cast<const OScalar<TScal>&>(mulNode2.left());
1680 const OLattice<TVec>& y = static_cast<const OLattice<TVec>&>(mulINode.child());
1681
1682
1683 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1684 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1685
1686 if( s.hasOrderedRep()) {
1687 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1688 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1689 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1690
1691
1692 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1693 int n_4vec = (s.end()-s.start()+1);
1694 qdp_vaxpbg5iy(zptr, aptr, xptr, bptr, yptr, n_4vec);
1695 }
1696 else {
1697 const int* tab = s.siteTable().slice();
1698 for(int j = 0; j < s.numSiteTable() ; j++) {
1699 int i = tab[j];
1700
1701 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1702 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1703 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1704
1705 qdp_vaxpbg5iy(zptr, aptr, xptr, bptr, yptr, 1);
1706 }
1707 }
1708
1709}
1710
1711// Vec = a*Vec - b*Gamma5*i*Vec
1712//
1713template<>
1714inline
1715void evaluate( OLattice< TVec > &d,
1716 const OpAssign &op,
1717 const QDPExpr<
1718 BinaryNode<OpSubtract,
1719 BinaryNode<OpMultiply,
1720 Reference< QDPType<TScal, OScalar< TScal > > >,
1721 Reference< QDPType<TVec, OLattice< TVec > > >
1722 >,
1723 BinaryNode<OpMultiply,
1724 Reference< QDPType< TScal, OScalar< TScal > > >,
1725 BinaryNode<
1726 OpGammaConstMultiply,
1727 GammaConst<Ns,Ns*Ns-1>,
1728 UnaryNode<FnTimesI,
1729 Reference< QDPType<TVec, OLattice< TVec > > >
1730 >
1731 >
1732 >
1733 >,
1734 OLattice< TVec > > &rhs,
1735 const Subset& s)
1736{
1737
1738#ifdef DEBUG_BLAS_G5
1739 QDPIO::cout << "z = a*x + b*GammaConst<Ns,Ns*Ns-1>()*timesI(y)" << endl;
1740#endif
1741
1742 typedef BinaryNode<OpMultiply,
1743 Reference< QDPType<TScal, OScalar< TScal > > >,
1744 Reference< QDPType<TVec, OLattice< TVec > > >
1745 > MN1;
1746
1747 typedef BinaryNode<OpMultiply,
1748 Reference< QDPType< TScal, OScalar< TScal > > >,
1749 BinaryNode<
1750 OpGammaConstMultiply,
1751 GammaConst<Ns,Ns*Ns-1>,
1752 UnaryNode< FnTimesI,
1753 Reference< QDPType< TVec, OLattice< TVec > > >
1754 >
1755 >
1756 > MN2;
1757
1758 typedef BinaryNode<
1759 OpGammaConstMultiply,
1760 GammaConst<Ns,Ns*Ns-1>,
1761 UnaryNode< FnTimesI,
1762 Reference< QDPType< TVec, OLattice< TVec > > >
1763 >
1764 > GN;
1765
1766 typedef UnaryNode< FnTimesI,
1767 Reference< QDPType< TVec, OLattice< TVec > > >
1768 > IN;
1769
1770 const MN1& mulNode1 = static_cast< const MN1& >(rhs.expression().left());
1771 const MN2& mulNode2 = static_cast< const MN2& >(rhs.expression().right());
1772 const GN& gammaNode = static_cast< const GN& >(mulNode2.right());
1773 const IN& mulINode = static_cast< const IN& >(gammaNode.right());
1774
1775 const OScalar<TScal>& a = static_cast<const OScalar<TScal>&>(mulNode1.left());
1776 const OLattice<TVec>& x = static_cast<const OLattice<TVec>&>(mulNode1.right());
1777
1778 const OScalar<TScal>& b = static_cast<const OScalar<TScal>&>(mulNode2.left());
1779 const OLattice<TVec>& y = static_cast<const OLattice<TVec>&>(mulINode.child());
1780
1781
1782 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1783 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1784
1785 if( s.hasOrderedRep() ) {
1786 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1787 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1788 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1789
1790 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1791 int n_4vec = (s.end()-s.start()+1);
1792 qdp_vaxmbg5iy(zptr, aptr, xptr, bptr, yptr, n_4vec);
1793 }
1794 else {
1795 const int* tab = s.siteTable().slice();
1796 for(int j = 0; j < s.numSiteTable() ; j++) {
1797 int i = tab[j];
1798
1799 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1800 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1801 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1802
1803 qdp_vaxmbg5iy(zptr, aptr, xptr, bptr, yptr, 1);
1804 }
1805 }
1806
1807}
1808
1809// Vec = Vec + a*Gamma5*i*Vec
1810//
1811template<>
1812inline
1813void evaluate( OLattice< TVec > &d,
1814 const OpAssign &op,
1815 const QDPExpr<
1816 BinaryNode<OpAdd,
1817 Reference< QDPType< TVec, OLattice< TVec > > >,
1818 BinaryNode<OpMultiply,
1819 Reference< QDPType< TScal, OScalar< TScal > > >,
1820 BinaryNode<
1821 OpGammaConstMultiply,
1822 GammaConst<Ns,Ns*Ns-1>,
1823 UnaryNode<FnTimesI,
1824 Reference< QDPType< TVec, OLattice< TVec > > >
1825 >
1826 >
1827 >
1828 >,
1829 OLattice< TVec > > &rhs,
1830 const Subset& s)
1831{
1832
1833#ifdef DEBUG_BLAS_G5
1834 QDPIO::cout << "z = x + a GammaConst<Ns,Ns*Ns-1>()*i*y" << endl;
1835#endif
1836 const OLattice<TVec>& x = static_cast<const OLattice<TVec>&>(rhs.expression().left());
1837
1838 typedef BinaryNode<OpMultiply,
1839 Reference< QDPType< TScal, OScalar< TScal > > >,
1840 BinaryNode<
1841 OpGammaConstMultiply,
1842 GammaConst<Ns,Ns*Ns-1>,
1843 UnaryNode<FnTimesI,
1844 Reference< QDPType< TVec, OLattice< TVec > > >
1845 >
1846 >
1847 > MN;
1848 const MN& mul_node = static_cast<const MN&>(rhs.expression().right());
1849
1850 const OScalar<TScal>& a = static_cast<const OScalar<TScal>& >(mul_node.left());
1851
1852 typedef BinaryNode<
1853 OpGammaConstMultiply,
1854 GammaConst<Ns,Ns*Ns-1>,
1855 UnaryNode<FnTimesI,
1856 Reference< QDPType< TVec, OLattice< TVec > > >
1857 >
1858 > GN;
1859
1860 typedef UnaryNode<FnTimesI,
1861 Reference< QDPType< TVec, OLattice< TVec > > >
1862 > IN;
1863
1864 const GN& gamma_node = static_cast<const GN&>(mul_node.right());
1865 const IN& timesI_node = static_cast<const IN&>(gamma_node.right());
1866
1867 const OLattice<TVec>& y = static_cast<const OLattice<TVec>&>(timesI_node.child());
1868
1869 REAL ar = a.elem().elem().elem().elem();
1870 REAL *aptr = (REAL *)&ar;
1871 if( s.hasOrderedRep() ) {
1872 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1873 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1874 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1875
1876
1877 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1878 int n_4vec = (s.end()-s.start()+1);
1879 qdp_vxpag5iy(zptr, aptr, xptr, yptr, n_4vec);
1880 }
1881 else {
1882 const int* tab = s.siteTable().slice();
1883 for(int j = 0; j < s.numSiteTable() ; j++) {
1884 int i = tab[j];
1885
1886 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1887 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1888 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1889
1890 qdp_vxpag5iy(zptr, aptr, xptr, yptr, 1);
1891
1892
1893 }
1894 }
1895
1896}
1897
1898// Vec = Vec - a*Gamma5*i*Vec
1899//
1900template<>
1901inline
1902void evaluate( OLattice< TVec > &d,
1903 const OpAssign &op,
1904 const QDPExpr<
1905 BinaryNode<OpSubtract,
1906 Reference< QDPType< TVec, OLattice< TVec > > >,
1907 BinaryNode<OpMultiply,
1908 Reference< QDPType< TScal, OScalar< TScal > > >,
1909 BinaryNode<
1910 OpGammaConstMultiply,
1911 GammaConst<Ns,Ns*Ns-1>,
1912 UnaryNode<FnTimesI,
1913 Reference< QDPType< TVec, OLattice< TVec > > >
1914 >
1915 >
1916 >
1917 >,
1918 OLattice< TVec > > &rhs,
1919 const Subset& s)
1920{
1921
1922#ifdef DEBUG_BLAS_G5
1923 QDPIO::cout << "z = x - a GammaConst<Ns,Ns*Ns-1>()*i*y" << endl;
1924#endif
1925 const OLattice<TVec>& x = static_cast<const OLattice<TVec>&>(rhs.expression().left());
1926
1927 typedef BinaryNode<OpMultiply,
1928 Reference< QDPType< TScal, OScalar< TScal > > >,
1929 BinaryNode<
1930 OpGammaConstMultiply,
1931 GammaConst<Ns,Ns*Ns-1>,
1932 UnaryNode<FnTimesI,
1933 Reference< QDPType< TVec, OLattice< TVec > > >
1934 >
1935 >
1936 > MN;
1937 const MN& mul_node = static_cast<const MN&>(rhs.expression().right());
1938
1939 const OScalar<TScal>& a = static_cast<const OScalar<TScal>& >(mul_node.left());
1940
1941 typedef BinaryNode<
1942 OpGammaConstMultiply,
1943 GammaConst<Ns,Ns*Ns-1>,
1944 UnaryNode<FnTimesI,
1945 Reference< QDPType< TVec, OLattice< TVec > > >
1946 >
1947 > GN;
1948
1949 typedef UnaryNode<FnTimesI,
1950 Reference< QDPType< TVec, OLattice< TVec > > >
1951 > IN;
1952
1953 const GN& gamma_node = static_cast<const GN&>(mul_node.right());
1954 const IN& timesI_node = static_cast<const IN&>(gamma_node.right());
1955
1956 const OLattice<TVec>& y = static_cast<const OLattice<TVec>&>(timesI_node.child());
1957
1958 REAL ar = a.elem().elem().elem().elem();
1959 REAL *aptr = (REAL *)&ar;
1960
1961 if( s.hasOrderedRep() ) {
1962 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1963 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1964 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1965
1966
1967 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1968 int n_4vec = (s.end()-s.start()+1);
1969 qdp_vxmag5iy(zptr, aptr, xptr, yptr, n_4vec);
1970 }
1971 else {
1972 const int* tab = s.siteTable().slice();
1973 for(int j = 0; j < s.numSiteTable() ; j++) {
1974 int i = tab[j];
1975
1976 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1977 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1978 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1979
1980 qdp_vxmag5iy(zptr, aptr, xptr, yptr, 1);
1981 }
1982 }
1983
1984}
1985
1986
1987
1988// Vec += a*Gamma5*i*Vec
1989//
1990template<>
1991inline
1992void evaluate( OLattice< TVec > &d,
1993 const OpAddAssign &op,
1994 const QDPExpr<
1995 BinaryNode<OpMultiply,
1996 Reference< QDPType< TScal, OScalar< TScal > > >,
1997 BinaryNode<
1998 OpGammaConstMultiply,
1999 GammaConst<Ns,Ns*Ns-1>,
2000 UnaryNode<FnTimesI,
2001 Reference< QDPType< TVec, OLattice< TVec > > >
2002 >
2003 >
2004 >,
2005 OLattice< TVec >
2006 > &rhs,
2007 const Subset& s)
2008{
2009
2010#ifdef DEBUG_BLAS_G5
2011 QDPIO::cout << "z += a GammaConst<Ns,Ns*Ns-1>()*i*y" << endl;
2012#endif
2013
2014#if 0
2015 typedef BinaryNode<OpMultiply,
2016 Reference< QDPType< TScal, OScalar< TScal > > >,
2017 BinaryNode<
2018 OpGammaConstMultiply,
2019 GammaConst<Ns,Ns*Ns-1>,
2020 UnaryNode<FnTimesI,
2021 Reference< QDPType< TVec, OLattice< TVec > > >
2022 >
2023 >
2024 > MN;
2025
2026 const MN& mul_node = static_cast<const MN&>(rhs.expression().right());
2027#endif
2028
2029 const OScalar<TScal>& a = static_cast<const OScalar<TScal>& >(rhs.expression().left());
2030
2031 typedef BinaryNode<
2032 OpGammaConstMultiply,
2033 GammaConst<Ns,Ns*Ns-1>,
2034 UnaryNode<FnTimesI,
2035 Reference< QDPType< TVec, OLattice< TVec > > >
2036 >
2037 > GN;
2038
2039 typedef UnaryNode<FnTimesI,
2040 Reference< QDPType< TVec, OLattice< TVec > > >
2041 > IN;
2042
2043 const GN& gamma_node = static_cast<const GN&>(rhs.expression().right());
2044 const IN& timesI_node = static_cast<const IN&>(gamma_node.right());
2045
2046 const OLattice<TVec>& x = static_cast<const OLattice<TVec>&>(timesI_node.child());
2047
2048 REAL ar = a.elem().elem().elem().elem();
2049 REAL *aptr = (REAL *)&ar;
2050
2051 if( s.hasOrderedRep() ) {
2052 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
2053 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
2054
2055
2056 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2057 int n_4vec = (s.end()-s.start()+1);
2058 qdp_vxpag5iy(zptr, aptr, zptr, xptr, n_4vec);
2059 }
2060 else {
2061 const int* tab = s.siteTable().slice();
2062 for(int j = 0; j < s.numSiteTable() ; j++) {
2063 int i = tab[j];
2064
2065 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
2066 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
2067
2068 qdp_vxpag5iy(zptr, aptr, zptr, xptr, 1);
2069
2070 }
2071 }
2072
2073}
2074
2075
2076// Vec -= a*Gamma5*i*Vec
2077//
2078template<>
2079inline
2080void evaluate( OLattice< TVec > &d,
2081 const OpSubtractAssign &op,
2082 const QDPExpr<
2083 BinaryNode<OpMultiply,
2084 Reference< QDPType< TScal, OScalar< TScal > > >,
2085 BinaryNode<
2086 OpGammaConstMultiply,
2087 GammaConst<Ns,Ns*Ns-1>,
2088 UnaryNode<FnTimesI,
2089 Reference< QDPType< TVec, OLattice< TVec > > >
2090 >
2091 >
2092 >,
2093 OLattice< TVec >
2094 > &rhs,
2095 const Subset& s)
2096{
2097
2098#ifdef DEBUG_BLAS_G5
2099 QDPIO::cout << "z -= a GammaConst<Ns,Ns*Ns-1>()*i*y" << endl;
2100#endif
2101
2102#if 0
2103 typedef BinaryNode<OpMultiply,
2104 Reference< QDPType< TScal, OScalar< TScal > > >,
2105 BinaryNode<
2106 OpGammaConstMultiply,
2107 GammaConst<Ns,Ns*Ns-1>,
2108 UnaryNode<FnTimesI,
2109 Reference< QDPType< TVec, OLattice< TVec > > >
2110 >
2111 >
2112 > MN;
2113
2114 const MN& mul_node = static_cast<const MN&>(rhs.expression().right());
2115#endif
2116
2117 const OScalar<TScal>& a = static_cast<const OScalar<TScal>& >(rhs.expression().left());
2118
2119 typedef BinaryNode<
2120 OpGammaConstMultiply,
2121 GammaConst<Ns,Ns*Ns-1>,
2122 UnaryNode<FnTimesI,
2123 Reference< QDPType< TVec, OLattice< TVec > > >
2124 >
2125 > GN;
2126
2127 typedef UnaryNode<FnTimesI,
2128 Reference< QDPType< TVec, OLattice< TVec > > >
2129 > IN;
2130
2131 const GN& gamma_node = static_cast<const GN&>(rhs.expression().right());
2132 const IN& timesI_node = static_cast<const IN&>(gamma_node.right());
2133
2134 const OLattice<TVec>& x = static_cast<const OLattice<TVec>&>(timesI_node.child());
2135
2136 REAL ar = a.elem().elem().elem().elem();
2137 REAL *aptr = (REAL *)&ar;
2138
2139 if( s.hasOrderedRep() ) {
2140 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
2141 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
2142
2143
2144 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2145 int n_4vec = (s.end()-s.start()+1);
2146 qdp_vxmag5iy(zptr, aptr, zptr, xptr, n_4vec);
2147 }
2148 else {
2149 const int* tab = s.siteTable().slice();
2150 for(int j = 0; j < s.numSiteTable() ; j++) {
2151 int i = tab[j];
2152
2153 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
2154 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
2155
2156 qdp_vxmag5iy(zptr, aptr, zptr, xptr, 1);
2157
2158 }
2159 }
2160
2161}
2162
2163} // namespace QDP;
2164
2165#endif // guard
2166
QDP::OLattice
Outer grid Lattice type.
Definition qdp_outer.h:264
QDP::OLattice::elem
T & elem(int i)
Definition qdp_outer.h:400
QDP::OScalar
Outer grid Scalar class *‍/.
Definition qdp_outer.h:37
QDP::PScalar
Primitive Scalar.
Definition qdp_primscalar.h:25
QDP::PSpinVector
Primitive spin Vector class.
Definition qdp_primspinvec.h:33
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::Subset::end
int end() const
Definition qdp_subset.h:81
QDP::Subset::siteTable
const multi1d< int > & siteTable() const
Definition qdp_subset.h:83
QDP::Subset::numSiteTable
int numSiteTable() const
Definition qdp_subset.h:84
QDP::Subset::hasOrderedRep
bool hasOrderedRep() const
Definition qdp_subset.h:79
QDP::Subset::start
int start() const
Definition qdp_subset.h:80
QDP::multi1d::slice
const T * slice() const
Return ref to a column slice.
Definition qdp_multi.h:225
QDP::REAL
REAL32 REAL
Definition qdp_scalarsite_defs.h:28
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::Ns
const int Ns
Definition qdp_params.h:26
QDP::QDPIO::cout
StandardOutputStream cout
Definition qdp_stdio.cc:21
QDP
Yet another random number generator.
Definition qdp_bgq_threadbind.cc:16
QDP::TScal
PScalar< PScalar< RScalar< REAL > > > TScal
Definition qdp_scalarsite_bagel_qdp_blas.h:20
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::TVec
PSpinVector< PColorVector< RComplex< REAL >, 3 >, Ns > TVec
Definition qdp_scalarsite_bagel_qdp_blas.h:19
Generated by doxygen 1.16.1