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