QDP++
qdp_scalarsite_generic_blas.h
Go to the documentation of this file.
1// $Id: qdp_scalarsite_generic_blas.h,v 1.26 2009-09-15 20:48:42 bjoo Exp $
2
8
9
10#ifndef QDP_SCALARSITE_GENERIC_BLAS_H
11#define QDP_SCALARSITE_GENERIC_BLAS_H
12
13#include "scalarsite_generic/generic_blas_vaxpy3.h"
14#include "scalarsite_generic/generic_blas_vaxmy3.h"
15#include "scalarsite_generic/generic_blas_vaxpby3.h"
22#include "scalarsite_generic/generic_blas_local_vcdot_real.h"
23
24namespace QDP {
25
26// Types needed for the expression templates.
27// TVec has outer Ns template so it ought to work for staggered as well
28typedef PSpinVector<PColorVector<RComplex<REAL>, 3>, 4> TVec;
29typedef PScalar<PScalar<RScalar<REAL> > > TScal;
30
32// Threading evaluates
33//
34// by Xu Guo, EPCC, 12 August, 2008
36
37// the wrappers for the functions to be threaded
38#include "qdp_scalarsite_generic_blas_wrapper.h"
39
40
41// #define DEBUG_BLAS
42// TVec is the LatticeFermion from qdp_dwdefs.h with the OLattice<> stripped
43// from around it
44
45// TScalar is the usual Real, with the OScalar<> stripped from it
46//
47// THis is simply to make the code more readable, and reduces < < s and > >s
48// in the template arguments
49
50// d += Scalar*Vec
51template<>
52inline
53void evaluate(OLattice< TVec >& d,
54 const OpAddAssign& op,
55 const QDPExpr<BinaryNode<OpMultiply,
56 Reference< QDPType< TScal, OScalar < TScal > > >,
57 Reference< QDPType< TVec, OLattice< TVec > > > >,
58 OLattice< TVec > > &rhs,
59 const Subset& s)
60{
61
62#ifdef DEBUG_BLAS
63 QDPIO::cout << "y += a*x" << endl;
64#endif
65
66 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().right());
67 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().left());
68
69 REAL ar = a.elem().elem().elem().elem();
70 REAL* aptr = &ar;
71
72 if( s.hasOrderedRep() ) {
73 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
74 REAL* yptr = &(d.elem(s.start()).elem(0).elem(0).real());
75 // cout << "Specialised axpy a ="<< ar << endl;
76
77 int total_n_3vec = (s.end()-s.start()+1);
78
79 ordered_vaxpy3_user_arg a = {yptr, aptr, xptr, yptr};
80
81 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
82
84 // Original code
86 //int n_3vec = (s.end()-s.start()+1)*Ns;
87 //vaxpy3(yptr, aptr, xptr, yptr, n_3vec);
88 }
89 else {
90
91 const int* tab = s.siteTable().slice();
92
93 int totalSize = s.numSiteTable();
94
95 unordered_vaxpy3_y_user_arg arg(x, d, aptr, tab, 1);
96
97 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_y_evaluate_function);
98
100 // Original code
101
102 /*
103 for(int j=0; j < s.numSiteTable(); j++) {
104 int i=tab[j];
105 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
106 REAL* yptr = &(d.elem(i).elem(0).elem(0).real());
107 vaxpy3(yptr, aptr, xptr, yptr, Ns);
108 }
109 */
110 }
111}
112
113// d -= Scalar*Vec
114template<>
115inline
116void evaluate(OLattice< TVec >& d,
117 const OpSubtractAssign& op,
118 const QDPExpr<BinaryNode<OpMultiply,
119 Reference< QDPType< TScal, OScalar < TScal > > >,
120 Reference< QDPType< TVec, OLattice< TVec > > > >,
121 OLattice< TVec > > &rhs,
122 const Subset& s)
123{
124
125#ifdef DEBUG_BLAS
126 QDPIO::cout << "y -= a*x" << endl;
127#endif
128
129 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().right());
130 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().left());
131
132 // - sign as y -= ax <=> y = y-ax = -ax + y = axpy with -a
133 REAL ar = -( a.elem().elem().elem().elem());
134 REAL* aptr = &ar;
135 if( s.hasOrderedRep() ) {
136
137 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
138 REAL* yptr = &(d.elem(s.start()).elem(0).elem(0).real());
139
140 int total_n_3vec = (s.end()-s.start()+1);
141
142 ordered_vaxpy3_user_arg a = {yptr, aptr, xptr, yptr};
143
144 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
146 // Original code
148 //int n_3vec = (s.end()-s.start()+1)*Ns;
149 //vaxpy3(yptr, aptr, xptr, yptr, n_3vec);
150 }
151 else {
152
153 const int* tab = s.siteTable().slice();
154
155 int totalSize = s.numSiteTable();
160
162 // Original code
164 //for(int j=0; j < s.numSiteTable(); j++) {
165 //int i=tab[j];
166 //REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
167 //REAL* yptr = &(d.elem(i).elem(0).elem(0).real());
168 //vaxpy3(yptr, aptr, xptr, yptr, Ns);
169 //}
170 }
171
172}
173
174// z = ax + y
175template<>
176inline
177void evaluate( OLattice< TVec > &d,
178 const OpAssign &op,
179 const QDPExpr<
180 BinaryNode<OpAdd,
181 BinaryNode<OpMultiply,
182 Reference< QDPType< TScal, OScalar< TScal > > >,
183 Reference< QDPType< TVec, OLattice< TVec > > > >,
184 Reference< QDPType< TVec, OLattice< TVec > > > >,
185 OLattice< TVec > > &rhs,
186 const Subset& s)
187{
188
189#ifdef DEBUG_BLAS
190 QDPIO::cout << "z = a*x + y" << endl;
191#endif
192
193 // Peel the stuff out of the expression
194 // y is the right side of rhs
195 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right());
196
197 // ax is the left side of rhs and is in a binary node
198 typedef BinaryNode<OpMultiply,
199 Reference< QDPType< TScal, OScalar< TScal > > >,
200 Reference< QDPType< TVec, OLattice< TVec > > > > BN;
201
202 // get the binary node
203 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
204
205 // get a and x out of the bynary node
206 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
207 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
208 // Set pointers
209 REAL ar = a.elem().elem().elem().elem();
210 REAL *aptr = (REAL *)&ar;
211 if( s.hasOrderedRep() ) {
212 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
213 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
214 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
215
216 int total_n_3vec = (s.end()-s.start()+1);
217
218 ordered_vaxpy3_user_arg a = {zptr, aptr, xptr, yptr};
219
220 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
221
223 // Original code
225 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
226 //int n_3vec = (s.end()-s.start()+1)*Ns;
227 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
228 }
229 else {
230 const int* tab = s.siteTable().slice();
231
232 int totalSize = s.numSiteTable();
233
234 unordered_vaxpy3_z_user_arg arg(x, y, d, aptr, tab);
235
236 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_z_evaluate_function);
237
239 // Original code
241 /*
242 for(int j=0; j < s.numSiteTable(); j++) {
243 int i=tab[j];
244 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
245 REAL* yptr = (REAL *)&(y.elem(i).elem(0).elem(0).real());
246 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
247 vaxpy3(zptr, aptr, xptr, yptr, Ns);
248 }
249 */
250 }
251
252}
253
254
255// Vec = Vec + Scal*Vec
256template<>
257inline
258void evaluate( OLattice< TVec > &d,
259 const OpAssign &op,
260 const QDPExpr<
261 BinaryNode<OpAdd,
262 Reference< QDPType< TVec, OLattice< TVec > > >,
263 BinaryNode<OpMultiply,
264 Reference< QDPType< TScal, OScalar< TScal > > >,
265 Reference< QDPType< TVec, OLattice< TVec > > > > >,
266 OLattice< TVec > > &rhs,
267 const Subset& s)
268{
269#ifdef DEBUG_BLAS
270 QDPIO::cout << "z = y + a*x" << endl;
271#endif
272
273
274 // Peel the stuff out of the expression
275
276 // y is the left side of rhs
277 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().left());
278
279 // ax is the right side of rhs and is in a binary node
280 typedef BinaryNode<OpMultiply,
281 Reference< QDPType< TScal, OScalar< TScal > > >,
282 Reference< QDPType< TVec, OLattice< TVec > > > > BN;
283
284 // get the binary node
285 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
286
287 // get a and x out of the bynary node
288 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
289 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
290 // Set pointers
291 REAL ar = a.elem().elem().elem().elem();
292 REAL *aptr = (REAL *)&ar;
293 if( s.hasOrderedRep() ) {
294
295 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
296 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
297 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
298
299 int total_n_3vec = (s.end()-s.start()+1);
300
301 ordered_vaxpy3_user_arg a = {zptr, aptr, xptr, yptr};
302
303 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
304
306 // Original code
308 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
309 //int n_3vec = (s.end()-s.start()+1)*Ns;
310 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
311 }
312 else {
313
314 const int* tab = s.siteTable().slice();
315
316 int totalSize = s.numSiteTable();
317
318 unordered_vaxpy3_z_user_arg arg(x, y, d, aptr, tab);
319
320 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_z_evaluate_function);
321
323 // Original code
325 /*
326 for(int j=0; j < s.numSiteTable(); j++) {
327 int i=tab[j];
328 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
329 REAL* yptr = (REAL *)&(y.elem(i).elem(0).elem(0).real());
330 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
331 vaxpy3(zptr, aptr, xptr, yptr, Ns);
332 }
333 */
334 }
335
336}
337
338// Vec = Scalar*Vec - Vec
339template<>
340inline
341void evaluate( OLattice< TVec > &d,
342 const OpAssign &op,
343 const QDPExpr<
344 BinaryNode<OpSubtract,
345 BinaryNode<OpMultiply,
346 Reference< QDPType< TScal, OScalar< TScal > > >,
347 Reference< QDPType< TVec, OLattice< TVec > > > >,
348 Reference< QDPType< TVec, OLattice< TVec > > > >,
349 OLattice< TVec > > &rhs,
350 const Subset& s)
351{
352#ifdef DEBUG_BLAS
353 QDPIO::cout << "z = a*x - y" << endl;
354#endif
355
356
357 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right());
358
359 // ax is the left side of rhs and is in a binary node
360 typedef BinaryNode<OpMultiply,
361 Reference< QDPType< TScal, OScalar< TScal > > >,
362 Reference< QDPType< TVec, OLattice< TVec > > > > BN;
363
364 // get the binary node
365 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
366
367 // get a and x out of the bynary node
368 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
369 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
370 // Set pointers
371 REAL ar = a.elem().elem().elem().elem();
372 REAL *aptr = (REAL *)&ar;
373 if( s.hasOrderedRep() ) {
374 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
375 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
376 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
377
378 int total_n_3vec = (s.end()-s.start()+1);
379
380 ordered_vaxmy3_user_arg a = {zptr, aptr, xptr, yptr};
381
382 dispatch_to_threads(total_n_3vec, a, ordered_vaxmy3_evaluate_function);
383
385 // Original code
387 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
388 //int n_3vec = (s.end()-s.start()+1)*Ns;
389 //vaxmy3(zptr, aptr, xptr, yptr, n_3vec);
390 }
391 else {
392 const int* tab = s.siteTable().slice();
393
394 int totalSize = s.numSiteTable();
395
396 unordered_vaxmy3_z_user_arg arg(x, y, d, aptr, tab);
397
398 dispatch_to_threads(totalSize, arg, unordered_vaxmy3_z_evaluate_function);
399
401 // Original code
402
403 /*
404 for(int j=0; j < s.numSiteTable(); j++) {
405 int i=tab[j];
406 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
407 REAL* yptr = (REAL *)&(y.elem(i).elem(0).elem(0).real());
408 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
409 vaxmy3(zptr, aptr, xptr, yptr, Ns);
410 }
411 */
412 }
413
414}
415
416template<>
417inline
418void evaluate( OLattice< TVec > &d,
419 const OpAssign &op,
420 const QDPExpr<
421 BinaryNode<OpSubtract,
422 Reference< QDPType< TVec, OLattice< TVec > > >,
423 BinaryNode<OpMultiply,
424 Reference< QDPType< TScal, OScalar< TScal > > >,
425 Reference< QDPType< TVec, OLattice< TVec > > > > >,
426 OLattice< TVec > > &rhs,
427 const Subset& s)
428{
429#ifdef DEBUG_BLAS
430 QDPIO::cout << "z = y - a*x" << endl;
431#endif
432
433 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().left());
434
435 // ax is the right side of rhs and is in a binary node
436 typedef BinaryNode<OpMultiply,
437 Reference< QDPType< TScal, OScalar< TScal > > >,
438 Reference< QDPType< TVec, OLattice< TVec > > > > BN;
439
440 // get the binary node
441 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
442
443 // get a and x out of the bynary node
444 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.left());
445 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.right());
446 // Set pointers etc.
447
448 // -ve sign as y - ax = -ax + y = axpy with -a.
449 REAL ar = -a.elem().elem().elem().elem();
450 REAL *aptr = (REAL *)&ar;
451 if( s.hasOrderedRep() ) {
452
453 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
454 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
455 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
456
457 int total_n_3vec = (s.end()-s.start()+1);
458
459 ordered_vaxpy3_user_arg a = {zptr, aptr, xptr, yptr};
460
461 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
462
464 // Original code
466 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
467 //int n_3vec = (s.end()-s.start()+1)*Ns;
468 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
469 }
470 else {
471
472 const int* tab = s.siteTable().slice();
473
474 int totalSize = s.numSiteTable();
475
476 unordered_vaxpy3_z_user_arg arg(x, y, d, aptr, tab);
477
478 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_z_evaluate_function);
479
481 // Original code
483 /*
484 for(int j=0; j < s.numSiteTable(); j++) {
485 int i=tab[j];
486 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
487 REAL* yptr = (REAL *)&(y.elem(i).elem(0).elem(0).real());
488 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
489 vaxpy3(zptr, aptr, xptr, yptr, Ns);
490 }
491 */
492 }
493
494}
495
496// Vec += Vec * Scalar (AXPY)
497template<>
498inline
499void evaluate(OLattice< TVec >& d,
500 const OpAddAssign& op,
501 const QDPExpr<BinaryNode<OpMultiply,
502 Reference< QDPType< TVec, OLattice< TVec > > >,
503 Reference< QDPType< TScal, OScalar < TScal > > > >,
504 OLattice< TVec > > &rhs,
505 const Subset& s)
506{
507
508#ifdef DEBUG_BLAS
509 QDPIO::cout << "y += x*a" << endl;
510#endif
511
512 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().left());
513 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().right());
514
515 REAL ar = a.elem().elem().elem().elem();
516 REAL* aptr = &ar;
517
518 if( s.hasOrderedRep() ) {
519
520 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
521 REAL* yptr = &(d.elem(s.start()).elem(0).elem(0).real());
522 // cout << "Specialised axpy a ="<< ar << endl;
523
524 int total_n_3vec = (s.end()-s.start()+1);
525
526 ordered_vaxpy3_user_arg a = {yptr, aptr, xptr, yptr};
527
528 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
529
531 // Original code
533 //int n_3vec = (s.end()-s.start()+1)*Ns;
534 //vaxpy3(yptr, aptr, xptr, yptr, n_3vec);
535 }
536 else {
537
538 const int* tab = s.siteTable().slice();
539
540 int totalSize = s.numSiteTable();
541
542 unordered_vaxpy3_y_user_arg arg(x, d, aptr, tab, 1);
543
544 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_y_evaluate_function);
545
547 // Original code
549 /*
550 for(int j=0; j < s.numSiteTable(); j++) {
551 int i=tab[j];
552 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
553 REAL* yptr = &(d.elem(i).elem(0).elem(0).real());
554
555 vaxpy3(yptr, aptr, xptr, yptr, Ns);
556 }
557 */
558 }
559
560
561}
562
563
564// Vec -= Vec *Scalar
565template<>
566inline
567void evaluate(OLattice< TVec >& d,
568 const OpSubtractAssign& op,
569 const QDPExpr<BinaryNode<OpMultiply,
570 Reference< QDPType< TVec, OLattice< TVec > > >,
571 Reference< QDPType< TScal, OScalar < TScal > > > >,
572 OLattice< TVec > > &rhs,
573 const Subset& s)
574{
575
576#ifdef DEBUG_BLAS
577 QDPIO::cout << "y -= x*a" << endl;
578#endif
579
580 const OLattice< TVec >& x = static_cast<const OLattice< TVec > &>(rhs.expression().left());
581 const OScalar< TScal >& a = static_cast<const OScalar< TScal > &> (rhs.expression().right());
582
583 // - sign as y -= ax <=> y = y-ax = -ax + y = axpy with -a
584 REAL ar = -( a.elem().elem().elem().elem());
585 REAL* aptr = &ar;
586
587 if( s.hasOrderedRep() ) {
588
589 REAL* xptr = (REAL *)&(x.elem(s.start()).elem(0).elem(0).real());
590 REAL* yptr = &(d.elem(s.start()).elem(0).elem(0).real());
591
592 int total_n_3vec = (s.end()-s.start()+1);
593
594 ordered_vaxpy3_user_arg a = {yptr, aptr, xptr, yptr};
595
596 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
597
599 // Original code
601 //int n_3vec = (s.end()-s.start()+1)*Ns;
602 //vaxpy3(yptr, aptr, xptr, yptr, n_3vec);
603 }
604 else {
605
606 const int* tab = s.siteTable().slice();
607
608 int totalSize = s.numSiteTable();
609
610 unordered_vaxpy3_y_user_arg arg(x, d, aptr, tab, 1);
611
612 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_y_evaluate_function);
613
615 // Original code
617 /*
618 for(int j=0; j < s.numSiteTable(); j++) {
619 int i=tab[j];
620 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
621 REAL* yptr = &(d.elem(i).elem(0).elem(0).real());
622
623 vaxpy3(yptr, aptr, xptr, yptr, Ns);
624 }
625 */
626 }
627
628}
629
630
631// Vec = Vec *Scalar + Vec (AXPY)
632template<>
633inline
634void evaluate( OLattice< TVec > &d,
635 const OpAssign &op,
636 const QDPExpr<
637 BinaryNode<OpAdd,
638 BinaryNode<OpMultiply,
639 Reference< QDPType< TVec, OLattice< TVec > > >,
640 Reference< QDPType< TScal, OScalar< TScal > > > >,
641 Reference< QDPType< TVec, OLattice< TVec > > > >,
642 OLattice< TVec > > &rhs,
643 const Subset& s)
644{
645
646#ifdef DEBUG_BLAS
647 QDPIO::cout << "z = x*a + y" << endl;
648#endif
649
650 // Peel the stuff out of the expression
651 // y is the right side of rhs
652 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right());
653
654 // ax is the right side of rhs and is in a binary node
655 typedef BinaryNode<OpMultiply,
656 Reference< QDPType< TVec, OLattice< TVec > > >,
657 Reference< QDPType< TScal, OScalar< TScal > > > > BN;
658
659 // get the binary node
660 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
661
662 // get a and x out of the bynary node
663 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.right());
664 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.left());
665 // Set pointers
666 REAL ar = a.elem().elem().elem().elem();
667 REAL *aptr = (REAL *)&ar;
668 if( s.hasOrderedRep() ) {
669
670 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
671 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
672 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
673
674 int total_n_3vec = (s.end()-s.start()+1);
675
676 ordered_vaxpy3_user_arg a = {zptr, aptr, xptr, yptr};
677
678 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
679
681 // Original code
683 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
684 //int n_3vec = (s.end()-s.start()+1)*Ns;
685 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
686 }
687 else {
688
689 const int* tab = s.siteTable().slice();
690
691 int totalSize = s.numSiteTable();
692
693 unordered_vaxpy3_z_user_arg arg(x, y, d, aptr, tab);
694
695 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_z_evaluate_function);
696
698 // Original code
700 /*
701 for(int j=0; j < s.numSiteTable(); j++) {
702 int i=tab[j];
703
704 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
705 REAL* yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
706 REAL* zptr = (REAL *) &(d.elem(i).elem(0).elem(0).real());
707 vaxpy3(zptr, aptr, xptr, yptr, Ns);
708 }
709 */
710 }
711
712}
713
714
715// Vec = Vec + Vec * Scalar (AXPY)
716template<>
717inline
718void evaluate( OLattice< TVec > &d,
719 const OpAssign &op,
720 const QDPExpr<
721 BinaryNode<OpAdd,
722 Reference< QDPType< TVec, OLattice< TVec > > >,
723 BinaryNode<OpMultiply,
724 Reference< QDPType< TVec, OLattice< TVec > > >,
725 Reference< QDPType< TScal, OScalar< TScal > > > > >,
726 OLattice< TVec > > &rhs,
727 const Subset& s)
728{
729#ifdef DEBUG_BLAS
730 QDPIO::cout << "z = y + x*a" << endl;
731#endif
732
733
734 // Peel the stuff out of the expression
735
736 // y is the left side of rhs
737 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().left());
738
739 // ax is the right side of rhs and is in a binary node
740 typedef BinaryNode<OpMultiply,
741 Reference< QDPType< TVec, OLattice< TVec > > >,
742 Reference< QDPType< TScal, OScalar< TScal > > > > BN;
743
744 // get the binary node
745 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
746
747 // get a and x out of the bynary node
748 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.right());
749 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.left());
750 // Set pointers
751 REAL ar = a.elem().elem().elem().elem();
752 REAL *aptr = (REAL *)&ar;
753
754 if( s.hasOrderedRep() ) {
755
756 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
757 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
758 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
759
760 int total_n_3vec = (s.end()-s.start()+1);
761
762 ordered_vaxpy3_user_arg a = {zptr, aptr, xptr, yptr};
763
764 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
765
767 // Original code
769 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
770 //int n_3vec = (s.end()-s.start()+1)*Ns;
771 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
772 }
773 else {
774
775 const int* tab = s.siteTable().slice();
776
777 int totalSize = s.numSiteTable();
778
779 unordered_vaxpy3_z_user_arg arg(x, y, d, aptr, tab);
780
781 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_z_evaluate_function);
782
784 // Original code
786 /*
787 for(int j=0; j < s.numSiteTable(); j++) {
788 int i=tab[j];
789
790 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
791 REAL* yptr = (REAL *)&(y.elem(i).elem(0).elem(0).real());
792 REAL* zptr = (REAL *)&(d.elem(i).elem(0).elem(0).real());
793 vaxpy3(zptr, aptr, xptr, yptr, Ns);
794 }
795 */
796 }
797
798}
799
800
801// Vec = Vec*Scalar - Vec (AXMY)
802template<>
803inline
804void evaluate( OLattice< TVec > &d,
805 const OpAssign &op,
806 const QDPExpr<
807 BinaryNode<OpSubtract,
808 BinaryNode<OpMultiply,
809 Reference< QDPType< TVec, OLattice< TVec > > >,
810 Reference< QDPType< TScal, OScalar< TScal > > > >,
811 Reference< QDPType< TVec, OLattice< TVec > > > >,
812 OLattice< TVec > > &rhs,
813 const Subset& s)
814{
815#ifdef DEBUG_BLAS
816 QDPIO::cout << "z = x*a - y" << endl;
817#endif
818
819 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().right());
820
821
822 typedef BinaryNode<OpMultiply,
823 Reference< QDPType< TVec, OLattice< TVec > > >,
824 Reference< QDPType< TScal, OScalar< TScal > > > > BN;
825
826 // get the binary node
827 const BN &mulNode = static_cast<const BN&> (rhs.expression().left());
828
829 // get a and x out of the bynary node
830 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.right());
831 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.left());
832 // Set pointers
833 REAL ar = a.elem().elem().elem().elem();
834 REAL *aptr = (REAL *)&ar;
835
836 if( s.hasOrderedRep() ) {
837 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
838 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
839 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
840
841 int total_n_3vec = (s.end()-s.start()+1);
842
843 ordered_vaxmy3_user_arg a = {zptr, aptr, xptr, yptr, };
844
845 dispatch_to_threads(total_n_3vec, a, ordered_vaxmy3_evaluate_function);
846
848 // Original code
850 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
851 //int n_3vec = (s.end()-s.start()+1)*Ns;
852 //vaxmy3(zptr, aptr, xptr, yptr, n_3vec);
853 }
854 else {
855 const int* tab = s.siteTable().slice();
856
857 int totalSize = s.numSiteTable();
858
859 unordered_vaxmy3_z_user_arg arg(x, y, d, aptr, tab);
860
861 dispatch_to_threads(totalSize, arg, unordered_vaxmy3_z_evaluate_function);
862
864 // Original code
866 /*for(int j=0; j < s.numSiteTable(); j++) {
867 int i=tab[j];
868
869 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
870 REAL* yptr = (REAL *)&(y.elem(i).elem(0).elem(0).real());
871 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
872 vaxmy3(zptr, aptr, xptr, yptr, Ns);
873 }
874 */
875 }
876
877}
878
879
880// Vec = Vec - Vec*Scalar (AXPY with -Scalar)
881template<>
882inline
883void evaluate( OLattice< TVec > &d,
884 const OpAssign &op,
885 const QDPExpr<
886 BinaryNode<OpSubtract,
887 Reference< QDPType< TVec, OLattice< TVec > > >,
888 BinaryNode<OpMultiply,
889 Reference< QDPType< TVec, OLattice< TVec > > >,
890 Reference< QDPType< TScal, OScalar< TScal > > > > >,
891 OLattice< TVec > > &rhs,
892 const Subset& s)
893{
894#ifdef DEBUG_BLAS
895 QDPIO::cout << "z = y - x*a" << endl;
896#endif
897
898 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&> (rhs.expression().left());
899
900 typedef BinaryNode<OpMultiply,
901 Reference< QDPType< TVec, OLattice< TVec > > >,
902 Reference< QDPType< TScal, OScalar< TScal > > > > BN;
903
904 // get the binary node
905 const BN &mulNode = static_cast<const BN&> (rhs.expression().right());
906
907 // get a and x out of the bynary node
908 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode.right());
909 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode.left());
910 // Set pointers etc.
911
912 // -ve sign as y - ax = -ax + y = axpy with -a.
913 REAL ar = -a.elem().elem().elem().elem();
914 REAL *aptr = (REAL *)&ar;
915
916 if( s.hasOrderedRep() ) {
917
918 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
919 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
920 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
921
922 int total_n_3vec = (s.end()-s.start()+1);
923
924 ordered_vaxpy3_user_arg a = {zptr, aptr, xptr, yptr};
925
926 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
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)*Ns;
933 //vaxpy3(zptr, aptr, xptr, yptr, n_3vec);
934 }
935 else {
936
937 const int* tab = s.siteTable().slice();
938
939 int totalSize = s.numSiteTable();
940
941 unordered_vaxpy3_z_user_arg arg(x, y, d, aptr, tab);
942
943 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_z_evaluate_function);
944
946 // Original code
948 /*
949 for(int j=0; j < s.numSiteTable(); j++) {
950 int i=tab[j];
951
952 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
953 REAL* yptr = (REAL *)&(y.elem(i).elem(0).elem(0).real());
954 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
955 vaxpy3(zptr, aptr, xptr, yptr, Ns);
956 }
957 */
958 }
959
960}
961
962
963template<>
964inline
965void evaluate( OLattice< TVec > &d,
966 const OpAssign &op,
967 const QDPExpr<
968 BinaryNode<OpAdd,
969 Reference< QDPType< TVec, OLattice< TVec > > >,
970 Reference< QDPType< TVec, OLattice< TVec > > > >,
971 OLattice< TVec > > &rhs,
972 const Subset& s)
973{
974#ifdef DEBUG_BLAS
975 cout << "BJ: v+v " << endl;
976#endif
977
978 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(rhs.expression().left());
979 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(rhs.expression().right());
980
981 REAL one = 1;
982
983 if( s.hasOrderedRep() ) {
984
985 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
986 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
987 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
988
989 int total_n_3vec = (s.end()-s.start()+1);
990
991 ordered_vaxpy3_user_arg a = {zptr, &one, xptr, yptr};
992
993 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
994
996 // Original code
998 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
999 //int n_3vec = (s.end()-s.start()+1)*Ns;
1000 //vaxpy3(zptr,&one, xptr, yptr, n_3vec);
1001 }
1002 else {
1003
1004 const int* tab = s.siteTable().slice();
1005
1006 int totalSize = s.numSiteTable();
1007
1008 unordered_vaxpy3_z_user_arg arg(x, y, d, &one, tab);
1009
1010 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_z_evaluate_function);
1011
1013 // Original code
1015 /*
1016 for(int j=0; j < s.numSiteTable(); j++) {
1017 int i=tab[j];
1018
1019 REAL* xptr = (REAL *)&(x.elem(i).elem(0).elem(0).real());
1020 REAL* yptr = (REAL *)&(y.elem(i).elem(0).elem(0).real());
1021 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1022 vaxpy3(zptr,&one, xptr, yptr, Ns);
1023 }
1024 */
1025 }
1026
1027
1028}
1029
1030template<>
1031inline
1032void evaluate( OLattice< TVec > &d,
1033 const OpAssign &op,
1034 const QDPExpr<
1035 BinaryNode<OpSubtract,
1036 Reference< QDPType< TVec, OLattice< TVec > > >,
1037 Reference< QDPType< TVec, OLattice< TVec > > > >,
1038 OLattice< TVec > > &rhs,
1039 const Subset& s)
1040{
1041#ifdef DEBUG_BLAS
1042 cout << "BJ: v-v " << endl;
1043#endif
1044
1045 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(rhs.expression().left());
1046 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(rhs.expression().right());
1047 REAL one=1;
1048
1049 if( s.hasOrderedRep() ) {
1050
1051 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1052 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1053 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1054
1055 int total_n_3vec = (s.end()-s.start()+1);
1056
1057 ordered_vaxmy3_user_arg a = {zptr, &one, xptr, yptr};
1058
1059 dispatch_to_threads(total_n_3vec, a, ordered_vaxmy3_evaluate_function);
1060
1062 // Original code
1064 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1065 //int n_3vec = (s.end()-s.start()+1)*Ns;
1066
1067 //vaxmy3(zptr,&one, xptr, yptr, n_3vec);
1068 }
1069 else {
1070 const int* tab = s.siteTable().slice();
1071
1072 int totalSize = s.numSiteTable();
1073
1074 unordered_vaxmy3_z_user_arg arg(x, y, d, &one, tab);
1075
1076 dispatch_to_threads(totalSize, arg, unordered_vaxmy3_z_evaluate_function);
1077
1079 // Original code
1081 /*for(int j=0; j < s.numSiteTable(); j++) {
1082 int i=tab[j];
1083 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1084 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1085 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1086
1087 vaxmy3(zptr,&one, xptr, yptr, Ns);
1088
1089 }*/
1090 }
1091
1092}
1093
1094// Vec = Scal * Vec
1095template<>
1096inline
1097void evaluate( OLattice< TVec > &d,
1098 const OpAssign &op,
1099 const QDPExpr<
1100 BinaryNode<OpMultiply,
1101 Reference< QDPType< TScal, OScalar< TScal > > >,
1102 Reference< QDPType< TVec, OLattice< TVec > > > >,
1103 OLattice< TVec > > &rhs,
1104 const Subset& s)
1105{
1106#ifdef DEBUG_BLAS
1107 cout << "BJ: v = a*v " << endl;
1108#endif
1109 const OLattice< TVec > &x = static_cast<const OLattice< TVec >&>(rhs.expression().right());
1110 const OScalar< TScal > &a = static_cast<const OScalar< TScal >&>(rhs.expression().left());
1111
1112 REAL ar = a.elem().elem().elem().elem();
1113 REAL *aptr = &ar;
1114
1115 if( s.hasOrderedRep() ) {
1116
1117 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1118 REAL *zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1119
1120 int total_n_3vec = (s.end()-s.start()+1);
1121
1122 ordered_vscal_user_arg a = {zptr, aptr, xptr};
1123
1124 dispatch_to_threads(total_n_3vec, a, ordered_vscal_evaluate_function);
1125
1127 // Original code
1129 //int n_3vec = (s.end()-s.start()+1)*Ns;
1130 //vscal(zptr, aptr, xptr, n_3vec);
1131 }
1132 else {
1133 const int* tab = s.siteTable().slice();
1134
1135 int totalSize = s.numSiteTable();
1136
1137 unordered_vscal_user_arg arg(x, d, aptr, tab);
1138
1139 dispatch_to_threads(totalSize, arg, unordered_vscal_evaluate_function);
1140
1142 // Original code
1144 /*for(int j=0; j < s.numSiteTable(); j++) {
1145 int i=tab[j];
1146 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1147 REAL *zptr = &(d.elem(i).elem(0).elem(0).real());
1148
1149 vscal(zptr, aptr, xptr, Ns);
1150 }*/
1151 }
1152
1153}
1154
1155template<>
1156inline
1157void evaluate( OLattice< TVec > &d,
1158 const OpAssign &op,
1159 const QDPExpr<
1160 BinaryNode<OpMultiply,
1161 Reference< QDPType< TVec, OLattice< TVec > > >,
1162 Reference< QDPType< TScal, OScalar< TScal > > > >,
1163 OLattice< TVec > > &rhs,
1164 const Subset& s)
1165{
1166#ifdef DEBUG_BLAS
1167 cout << "BJ: v = v*a " << endl;
1168#endif
1169
1170 const OLattice< TVec > &x = static_cast<const OLattice< TVec >&>(rhs.expression().left());
1171 const OScalar< TScal > &a = static_cast<const OScalar< TScal >&>(rhs.expression().right());
1172
1173 REAL ar = a.elem().elem().elem().elem();
1174 REAL *aptr = &ar;
1175
1176 if( s.hasOrderedRep() ) {
1177 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1178 REAL *zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1179
1180 int total_n_3vec = (s.end()-s.start()+1);
1181
1182 ordered_vscal_user_arg a = {zptr, aptr, xptr};
1183
1184 dispatch_to_threads(total_n_3vec, a, ordered_vscal_evaluate_function);
1185
1187 // Original code
1189 //int n_3vec = (s.end()-s.start()+1)*Ns;
1190 //vscal(zptr, aptr, xptr, n_3vec);
1191 }
1192 else {
1193 const int* tab = s.siteTable().slice();
1194
1195 int totalSize = s.numSiteTable();
1196
1197 unordered_vscal_user_arg arg(x, d, aptr, tab);
1198
1199 dispatch_to_threads(totalSize, arg, unordered_vscal_evaluate_function);
1200
1202 // Original code
1204 /*for(int j=0; j < s.numSiteTable(); j++) {
1205 int i=tab[j];
1206 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1207 REAL *zptr = &(d.elem(i).elem(0).elem(0).real());
1208
1209 vscal(zptr, aptr, xptr, Ns);
1210 }*/
1211 }
1212}
1213
1214// v *= a
1215template<>
1216inline
1217void evaluate( OLattice< TVec > &d,
1218 const OpMultiplyAssign &op,
1219 const QDPExpr<
1220 UnaryNode<OpIdentity,
1221 Reference< QDPType< TScal, OScalar< TScal > > > >,
1222 OScalar< TScal > > &rhs,
1223 const Subset& s)
1224{
1225 const OScalar< TScal >& a = static_cast< const OScalar<TScal >&>(rhs.expression().child());
1226
1227
1228#ifdef DEBUG_BLAS
1229 QDPIO::cout << "BJ: v *= a, a = " << a << endl;
1230#endif
1231
1232 REAL ar = a.elem().elem().elem().elem();
1233 if( s.hasOrderedRep() ) {
1234
1235 REAL* xptr = &(d.elem(s.start()).elem(0).elem(0).real());
1236 REAL* zptr = xptr;
1237
1238 int total_n_3vec = (s.end()-s.start()+1);
1239
1240 ordered_vscal_user_arg a = {zptr, &ar, xptr};
1241
1242 dispatch_to_threads(total_n_3vec, a, ordered_vscal_evaluate_function);
1243
1245 // Original code
1247 //int n_3vec = (s.end()-s.start()+1)*Ns;
1248 //vscal(zptr,&ar, xptr, n_3vec);
1249 }
1250 else {
1251 const int* tab = s.siteTable().slice();
1252
1253 int totalSize = s.numSiteTable();
1254
1255 unordered_vscal_user_arg arg(d, d, &ar, tab);
1256
1257 dispatch_to_threads(totalSize, arg, unordered_vscal_evaluate_function);
1258
1260 // Original code
1262 /*for(int j=0; j < s.numSiteTable(); j++) {
1263 int i=tab[j];
1264
1265 REAL* xptr = &(d.elem(i).elem(0).elem(0).real());
1266 REAL* zptr = xptr;
1267
1268 vscal(zptr,&ar, xptr, Ns);
1269 }*/
1270 }
1271}
1272
1273// v /= a
1274template<>
1275inline
1276void evaluate( OLattice< TVec > &d,
1277 const OpDivideAssign &op,
1278 const QDPExpr<
1279 UnaryNode<OpIdentity,
1280 Reference< QDPType< TScal, OScalar< TScal > > > >,
1281 OScalar< TScal > > &rhs,
1282 const Subset& s)
1283{
1284 const OScalar< TScal >& a = static_cast< const OScalar<TScal >&>(rhs.expression().child());
1285
1286
1287#ifdef DEBUG_BLAS
1288 QDPIO::cout << "BJ: v /= a, a = " << a << endl;
1289#endif
1290
1291 REAL ar = (REAL)1/a.elem().elem().elem().elem();
1292 if( s.hasOrderedRep() ) {
1293 REAL* xptr = &(d.elem(s.start()).elem(0).elem(0).real());
1294 REAL* zptr = xptr;
1295
1296 int total_n_3vec = (s.end()-s.start()+1);
1297
1298 ordered_vscal_user_arg a = {zptr, &ar, xptr};
1299
1300 dispatch_to_threads(total_n_3vec, a, ordered_vscal_evaluate_function);
1301
1303 // Original code
1305 //int n_3vec = (s.end()-s.start()+1)*Ns;
1306 //vscal(zptr,&ar, xptr, n_3vec);
1307 }
1308 else {
1309 const int* tab = s.siteTable().slice();
1310
1311 int totalSize = s.numSiteTable();
1312
1313 unordered_vscal_user_arg arg(d, d, &ar, tab);
1314
1315 dispatch_to_threads(totalSize, arg, unordered_vscal_evaluate_function);
1316
1318 // Original code
1320 /*for(int j=0; j < s.numSiteTable(); j++) {
1321 int i=tab[j];
1322
1323 REAL* xptr = &(d.elem(i).elem(0).elem(0).real());
1324 REAL* zptr = xptr;
1325
1326 vscal(zptr,&ar, xptr, Ns);
1327 }*/
1328 }
1329}
1330
1331// v += v
1332template<>
1333inline
1334void evaluate( OLattice< TVec > &d,
1335 const OpAddAssign &op,
1336 const QDPExpr<
1337 UnaryNode<OpIdentity,
1338 Reference< QDPType< TVec, OLattice< TVec > > > >,
1339 OLattice< TVec > > &rhs,
1340 const Subset& s)
1341{
1342 const OLattice< TVec >& x = static_cast< const OLattice<TVec >&>(rhs.expression().child());
1343
1344
1345
1346#ifdef DEBUG_BLAS
1347 QDPIO::cout << "BJ: v += v" << endl;
1348#endif
1349 REAL one = 1;
1350
1351 if( s.hasOrderedRep() ) {
1352 //int n_3vec = (s.end() - s.start()+1)*Ns;
1353 REAL *xptr = (REAL *)(&x.elem(s.start()).elem(0).elem(0).real());
1354 REAL *yptr = (REAL *)(&d.elem(s.start()).elem(0).elem(0).real());
1355
1356
1357 int total_n_3vec = (s.end()-s.start()+1);
1358
1359 ordered_vaxpy3_user_arg a = {yptr, &one, yptr, xptr};
1360
1361 dispatch_to_threads(total_n_3vec, a, ordered_vaxpy3_evaluate_function);
1362
1364 // Original code
1366 //int n_3vec = (s.end() - s.start()+1)*Ns;
1367 //vaxpy3(yptr, &one, yptr, xptr,n_3vec);
1368 }
1369 else {
1370
1371 const int* tab = s.siteTable().slice();
1372
1373 int totalSize = s.numSiteTable();
1374
1375 unordered_vaxpy3_y_user_arg arg(x, d, &one, tab, 0);
1376
1377 dispatch_to_threads(totalSize, arg, unordered_vaxpy3_y_evaluate_function);
1378
1380 // Original code
1382 /*
1383 for(int j=0; j < s.numSiteTable(); j++) {
1384 int i=tab[j];
1385
1386 REAL *xptr = (REAL *)(&x.elem(i).elem(0).elem(0).real());
1387 REAL *yptr = (REAL *)(&d.elem(i).elem(0).elem(0).real());
1388
1389 vaxpy3(yptr, &one, yptr, xptr,Ns);// yptr and xptr change place
1390
1391 }
1392 */
1393 }
1394
1395}
1396
1397// v -= v
1398template<>
1399inline
1400void evaluate( OLattice< TVec > &d,
1401 const OpSubtractAssign &op,
1402 const QDPExpr<
1403 UnaryNode<OpIdentity,
1404 Reference< QDPType< TVec, OLattice< TVec > > > >,
1405 OLattice< TVec > > &rhs,
1406 const Subset& s)
1407{
1408 const OLattice< TVec >& x = static_cast< const OLattice<TVec >&>(rhs.expression().child());
1409
1410
1411
1412#ifdef DEBUG_BLAS
1413 QDPIO::cout << "BJ: v -= v" << endl;
1414#endif
1415 REAL one = 1;
1416
1417 if( s.hasOrderedRep() ) {
1418
1419 REAL *xptr = (REAL *)(&x.elem(s.start()).elem(0).elem(0).real());
1420 REAL *yptr = (REAL *)(&d.elem(s.start()).elem(0).elem(0).real());
1421
1422 int total_n_3vec = (s.end()-s.start()+1);
1423
1424 ordered_vaxmy3_user_arg a = {yptr, &one, yptr, xptr};
1425
1426 dispatch_to_threads(total_n_3vec, a, ordered_vaxmy3_evaluate_function);
1427
1429 // Original code
1431 //int n_3vec = (s.end() - s.start()+1)*Ns;
1432 //vaxmy3(yptr, &one, yptr, xptr, n_3vec);
1433 }
1434 else {
1435 const int* tab = s.siteTable().slice();
1436
1437 int totalSize = s.numSiteTable();
1438
1439 unordered_vaxmy3_y_user_arg arg(x, d, &one, tab);
1440
1441 dispatch_to_threads(totalSize, arg, unordered_vaxmy3_y_evaluate_function);
1442
1444 // Original code
1446 /*for(int j=0; j < s.numSiteTable(); j++) {
1447 int i=tab[j];
1448 REAL *xptr = (REAL *)(&x.elem(i).elem(0).elem(0).real());
1449 REAL *yptr = (REAL *)(&d.elem(i).elem(0).elem(0).real());
1450
1451 vaxmy3(yptr, &one, yptr, xptr, Ns);
1452
1453 }*/
1454 }
1455
1456}
1457
1458
1459// z = ax + by
1460template<>
1461inline
1462void evaluate( OLattice< TVec > &d,
1463 const OpAssign &op,
1464 const QDPExpr<
1465 BinaryNode<OpAdd,
1466 BinaryNode<OpMultiply,
1467 Reference< QDPType< TScal, OScalar< TScal > > >,
1468 Reference< QDPType< TVec, OLattice< TVec > > > >,
1469 BinaryNode<OpMultiply,
1470 Reference< QDPType< TScal, OScalar< TScal > > >,
1471 Reference< QDPType< TVec, OLattice< TVec > > > > >,
1472 OLattice< TVec > > &rhs,
1473 const Subset& s)
1474{
1475
1476#ifdef DEBUG_BLAS
1477 QDPIO::cout << "z = a*x + b*y" << endl;
1478#endif
1479
1480 // Peel the stuff out of the expression
1481 // y is the right side of rhs
1482
1483 // ax is the left side of rhs and is in a binary node
1484 typedef BinaryNode<OpMultiply,
1485 Reference< QDPType< TScal, OScalar< TScal > > >,
1486 Reference< QDPType< TVec, OLattice< TVec > > > > BN;
1487
1488 // get the binary node
1489 const BN &mulNode1 = static_cast<const BN&> (rhs.expression().left());
1490 const BN &mulNode2 = static_cast<const BN&> (rhs.expression().right());
1491
1492 // get a and x out of the binary node
1493 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
1494 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
1495
1496 // get b and y out of the binary node
1497 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1498 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right());
1499
1500
1501 // Set pointers
1502 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1503 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1504
1505 if( s.hasOrderedRep() ) {
1506 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1507 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1508 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1509
1510 int total_n_3vec = (s.end()-s.start()+1);
1511
1512 ordered_vaxpby3_user_arg a = {zptr, aptr, xptr, bptr, yptr};
1513
1514 dispatch_to_threads(total_n_3vec, a, ordered_vaxpby3_evaluate_function);
1515
1517 // Original code
1519 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1520 //int n_3vec = (s.end()-s.start()+1)*Ns;
1521 //vaxpby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1522 }
1523 else {
1524 const int* tab = s.siteTable().slice();
1525
1526 int totalSize = s.numSiteTable();
1527
1528 unordered_vaxpby3_user_arg arg(x, y, d, aptr, bptr, tab);
1529
1530 dispatch_to_threads(totalSize, arg, unordered_vaxpby3_evaluate_function);
1531
1533 // Original code
1535 /*for(int j=0; j < s.numSiteTable(); j++) {
1536 int i=tab[j];
1537 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1538 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1539 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1540
1541 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1542 vaxpby3(zptr, aptr, xptr, bptr, yptr, Ns);
1543
1544 }*/
1545 }
1546
1547}
1548
1549
1550// z = xa + by
1551template<>
1552inline
1553void evaluate( OLattice< TVec > &d,
1554 const OpAssign &op,
1555 const QDPExpr<
1556 BinaryNode<OpAdd,
1557 BinaryNode<OpMultiply,
1558 Reference< QDPType< TVec, OLattice< TVec > > >,
1559 Reference< QDPType< TScal, OScalar< TScal > > > >,
1560 BinaryNode<OpMultiply,
1561 Reference< QDPType< TScal, OScalar< TScal > > >,
1562 Reference< QDPType< TVec, OLattice< TVec > > > > >,
1563 OLattice< TVec > > &rhs,
1564 const Subset& s)
1565{
1566
1567#ifdef DEBUG_BLAS
1568 QDPIO::cout << "z = x*a + b*y" << endl;
1569#endif
1570
1571 // Peel the stuff out of the expression
1572 // y is the right side of rhs
1573
1574 // ax is the left side of rhs and is in a binary node
1575 typedef BinaryNode<OpMultiply,
1576 Reference< QDPType< TVec, OLattice< TVec > > >,
1577 Reference< QDPType< TScal, OScalar< TScal > > > > BN1;
1578
1579 typedef BinaryNode<OpMultiply,
1580 Reference< QDPType< TScal, OScalar< TScal > > >,
1581 Reference< QDPType< TVec, OLattice< TVec > > > > BN2;
1582
1583
1584
1585 // get the binary node
1586 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
1587 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
1588
1589 // get a and x out of the binary node
1590 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.left());
1591
1592 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.right());
1593
1594 // get b and y out of the binary node
1595 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1596 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right());
1597
1598
1599 // Set pointers
1600 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1601 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1602
1603 if( s.hasOrderedRep() ) {
1604 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1605 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1606 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1607
1608 int total_n_3vec = (s.end()-s.start()+1);
1609
1610 ordered_vaxpby3_user_arg a = {zptr, aptr, xptr, bptr, yptr};
1611
1612 dispatch_to_threads(total_n_3vec, a, ordered_vaxpby3_evaluate_function);
1613
1615 // Original code
1617 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1618 //int n_3vec = (s.end()-s.start()+1)*Ns;
1619 //vaxpby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1620 }
1621 else {
1622 const int* tab = s.siteTable().slice();
1623
1624 int totalSize = s.numSiteTable();
1625
1626 unordered_vaxpby3_user_arg arg(x, y, d, aptr, bptr, tab);
1627
1628 dispatch_to_threads(totalSize, arg, unordered_vaxpby3_evaluate_function);
1629
1631 // Original code
1633 /*for(int j=0; j < s.numSiteTable(); j++) {
1634 int i=tab[j];
1635
1636 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1637 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1638 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1639
1640 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1641 vaxpby3(zptr, aptr, xptr, bptr, yptr, Ns);
1642
1643 }*/
1644 }
1645}
1646
1647// z = ax + yb
1648template<>
1649inline
1650void evaluate( OLattice< TVec > &d,
1651 const OpAssign &op,
1652 const QDPExpr<
1653 BinaryNode<OpAdd,
1654 BinaryNode<OpMultiply,
1655 Reference< QDPType< TScal, OScalar< TScal > > >,
1656 Reference< QDPType< TVec, OLattice< TVec > > > >,
1657 BinaryNode<OpMultiply,
1658 Reference< QDPType< TVec, OLattice< TVec > > >,
1659 Reference< QDPType< TScal, OScalar< TScal > > > > >,
1660 OLattice< TVec > > &rhs,
1661 const Subset& s)
1662{
1663
1664#ifdef DEBUG_BLAS
1665 QDPIO::cout << "z = a*x + y*b" << endl;
1666#endif
1667
1668 // Peel the stuff out of the expression
1669 // y is the right side of rhs
1670
1671 // type of a*x
1672 typedef BinaryNode<OpMultiply,
1673 Reference< QDPType< TScal, OScalar< TScal > > >,
1674 Reference< QDPType< TVec, OLattice< TVec > > > > BN1;
1675
1676 // type of y*b
1677 typedef BinaryNode<OpMultiply,
1678 Reference< QDPType< TVec, OLattice< TVec > > >,
1679 Reference< QDPType< TScal, OScalar< TScal > > > > BN2;
1680
1681
1682
1683 // get the binary nodes
1684 // a*x node
1685 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
1686
1687 // y*b node
1688 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
1689
1690 // get a and x out of the binary node
1691 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
1692
1693 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
1694
1695
1696 // get b and y out of the binary node
1697 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.left());
1698
1699 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.right());
1700
1701
1702 // Set pointers
1703 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1704 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1705
1706 if( s.hasOrderedRep() ) {
1707 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1708 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1709 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1710
1711 int total_n_3vec = (s.end()-s.start()+1);
1712
1713 ordered_vaxpby3_user_arg a = {zptr, aptr, xptr, bptr, yptr};
1714
1715 dispatch_to_threads(total_n_3vec, a, ordered_vaxpby3_evaluate_function);
1716
1718 // Original code
1720 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1721 //int n_3vec = (s.end()-s.start()+1)*Ns;
1722 //vaxpby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1723 }
1724 else {
1725 const int* tab = s.siteTable().slice();
1726
1727 int totalSize = s.numSiteTable();
1728
1729 unordered_vaxpby3_user_arg arg(x, y, d, aptr, bptr, tab);
1730
1731 dispatch_to_threads(totalSize, arg, unordered_vaxpby3_evaluate_function);
1732
1734 // Original code
1736 /*for(int j=0; j < s.numSiteTable(); j++) {
1737 int i=tab[j];
1738
1739 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1740 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1741 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1742
1743 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1744 vaxpby3(zptr, aptr, xptr, bptr, yptr, Ns);
1745
1746 }*/
1747 }
1748}
1749
1750// z = xa + yb
1751template<>
1752inline
1753void evaluate( OLattice< TVec > &d,
1754 const OpAssign &op,
1755 const QDPExpr<
1756 BinaryNode<OpAdd,
1757 BinaryNode<OpMultiply,
1758 Reference< QDPType< TVec, OLattice< TVec > > >,
1759 Reference< QDPType< TScal, OScalar< TScal > > > >,
1760 BinaryNode<OpMultiply,
1761 Reference< QDPType< TVec, OLattice< TVec > > >,
1762 Reference< QDPType< TScal, OScalar< TScal > > > > >,
1763 OLattice< TVec > > &rhs,
1764 const Subset& s)
1765{
1766
1767#ifdef DEBUG_BLAS
1768 QDPIO::cout << "z = x*a + y*b" << endl;
1769#endif
1770
1771 // Peel the stuff out of the expression
1772 // y is the right side of rhs
1773
1774 // ax is the left side of rhs and is in a binary node
1775 typedef BinaryNode<OpMultiply,
1776 Reference< QDPType< TVec, OLattice< TVec > > >,
1777 Reference< QDPType< TScal, OScalar< TScal > > > > BN;
1778
1779 // get the binary node
1780 const BN &mulNode1 = static_cast<const BN&> (rhs.expression().left());
1781 const BN &mulNode2 = static_cast<const BN&> (rhs.expression().right());
1782
1783 // get a and x out of the binary node
1784 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.left());
1785 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.right());
1786
1787 // get b and y out of the binary node
1788 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.left());
1789
1790 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.right());
1791
1792 // Set pointers
1793 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1794 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1795
1796 if( s.hasOrderedRep() ) {
1797 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1798 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1799 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1800
1801 int total_n_3vec = (s.end()-s.start()+1);
1802
1803 ordered_vaxpby3_user_arg a = {zptr, aptr, xptr, bptr, yptr};
1804
1805 dispatch_to_threads(total_n_3vec, a, ordered_vaxpby3_evaluate_function);
1806
1808 // Original code
1810 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1811 //int n_3vec = (s.end()-s.start()+1)*Ns;
1812 //vaxpby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1813 }
1814 else {
1815 const int* tab = s.siteTable().slice();
1816
1817 int totalSize = s.numSiteTable();
1818
1819 unordered_vaxpby3_user_arg arg(x, y, d, aptr, bptr, tab);
1820
1821 dispatch_to_threads(totalSize, arg, unordered_vaxpby3_evaluate_function);
1822
1824 // Original code
1826 /*for(int j=0; j < s.numSiteTable(); j++) {
1827 int i=tab[j];
1828
1829 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1830 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1831 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1832
1833 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1834 vaxpby3(zptr, aptr, xptr, bptr, yptr, Ns);
1835
1836 }*/
1837 }
1838}
1839
1840// z = ax - by
1841template<>
1842inline
1843void evaluate( OLattice< TVec > &d,
1844 const OpAssign &op,
1845 const QDPExpr<
1846 BinaryNode<OpSubtract,
1847 BinaryNode<OpMultiply,
1848 Reference< QDPType< TScal, OScalar< TScal > > >,
1849 Reference< QDPType< TVec, OLattice< TVec > > > >,
1850 BinaryNode<OpMultiply,
1851 Reference< QDPType< TScal, OScalar< TScal > > >,
1852 Reference< QDPType< TVec, OLattice< TVec > > > > >,
1853 OLattice< TVec > > &rhs,
1854 const Subset& s)
1855{
1856
1857#ifdef DEBUG_BLAS
1858 QDPIO::cout << "z = a*x - b*y" << endl;
1859#endif
1860
1861 // Peel the stuff out of the expression
1862 // y is the right side of rhs
1863
1864 // ax is the left side of rhs and is in a binary node
1865 typedef BinaryNode<OpMultiply,
1866 Reference< QDPType< TScal, OScalar< TScal > > >,
1867 Reference< QDPType< TVec, OLattice< TVec > > > > BN;
1868
1869 // get the binary node
1870 const BN &mulNode1 = static_cast<const BN&> (rhs.expression().left());
1871 const BN &mulNode2 = static_cast<const BN&> (rhs.expression().right());
1872
1873 // get a and x out of the binary node
1874 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
1875 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
1876
1877 // get b and y out of the binary node
1878 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1879 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right());
1880
1881
1882 // Set pointers
1883 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1884 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1885
1886 if( s.hasOrderedRep() ) {
1887 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1888 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1889 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1890
1891 int total_n_3vec = (s.end()-s.start()+1);
1892
1893 ordered_vaxmby3_user_arg a = {zptr, aptr, xptr, bptr, yptr};
1894
1895 dispatch_to_threads(total_n_3vec, a, ordered_vaxmby3_evaluate_function);
1896
1898 // Original code
1900 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1901 //int n_3vec = (s.end()-s.start()+1)*Ns;
1902 //vaxmby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
1903 }
1904 else {
1905 const int* tab = s.siteTable().slice();
1906
1907 int totalSize = s.numSiteTable();
1908
1909 unordered_vaxmby3_user_arg arg(x, y, d, aptr, bptr, tab);
1910
1911 dispatch_to_threads(totalSize, arg, unordered_vaxmby3_evaluate_function);
1912
1914 // Original code
1916 /*for(int j=0; j < s.numSiteTable(); j++) {
1917 int i=tab[j];
1918
1919 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
1920 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
1921 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
1922
1923 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1924 vaxmby3(zptr, aptr, xptr, bptr, yptr, Ns);
1925
1926 }*/
1927 }
1928}
1929
1930
1931// z = xa - by
1932template<>
1933inline
1934void evaluate( OLattice< TVec > &d,
1935 const OpAssign &op,
1936 const QDPExpr<
1937 BinaryNode<OpSubtract,
1938 BinaryNode<OpMultiply,
1939 Reference< QDPType< TVec, OLattice< TVec > > >,
1940 Reference< QDPType< TScal, OScalar< TScal > > > >,
1941 BinaryNode<OpMultiply,
1942 Reference< QDPType< TScal, OScalar< TScal > > >,
1943 Reference< QDPType< TVec, OLattice< TVec > > > > >,
1944 OLattice< TVec > > &rhs,
1945 const Subset& s)
1946{
1947
1948#ifdef DEBUG_BLAS
1949 QDPIO::cout << "z = x*a - b*y" << endl;
1950#endif
1951
1952 // Peel the stuff out of the expression
1953 // y is the right side of rhs
1954
1955 // ax is the left side of rhs and is in a binary node
1956 typedef BinaryNode<OpMultiply,
1957 Reference< QDPType< TVec, OLattice< TVec > > >,
1958 Reference< QDPType< TScal, OScalar< TScal > > > > BN1;
1959
1960 typedef BinaryNode<OpMultiply,
1961 Reference< QDPType< TScal, OScalar< TScal > > >,
1962 Reference< QDPType< TVec, OLattice< TVec > > > > BN2;
1963
1964
1965
1966 // get the binary node
1967 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
1968 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
1969
1970 // get a and x out of the binary node
1971 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.left());
1972
1973 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.right());
1974
1975 // get b and y out of the binary node
1976 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.left());
1977 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.right());
1978
1979
1980 // Set pointers
1981 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
1982 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
1983
1984 if( s.hasOrderedRep() ) {
1985 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
1986 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
1987 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
1988
1989 int total_n_3vec = (s.end()-s.start()+1);
1990
1991 ordered_vaxmby3_user_arg a = {zptr, aptr, xptr, bptr, yptr};
1992
1993 dispatch_to_threads(total_n_3vec, a, ordered_vaxmby3_evaluate_function);
1994
1996 // Original code
1998 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
1999 //int n_3vec = (s.end()-s.start()+1)*Ns;
2000 //vaxmby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
2001 }
2002 else {
2003 const int* tab = s.siteTable().slice();
2004
2005 int totalSize = s.numSiteTable();
2006
2007 unordered_vaxmby3_user_arg arg(x, y, d, aptr, bptr, tab);
2008
2009 dispatch_to_threads(totalSize, arg, unordered_vaxmby3_evaluate_function);
2010
2012 // Original code
2014 /*for(int j=0; j < s.numSiteTable(); j++) {
2015 int i=tab[j];
2016
2017 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
2018 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
2019 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
2020
2021 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2022 vaxmby3(zptr, aptr, xptr, bptr, yptr, Ns);
2023
2024 }*/
2025 }
2026}
2027
2028// z = ax - yb
2029template<>
2030inline
2031void evaluate( OLattice< TVec > &d,
2032 const OpAssign &op,
2033 const QDPExpr<
2034 BinaryNode<OpSubtract,
2035 BinaryNode<OpMultiply,
2036 Reference< QDPType< TScal, OScalar< TScal > > >,
2037 Reference< QDPType< TVec, OLattice< TVec > > > >,
2038 BinaryNode<OpMultiply,
2039 Reference< QDPType< TVec, OLattice< TVec > > >,
2040 Reference< QDPType< TScal, OScalar< TScal > > > > >,
2041 OLattice< TVec > > &rhs,
2042 const Subset& s)
2043{
2044
2045#ifdef DEBUG_BLAS
2046 QDPIO::cout << "z = a*x - y*b" << endl;
2047#endif
2048
2049 // Peel the stuff out of the expression
2050 // y is the right side of rhs
2051
2052 // type of a*x
2053 typedef BinaryNode<OpMultiply,
2054 Reference< QDPType< TScal, OScalar< TScal > > >,
2055 Reference< QDPType< TVec, OLattice< TVec > > > > BN1;
2056
2057 // type of y*b
2058 typedef BinaryNode<OpMultiply,
2059 Reference< QDPType< TVec, OLattice< TVec > > >,
2060 Reference< QDPType< TScal, OScalar< TScal > > > > BN2;
2061
2062
2063
2064 // get the binary nodes
2065 // a*x node
2066 const BN1 &mulNode1 = static_cast<const BN1&> (rhs.expression().left());
2067
2068 // y*b node
2069 const BN2 &mulNode2 = static_cast<const BN2&> (rhs.expression().right());
2070
2071 // get a and x out of the binary node
2072 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.left());
2073
2074 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.right());
2075
2076
2077 // get b and y out of the binary node
2078 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.left());
2079
2080 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.right());
2081
2082
2083 // Set pointers
2084 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
2085 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
2086
2087 if( s.hasOrderedRep() ) {
2088 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
2089 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
2090 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
2091
2092 int total_n_3vec = (s.end()-s.start()+1);
2093
2094 ordered_vaxmby3_user_arg a = {zptr, aptr, xptr, bptr, yptr};
2095
2096 dispatch_to_threads(total_n_3vec, a, ordered_vaxmby3_evaluate_function);
2097
2099 // Original code
2101 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2102 //int n_3vec = (s.end()-s.start()+1)*Ns;
2103 //vaxmby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
2104 }
2105 else {
2106 const int* tab = s.siteTable().slice();
2107
2108 int totalSize = s.numSiteTable();
2109
2110 unordered_vaxmby3_user_arg arg(x, y, d, aptr, bptr, tab);
2111
2112 dispatch_to_threads(totalSize, arg, unordered_vaxmby3_evaluate_function);
2113
2115 // Original code
2117 /*for(int j=0; j < s.numSiteTable(); j++) {
2118 int i=tab[j];
2119
2120 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
2121 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
2122 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
2123
2124 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2125 vaxmby3(zptr, aptr, xptr, bptr, yptr, Ns);
2126
2127 }*/
2128 }
2129}
2130
2131// z = xa - yb
2132template<>
2133inline
2134void evaluate( OLattice< TVec > &d,
2135 const OpAssign &op,
2136 const QDPExpr<
2137 BinaryNode<OpSubtract,
2138 BinaryNode<OpMultiply,
2139 Reference< QDPType< TVec, OLattice< TVec > > >,
2140 Reference< QDPType< TScal, OScalar< TScal > > > >,
2141 BinaryNode<OpMultiply,
2142 Reference< QDPType< TVec, OLattice< TVec > > >,
2143 Reference< QDPType< TScal, OScalar< TScal > > > > >,
2144 OLattice< TVec > > &rhs,
2145 const Subset& s)
2146{
2147
2148#ifdef DEBUG_BLAS
2149 QDPIO::cout << "z = x*a - y*b" << endl;
2150#endif
2151
2152 // Peel the stuff out of the expression
2153 // y is the right side of rhs
2154
2155 // ax is the left side of rhs and is in a binary node
2156 typedef BinaryNode<OpMultiply,
2157 Reference< QDPType< TVec, OLattice< TVec > > >,
2158 Reference< QDPType< TScal, OScalar< TScal > > > > BN;
2159
2160 // get the binary node
2161 const BN &mulNode1 = static_cast<const BN&> (rhs.expression().left());
2162 const BN &mulNode2 = static_cast<const BN&> (rhs.expression().right());
2163
2164 // get a and x out of the binary node
2165 const OLattice< TVec >& x = static_cast<const OLattice< TVec >&>(mulNode1.left());
2166 const OScalar< TScal >& a = static_cast<const OScalar< TScal >&>(mulNode1.right());
2167
2168 // get b and y out of the binary node
2169 const OLattice< TVec >& y = static_cast<const OLattice< TVec >&>(mulNode2.left());
2170
2171 const OScalar< TScal >& b = static_cast<const OScalar< TScal >&>(mulNode2.right());
2172
2173 // Set pointers
2174 REAL *aptr = (REAL *)&(a.elem().elem().elem().elem());
2175 REAL *bptr = (REAL *)&(b.elem().elem().elem().elem());
2176 if( s.hasOrderedRep() ) {
2177
2178 REAL *xptr = (REAL *) &(x.elem(s.start()).elem(0).elem(0).real());
2179 REAL *yptr = (REAL *) &(y.elem(s.start()).elem(0).elem(0).real());
2180 REAL* zptr = &(d.elem(s.start()).elem(0).elem(0).real());
2181
2182 int total_n_3vec = (s.end()-s.start()+1);
2183
2184 ordered_vaxmby3_user_arg a = {zptr, aptr, xptr, bptr, yptr};
2185
2186 dispatch_to_threads(total_n_3vec, a, ordered_vaxmby3_evaluate_function);
2187
2189 // Original code
2191 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2192 //int n_3vec = (s.end()-s.start()+1)*Ns;
2193 //vaxmby3(zptr, aptr, xptr, bptr, yptr, n_3vec);
2194 }
2195 else {
2196 const int* tab = s.siteTable().slice();
2197
2198 int totalSize = s.numSiteTable();
2199
2200 unordered_vaxmby3_user_arg arg(x, y, d, aptr, bptr, tab);
2201
2202 dispatch_to_threads(totalSize, arg, unordered_vaxmby3_evaluate_function);
2203
2205 // Original code
2207 /*for(int j=0; j < s.numSiteTable(); j++) {
2208 int i=tab[j];
2209
2210 REAL *xptr = (REAL *) &(x.elem(i).elem(0).elem(0).real());
2211 REAL *yptr = (REAL *) &(y.elem(i).elem(0).elem(0).real());
2212 REAL* zptr = &(d.elem(i).elem(0).elem(0).real());
2213
2214 // Get the no of 3vecs. s.start() and s.end() are inclusive so add +1
2215 vaxmby3(zptr, aptr, xptr, bptr, yptr, Ns);
2216
2217 }*/
2218 }
2219}
2220
2221
2222// Global norm squared of a vector...
2223template<>
2224inline UnaryReturn<OLattice< TVec >, FnNorm2>::Type_t
2225norm2(const QDPType<TVec ,OLattice< TVec > >& s1, const Subset& s)
2226{
2227#ifdef DEBUG_BLAS
2228 QDPIO::cout << "Using BJ sumsq" << endl;
2229#endif
2230
2231 if ( s.hasOrderedRep() ) {
2232
2233#ifdef DEBUG_BLAS
2234 QDPIO::cout << "BJ sumsq " << endl;
2235#endif
2236 int n_3vec = (s.end() - s.start() + 1);
2237 const REAL *s1ptr = &(s1.elem(s.start()).elem(0).elem(0).real());
2238
2239 // Has Type OScalar< PScalar < PScalar < RScalar < REAL > > > >
2240
2241 DOUBLE lsum =(double)0;
2242
2243 local_sumsq(&lsum,(REAL *)s1ptr, n_3vec);
2244 UnaryReturn< OLattice< TVec >, FnNorm2>::Type_t gsum(lsum);
2245 QDPInternal::globalSum(gsum);
2246 return gsum;
2247 }
2248 else {
2249
2250 // Has Type OScalar< PScalar < PScalar < RScalar < REAL > > > >
2251 DOUBLE lsum =(DOUBLE)0;
2252 DOUBLE ltmp =(DOUBLE)0;
2253
2254 const int* tab=s.siteTable().slice();
2255 for(int j=0; j < s.numSiteTable(); j++) {
2256 int i=tab[j];
2257 REAL* s1ptr = (REAL *)&(s1.elem(i).elem(0).elem(0).real());
2258 local_sumsq(&ltmp,s1ptr,1);
2259 lsum +=ltmp;
2260 }
2261
2262 UnaryReturn< OLattice< TVec >, FnNorm2>::Type_t gsum(lsum);
2263 QDPInternal::globalSum(gsum);
2264 return gsum;
2265 }
2266}
2267
2268
2269template<>
2270inline UnaryReturn<OLattice< TVec >, FnNorm2>::Type_t
2271norm2(const QDPType<TVec ,OLattice< TVec > >& s1)
2272{
2273#ifdef DEBUG_BLAS
2274 QDPIO::cout << "Using BJ sumsq all" << endl;
2275#endif
2276
2277 int n_3vec = (all.end() - all.start() + 1);
2278 const REAL *s1ptr = &(s1.elem(all.start()).elem(0).elem(0).real());
2279
2280
2281
2282 DOUBLE lsum = 0;
2283 local_sumsq(&lsum, (REAL *)s1ptr, n_3vec);
2284 UnaryReturn< OLattice< TVec >, FnNorm2>::Type_t gsum(lsum);
2285 QDPInternal::globalSum(gsum);
2286 return gsum;
2287}
2288
2289
2290
2291template<>
2292inline BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t
2293innerProduct(const QDPType< TVec, OLattice<TVec> > &v1,
2294 const QDPType< TVec, OLattice<TVec> > &v2)
2295{
2296#ifdef DEBUG_BLAS
2297 QDPIO::cout << "BJ: innerProduct all" << endl;
2298#endif
2299
2300 // This BinaryReturn has Type_t
2301 // OScalar<OScalar<OScalar<RComplex<PScalar<REAL> > > > >
2302 BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t lprod;
2303 // Inner product is accumulated internally in DOUBLE
2304 DOUBLE ip[2];
2305 ip[0]=0;
2306 ip[1]=0;
2307
2308 // Length of subset
2309 unsigned long n_3vec = (all.end() - all.start() + 1);
2310
2311 // Call My CDOT
2312 l_vcdot(&(ip[0]), &(ip[1]),
2313 (REAL *)&(v1.elem(all.start()).elem(0).elem(0).real()),
2314 (REAL *)&(v2.elem(all.start()).elem(0).elem(0).real()),
2315 n_3vec);
2316
2317
2318 // Global sum -- still on a vector of doubles
2319 QDPInternal::globalSumArray(ip,2);
2320
2321 // Downcast (and possibly lose precision) here
2322 lprod.elem().elem().elem().real() = ip[0];
2323 lprod.elem().elem().elem().imag() = ip[1];
2324
2325 // Return
2326 return lprod;
2327}
2328
2329template<>
2330inline BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t
2331innerProduct(const QDPType< TVec, OLattice<TVec> > &v1,
2332 const QDPType< TVec, OLattice<TVec> > &v2,
2333 const Subset& s)
2334{
2335
2336 if( s.hasOrderedRep() ) {
2337#ifdef DEBUG_BLAS
2338 QDPIO::cout << "BJ: innerProduct s" << endl;
2339#endif
2340
2341 // This BinaryReturn has Type_t
2342 // OScalar<OScalar<OScalar<RComplex<PScalar<REAL> > > > >
2343 BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t lprod;
2344 DOUBLE ip[2];
2345 ip[0] = 0;
2346 ip[1] = 0;
2347
2348 unsigned long n_3vec = (s.end() - s.start() + 1);
2349 l_vcdot(&(ip[0]), &(ip[1]),
2350 (REAL *)&(v1.elem(s.start()).elem(0).elem(0).real()),
2351 (REAL *)&(v2.elem(s.start()).elem(0).elem(0).real()),
2352 n_3vec);
2353
2354
2355 QDPInternal::globalSumArray(ip,2);
2356
2357 lprod.elem().elem().elem().real() = ip[0];
2358 lprod.elem().elem().elem().imag() = ip[1];
2359
2360
2361 return lprod;
2362 }
2363 else {
2364
2365 BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t lprod;
2366 DOUBLE ip[2], ip_tmp[2];
2367 ip[0] = 0;
2368 ip[1] = 0;
2369
2370 const int *tab = s.siteTable().slice();
2371 for(int j=0; j < s.numSiteTable(); j++) {
2372
2373 int i=tab[j];
2374
2375 l_vcdot(&(ip_tmp[0]), &(ip_tmp[1]),
2376 (REAL *)&(v1.elem(i).elem(0).elem(0).real()),
2377 (REAL *)&(v2.elem(i).elem(0).elem(0).real()),
2378 1);
2379
2380 ip[0] += ip_tmp[0];
2381 ip[1] += ip_tmp[1];
2382 }
2383
2384 QDPInternal::globalSumArray(ip,2);
2385
2386 lprod.elem().elem().elem().real() = ip[0];
2387 lprod.elem().elem().elem().imag() = ip[1];
2388
2389
2390 return lprod;
2391
2392 }
2393}
2394
2395
2396
2397// Inner Product Real
2398template<>
2399inline
2400BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t
2401innerProductReal(const QDPType< TVec, OLattice<TVec> > &v1,
2402 const QDPType< TVec, OLattice<TVec> > &v2)
2403{
2404#ifdef DEBUG_BLAS
2405 QDPIO::cout << "BJ: innerProductReal all" << endl;
2406#endif
2407
2408 // This BinaryReturn has Type_t
2409 // OScalar<OScalar<OScalar<RScalar<PScalar<REAL> > > > >
2410 BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t lprod;
2411 // Inner product is accumulated internally in DOUBLE
2412 DOUBLE ip_re=0;
2413
2414 // Length of subset
2415 unsigned long n_3vec = (all.end() - all.start() + 1);
2416
2417 // Call My CDOT
2418 l_vcdot_real(&ip_re,
2419 (REAL *)&(v1.elem(all.start()).elem(0).elem(0).real()),
2420 (REAL *)&(v2.elem(all.start()).elem(0).elem(0).real()),
2421 n_3vec);
2422
2423 // Global sum
2424 QDPInternal::globalSum(ip_re);
2425
2426 // Whether CDOT did anything or not ip_re and ip_im should
2427 // now be right. Assign them to the ReturnType
2428 lprod.elem().elem().elem().elem() = ip_re;
2429
2430
2431 // Return
2432 return lprod;
2433}
2434
2435
2436template<>
2437inline
2438BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t
2439innerProductReal(const QDPType< TVec, OLattice<TVec> > &v1,
2440 const QDPType< TVec, OLattice<TVec> > &v2,
2441 const Subset& s)
2442{
2443 if( s.hasOrderedRep() ) {
2444#ifdef DEBUG_BLAS
2445 QDPIO::cout << "BJ: innerProductReal s" << endl;
2446#endif
2447
2448 // This BinaryReturn has Type_t
2449 // OScalar<OScalar<OScalar<RScalar<PScalar<REAL> > > > >
2450 BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t lprod;
2451 DOUBLE ip_re=0;
2452
2453 unsigned long n_3vec = (s.end() - s.start() + 1);
2454 l_vcdot_real(&ip_re,
2455 (REAL *)&(v1.elem(s.start()).elem(0).elem(0).real()),
2456 (REAL *)&(v2.elem(s.start()).elem(0).elem(0).real()),
2457 n_3vec);
2458
2459 QDPInternal::globalSum(ip_re);
2460 lprod.elem().elem().elem().elem() = ip_re;
2461
2462
2463 return lprod;
2464 }
2465 else {
2466
2467
2468 BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t lprod;
2469 DOUBLE ip_re=0, ip_re_tmp;
2470
2471
2472 const int *tab = s.siteTable().slice();
2473 for(int j=0; j < s.numSiteTable(); j++) {
2474
2475 int i=tab[j];
2476
2477 l_vcdot_real(&ip_re_tmp,
2478 (REAL *)&(v1.elem(i).elem(0).elem(0).real()),
2479 (REAL *)&(v2.elem(i).elem(0).elem(0).real()),
2480 1);
2481
2482 ip_re += ip_re_tmp;
2483 }
2484 QDPInternal::globalSum(ip_re);
2485 lprod.elem().elem().elem().elem() = ip_re;
2486 return lprod;
2487 }
2488}
2489
2490
2491template<>
2492inline UnaryReturn<OLattice< TVec >, FnNorm2>::Type_t
2493norm2(const multi1d< OLattice< TVec > >& s1)
2494{
2495#ifdef DEBUG_BLAS
2496 QDPIO::cout << "Using SSE multi1d sumsq all" << endl;
2497#endif
2498
2499 int n_3vec = (all.end() - all.start() + 1);
2500 DOUBLE ltmp = 0;
2501 for(int n=0; n < s1.size(); ++n)
2502 {
2503 const REAL* s1ptr = &(s1[n].elem(all.start()).elem(0).elem(0).real());
2504
2505 // I am relying on this being a Double here
2506 DOUBLE lltmp;
2507 local_sumsq(&lltmp, (REAL*)s1ptr, n_3vec);
2508
2509 ltmp += lltmp;
2510 }
2511
2512 UnaryReturn< OLattice< TVec >, FnNorm2>::Type_t lsum(ltmp);
2513 QDPInternal::globalSum(lsum);
2514 return lsum;
2515}
2516
2517
2518template<>
2519inline BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t
2520innerProduct(const multi1d< OLattice<TVec> > &v1,
2521 const multi1d< OLattice<TVec> > &v2)
2522{
2523#ifdef DEBUG_BLAS
2524 QDPIO::cout << "BJ: multi1d innerProduct all" << endl;
2525#endif
2526
2527 // This BinaryReturn has Type_t
2528 // OScalar<OScalar<OScalar<RComplex<PScalar<REAL> > > > >
2529 BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProduct>::Type_t lprod;
2530 // Inner product is accumulated internally in DOUBLE
2531 DOUBLE ip[2];
2532 ip[0]=0;
2533 ip[1]=0;
2534
2535 // Length of subset
2536 unsigned long n_3vec = (all.end() - all.start() + 1);
2537
2538 for(int n=0; n < v1.size(); ++n)
2539 {
2540 DOUBLE iip[2];
2541 iip[0]=0;
2542 iip[1]=0;
2543
2544 // Call My CDOT
2545 l_vcdot(&(iip[0]), &(iip[1]),
2546 (REAL *)&(v1[n].elem(all.start()).elem(0).elem(0).real()),
2547 (REAL *)&(v2[n].elem(all.start()).elem(0).elem(0).real()),
2548 n_3vec);
2549
2550 ip[0] += iip[0];
2551 ip[1] += iip[1];
2552 }
2553
2554 // Global sum -- still on a vector of doubles
2555 QDPInternal::globalSumArray(ip,2);
2556
2557 // Downcast (and possibly lose precision) here
2558 lprod.elem().elem().elem().real() = ip[0];
2559 lprod.elem().elem().elem().imag() = ip[1];
2560
2561 // Return
2562 return lprod;
2563}
2564
2565
2566// Inner Product Real
2567template<>
2568inline
2569BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t
2570innerProductReal(const multi1d< OLattice<TVec> > &v1,
2571 const multi1d< OLattice<TVec> > &v2)
2572{
2573#ifdef DEBUG_BLAS
2574 QDPIO::cout << "BJ: innerProductReal(multi1d) all" << endl;
2575#endif
2576
2577 // This BinaryReturn has Type_t
2578 // OScalar<OScalar<OScalar<RScalar<PScalar<REAL> > > > >
2579 BinaryReturn< OLattice<TVec>, OLattice<TVec>, FnInnerProductReal>::Type_t lprod;
2580 // Inner product is accumulated internally in DOUBLE
2581 DOUBLE ip_re=0;
2582
2583 // Length of subset
2584 unsigned long n_3vec = (all.end() - all.start() + 1);
2585
2586 for(int n=0; n < v1.size(); ++n)
2587 {
2588 DOUBLE iip_re=0;
2589
2590 // Call My CDOT
2591 l_vcdot_real(&iip_re,
2592 (REAL *)&(v1[n].elem(all.start()).elem(0).elem(0).real()),
2593 (REAL *)&(v2[n].elem(all.start()).elem(0).elem(0).real()),
2594 n_3vec);
2595
2596 ip_re += iip_re;
2597 }
2598
2599 // Global sum
2600 QDPInternal::globalSum(ip_re);
2601
2602 // Whether CDOT did anything or not ip_re and ip_im should
2603 // now be right. Assign them to the ReturnType
2604 lprod.elem().elem().elem().elem() = ip_re;
2605
2606
2607 // Return
2608 return lprod;
2609}
2610
2611
2612} // namespace QDP;
2613
2614#endif // guard
2615
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
Container for a multi-dimensional 1D array.
Definition qdp_multi.h:25
QDP::multi1d::slice
const T * slice() const
Return ref to a column slice.
Definition qdp_multi.h:225
generic_blas_local_sumsq.h
Generic Scalar, local sum squared routine.
generic_blas_local_vcdot.h
Generic Scalar, CDOT routine.
generic_blas_local_vcdot_real.h
Generic Scalar, CDOT routine.
generic_blas_vaxmby3.h
Generic Scalar VAXMBY routine.
generic_blas_vaxmby3_norm.h
Generic Scalar VAXPY NORM routine.
generic_blas_vaxmy3.h
Generic Scalar VAXMY routine.
generic_blas_vaxpby3.h
Generic Scalar VAXPBY routine.
generic_blas_vaxpy3.h
Generic Scalar VAXPY routine.
generic_blas_vaxpy3_norm.h
Generic Scalar VAXPY NORM routine.
generic_blas_vscal.h
Generic Scalar VSCAL routine.
QDP::REAL
REAL32 REAL
Definition qdp_scalarsite_defs.h:28
QDP::DOUBLE
REAL64 DOUBLE
Definition qdp_scalarsite_defs.h:29
QDP::innerProductReal
BinaryReturn< C1, C2, FnInnerProductReal >::Type_t innerProductReal(const QDPType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
OScalar = innerProductReal(adj(source1)*source2).
Definition qdp_globalfuncs.h:172
QDP::norm2
UnaryReturn< C, FnNorm2 >::Type_t norm2(const QDPType< T, C > &s1)
OScalar = norm2(trace(adj(source)*source)).
Definition qdp_globalfuncs.h:55
QDP::innerProduct
BinaryReturn< C1, C2, FnInnerProduct >::Type_t innerProduct(const QDPType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
OScalar = innerProduct(adj(source1)*source2).
Definition qdp_globalfuncs.h:98
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::all
Subset all
Default all subset.
Definition qdp_subset.cc:16
QDP::QDPIO::cout
StandardOutputStream cout
Definition qdp_stdio.cc:21
QDP::QDPInternal::globalSum
void globalSum(T &dest)
Sum across all nodes.
Definition qdp_parscalar_specific.h:125
QDP::QDPInternal::globalSumArray
void globalSumArray(unsigned int *dest, int len)
Wrapper to get a functional unsigned global sum.
Definition qdp_parscalar_specific.h:65
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::unordered_vaxpy3_z_evaluate_function
void unordered_vaxpy3_z_evaluate_function(int lo, int hi, int myId, unordered_vaxpy3_z_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:104
QDP::ordered_vaxmy3_evaluate_function
void ordered_vaxmy3_evaluate_function(int lo, int hi, int myId, ordered_vaxmy3_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:137
QDP::dispatch_to_threads
void dispatch_to_threads(int numSiteTable, Arg a, void(*func)(int, int, int, Arg *))
Definition qdp_dispatch.h:105
QDP::ordered_vaxpby3_evaluate_function
void ordered_vaxpby3_evaluate_function(int lo, int hi, int myId, ordered_vaxpby3_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:297
QDP::unordered_vscal_evaluate_function
void unordered_vscal_evaluate_function(int lo, int hi, int myId, unordered_vscal_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:264
QDP::l_vcdot
void l_vcdot(DOUBLE *Out_re, DOUBLE *Out_im, REAL *V1, REAL *V2, int n_3vec)
Definition generic_blas_local_vcdot.h:20
QDP::ordered_vaxmby3_evaluate_function
void ordered_vaxmby3_evaluate_function(int lo, int hi, int myId, ordered_vaxmby3_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:370
QDP::ordered_vaxpy3_evaluate_function
void ordered_vaxpy3_evaluate_function(int lo, int hi, int myId, ordered_vaxpy3_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:25
QDP::unordered_vaxmby3_evaluate_function
void unordered_vaxmby3_evaluate_function(int lo, int hi, int myId, unordered_vaxmby3_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:405
QDP::TVec
PSpinVector< PColorVector< RComplex< REAL >, 3 >, Ns > TVec
Definition qdp_scalarsite_bagel_qdp_blas.h:19
QDP::unordered_vaxmy3_z_evaluate_function
void unordered_vaxmy3_z_evaluate_function(int lo, int hi, int myId, unordered_vaxmy3_z_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:202
QDP::l_vcdot_real
void l_vcdot_real(DOUBLE *Out, REAL *V1, REAL *V2, int n_3vec)
Definition generic_blas_local_vcdot_real.h:19
QDP::unordered_vaxmy3_y_evaluate_function
void unordered_vaxmy3_y_evaluate_function(int lo, int hi, int myId, unordered_vaxmy3_y_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:169
QDP::local_sumsq
void local_sumsq(DOUBLE *Out, REAL *In, int n_3vec)
Definition generic_blas_local_sumsq.h:16
QDP::ordered_vscal_evaluate_function
void ordered_vscal_evaluate_function(int lo, int hi, int myId, ordered_vscal_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:234
QDP::unordered_vaxpy3_y_evaluate_function
void unordered_vaxpy3_y_evaluate_function(int lo, int hi, int myId, unordered_vaxpy3_y_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:59
QDP::unordered_vaxpby3_evaluate_function
void unordered_vaxpby3_evaluate_function(int lo, int hi, int myId, unordered_vaxpby3_user_arg *a)
Definition qdp_scalarsite_generic_blas.h:331
qdp_scalarsite_generic_blas_wrapper.h
ordered_vaxpy3_evaluate_function
void ordered_vaxpy3_evaluate_function(int lo, int hi, int myId, ordered_vaxpy3_user_arg *a)
Definition qdp_scalarsite_generic_blas_wrapper.h:24
QDP::ordered_vaxmby3_user_arg
Definition qdp_scalarsite_generic_blas.h:360
QDP::ordered_vaxmy3_user_arg
Definition qdp_scalarsite_generic_blas.h:128
QDP::ordered_vaxpby3_user_arg
Definition qdp_scalarsite_generic_blas.h:287
QDP::ordered_vaxpy3_user_arg
Definition qdp_scalarsite_generic_blas.h:16
QDP::ordered_vscal_user_arg
Definition qdp_scalarsite_generic_blas.h:226
QDP::unordered_vaxmby3_user_arg
Definition qdp_scalarsite_generic_blas.h:387
QDP::unordered_vaxmby3_user_arg::tab
const int * tab
Definition qdp_scalarsite_generic_blas.h:400
QDP::unordered_vaxmby3_user_arg::aptr
REAL * aptr
Definition qdp_scalarsite_generic_blas.h:398
QDP::unordered_vaxmby3_user_arg::d
OLattice< TVec > & d
Definition qdp_scalarsite_generic_blas.h:397
QDP::unordered_vaxmby3_user_arg::y
const OLattice< TVec > & y
Definition qdp_scalarsite_generic_blas.h:396
QDP::unordered_vaxmby3_user_arg::x
const OLattice< TVec > & x
Definition qdp_scalarsite_generic_blas.h:395
QDP::unordered_vaxmy3_y_user_arg
Definition qdp_scalarsite_generic_blas.h:155
QDP::unordered_vaxmy3_y_user_arg::tab
const int * tab
Definition qdp_scalarsite_generic_blas.h:164
QDP::unordered_vaxmy3_y_user_arg::d
OLattice< TVec > & d
Definition qdp_scalarsite_generic_blas.h:162
QDP::unordered_vaxmy3_y_user_arg::x
const OLattice< TVec > & x
Definition qdp_scalarsite_generic_blas.h:161
QDP::unordered_vaxmy3_z_user_arg
Definition qdp_scalarsite_generic_blas.h:185
QDP::unordered_vaxmy3_z_user_arg::y
const OLattice< TVec > & y
Definition qdp_scalarsite_generic_blas.h:194
QDP::unordered_vaxpby3_user_arg
Definition qdp_scalarsite_generic_blas.h:314
QDP::unordered_vaxpby3_user_arg::d
OLattice< TVec > & d
Definition qdp_scalarsite_generic_blas.h:323
QDP::unordered_vaxpby3_user_arg::x
const OLattice< TVec > & x
Definition qdp_scalarsite_generic_blas.h:321
QDP::unordered_vaxpby3_user_arg::y
const OLattice< TVec > & y
Definition qdp_scalarsite_generic_blas.h:322
QDP::unordered_vaxpby3_user_arg::aptr
REAL * aptr
Definition qdp_scalarsite_generic_blas.h:324
QDP::unordered_vaxpby3_user_arg::tab
const int * tab
Definition qdp_scalarsite_generic_blas.h:326
QDP::unordered_vaxpy3_y_user_arg
Definition qdp_scalarsite_generic_blas.h:43
QDP::unordered_vaxpy3_y_user_arg::d
OLattice< TVec > & d
Definition qdp_scalarsite_generic_blas.h:51
QDP::unordered_vaxpy3_z_user_arg
Definition qdp_scalarsite_generic_blas.h:88
QDP::unordered_vaxpy3_z_user_arg::d
OLattice< TVec > & d
Definition qdp_scalarsite_generic_blas.h:97
QDP::unordered_vaxpy3_z_user_arg::x
const OLattice< TVec > & x
Definition qdp_scalarsite_generic_blas.h:95
QDP::unordered_vaxpy3_z_user_arg::tab
const int * tab
Definition qdp_scalarsite_generic_blas.h:99
QDP::unordered_vscal_user_arg
Definition qdp_scalarsite_generic_blas.h:248
QDP::unordered_vscal_user_arg::d
OLattice< TVec > & d
Definition qdp_scalarsite_generic_blas.h:257
Generated by doxygen 1.16.1