QDP++
qdp_scalarsite_bagel_qdp_linalg.h
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id: qdp_scalarsite_bagel_qdp_linalg.h,v 1.5 2008-05-09 17:45:20 bjoo Exp $
3
9
10#ifndef QDP_SCALARSITE_BAGEL_QDP_LINALG_H
11#define QDP_SCALARSITE_BAGEL_QDP_LINALG_H
12
13namespace QDP {
14
21
22// Use this def just to safe some typing later on in the file
23
24
25
26#include "bagel_qdp.h"
27
28#if 1
29typedef RComplex<BAGELQDPFloat> RComplexFloat;
30
31
32template<>
33inline
34void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
35 const OpAssign& op,
36 const QDPExpr<BinaryNode<OpMultiply,
37 Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
38 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > >,
39 Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
40 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
41 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > >& rhs,
42 const Subset& s) {
43
44 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
45 const C& l = static_cast<const C&>(rhs.expression().left());
46 const C& r = static_cast<const C&>(rhs.expression().right());
47
48#ifdef DEBUG_BAGELQDP_LINALG
49 QDPIO::cout << "evaluate(M*M) subset = s " << endl;
50#endif
51
52 if( s.hasOrderedRep() ) {
53 // Do whole subset
54 unsigned int start = s.start();
55 unsigned int end = s.end();
56
57 unsigned long num_sites = end - start + 1;
58
59 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
60 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
61 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
62
63
64 qdp_su3_mm(resptr, lptr, rptr, num_sites, (unsigned long)0);
65
66 }
67 else {
68 // Do site by site
69 const int* tab = s.siteTable().slice();
70 unsigned int num_sites = s.numSiteTable();
71 for(unsigned int j=0; j < num_sites; j++) {
72 int i = tab[j];
73 unsigned long one_site = 1;
74 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
75
76 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
77 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
78
79 qdp_su3_mm(resptr, lptr, rptr, one_site, (unsigned long)0);
80
81 }
82 }
83
84}
85
86template<>
87inline
88void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
89 const OpAssign& op,
90 const QDPExpr<BinaryNode<OpAdjMultiply,
91 UnaryNode<OpIdentity, Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
92 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
93 Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
94 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
95 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > >& rhs,
96 const Subset& s) {
97
98 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
99 const C& l = static_cast<const C&>(rhs.expression().left().child());
100 const C& r = static_cast<const C&>(rhs.expression().right());
101
102#ifdef DEBUG_BAGELQDP_LINALG
103 QDPIO::cout << "evaluate(A*M) subset = s " << endl;
104#endif
105
106 if( s.hasOrderedRep() ) {
107 // Do whole subset
108 unsigned int start = s.start();
109 unsigned int end = s.end();
110
111 unsigned long num_sites = end - start + 1;
112
113 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
114 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
115 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
116
117
118 qdp_su3_am(resptr, lptr, rptr, num_sites, (unsigned long)0);
119
120 }
121 else {
122 // Do site by site
123 const int* tab = s.siteTable().slice();
124 unsigned int num_sites = s.numSiteTable();
125 for(unsigned int j=0; j < num_sites; j++) {
126 int i = tab[j];
127 unsigned long one_site = 1;
128 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
129
130 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
131 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
132
133 qdp_su3_am(resptr, lptr, rptr, one_site, (unsigned long)0);
134
135 }
136 }
137
138}
139
140template<>
141inline
142void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
143 const OpAssign& op,
144 const QDPExpr<BinaryNode<OpMultiplyAdj,
145 Reference<
146 QDPType<
147 PScalar<PColorMatrix<RComplexFloat, 3> >,
148 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
149 >
150 >,
151
152 UnaryNode< OpIdentity,
153 Reference<
154 QDPType<
155 PScalar<PColorMatrix<RComplexFloat, 3> >,
156 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
157 >
158 >
159 >
160
161 >,
162
163 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
164 >& rhs,
165 const Subset& s) {
166
167 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
168 const C& l = static_cast<const C&>(rhs.expression().left());
169 const C& r = static_cast<const C&>(rhs.expression().right().child());
170
171#ifdef DEBUG_BAGELQDP_LINALG
172 QDPIO::cout << "evaluate(M*A) subset = s " << endl;
173#endif
174
175 if( s.hasOrderedRep() ) {
176 // Do whole subset
177 unsigned int start = s.start();
178 unsigned int end = s.end();
179
180 unsigned long num_sites = end - start + 1;
181
182 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
183 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
184 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
185
186
187 qdp_su3_ma(resptr, lptr, rptr, num_sites, (unsigned long)0);
188
189 }
190 else {
191 // Do site by site
192 const int* tab = s.siteTable().slice();
193 unsigned int num_sites = s.numSiteTable();
194 for(unsigned int j=0; j < num_sites; j++) {
195 int i = tab[j];
196 unsigned long one_site = 1;
197 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
198
199 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
200 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
201
202 qdp_su3_ma(resptr, lptr, rptr, one_site, (unsigned long)0);
203
204 }
205 }
206
207}
208
209
210template<>
211inline
212void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
213 const OpAssign& op,
214 const QDPExpr<BinaryNode<OpAdjMultiplyAdj,
215 UnaryNode<OpIdentity, Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
216 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
217 UnaryNode<OpIdentity, Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
218 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > > >,
219 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > >& rhs,
220 const Subset& s) {
221
222 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
223 const C& l = static_cast<const C&>(rhs.expression().left().child());
224 const C& r = static_cast<const C&>(rhs.expression().right().child());
225
226#ifdef DEBUG_BAGELQDP_LINALG
227 QDPIO::cout << "evaluate(A*A) subset = s " << endl;
228#endif
229 BAGELQDPFloat one_minus_i[2] QDP_ALIGN16;
230 one_minus_i[0] = (BAGELQDPFloat)1;
231 one_minus_i[1] = (BAGELQDPFloat)(-1);
232
233 if( s.hasOrderedRep() ) {
234 // Do whole subset
235 unsigned int start = s.start();
236 unsigned int end = s.end();
237
238 unsigned long num_sites = end - start + 1;
239
240 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
241 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
242 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
243
244
245 qdp_su3_aa(resptr, lptr, rptr, num_sites, (unsigned long)one_minus_i);
246
247 }
248 else {
249 // Do site by site
250 const int* tab = s.siteTable().slice();
251 unsigned int num_sites = s.numSiteTable();
252 for(unsigned int j=0; j < num_sites; j++) {
253 int i = tab[j];
254 unsigned long one_site = 1;
255 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
256
257 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
258 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
259
260 qdp_su3_aa(resptr, lptr, rptr, one_site, (unsigned long)one_minus_i);
261
262 }
263 }
264
265}
266
267
268 // +=
269template<>
270inline
271void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
272 const OpAddAssign& op,
273 const QDPExpr<BinaryNode<OpMultiply,
274 Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
275 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > >,
276 Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
277 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
278 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > >& rhs,
279 const Subset& s) {
280
281 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
282 const C& l = static_cast<const C&>(rhs.expression().left());
283 const C& r = static_cast<const C&>(rhs.expression().right());
284
285#ifdef DEBUG_BAGELQDP_LINALG
286 QDPIO::cout << "evaluate(M*M) subset = s " << endl;
287#endif
288
289 BAGELQDPFloat plus_one[2] QDP_ALIGN16;
290 plus_one[0] = (BAGELQDPFloat)1;
291 plus_one[1] = (BAGELQDPFloat)0;
292
293 if( s.hasOrderedRep() ) {
294 // Do whole subset
295 unsigned int start = s.start();
296 unsigned int end = s.end();
297
298 unsigned long num_sites = end - start + 1;
299
300 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
301 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
302 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
303
304 qdp_su3_mm_peq(resptr, plus_one, lptr, rptr, num_sites, (unsigned long)0);
305
306 }
307 else {
308 // Do site by site
309 const int* tab = s.siteTable().slice();
310 unsigned int num_sites = s.numSiteTable();
311 for(unsigned int j=0; j < num_sites; j++) {
312 int i = tab[j];
313 unsigned long one_site = 1;
314 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
315
316 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
317 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
318
319 qdp_su3_mm_peq(resptr, plus_one, lptr, rptr, one_site, (unsigned long)0);
320
321 }
322 }
323
324}
325
326template<>
327inline
328void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
329 const OpAddAssign& op,
330 const QDPExpr<BinaryNode<OpAdjMultiply,
331 UnaryNode<OpIdentity, Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
332 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
333 Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
334 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
335 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > >& rhs,
336 const Subset& s) {
337
338 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
339 const C& l = static_cast<const C&>(rhs.expression().left().child());
340 const C& r = static_cast<const C&>(rhs.expression().right());
341
342#ifdef DEBUG_BAGELQDP_LINALG
343 QDPIO::cout << "evaluate(A*M) subset = s " << endl;
344#endif
345
346 BAGELQDPFloat plus_one[2] QDP_ALIGN16;
347 plus_one[0] = (BAGELQDPFloat)1;
348 plus_one[1] = (BAGELQDPFloat)0;
349
350 if( s.hasOrderedRep() ) {
351 // Do whole subset
352 unsigned int start = s.start();
353 unsigned int end = s.end();
354
355 unsigned long num_sites = end - start + 1;
356
357 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
358 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
359 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
360
361
362 qdp_su3_am_peq(resptr, plus_one, lptr, rptr, num_sites, (unsigned long)0);
363
364 }
365 else {
366 // Do site by site
367 const int* tab = s.siteTable().slice();
368 unsigned int num_sites = s.numSiteTable();
369 for(unsigned int j=0; j < num_sites; j++) {
370 int i = tab[j];
371 unsigned long one_site = 1;
372 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
373
374 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
375 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
376
377 qdp_su3_am_peq(resptr, plus_one, lptr, rptr, one_site, (unsigned long)0);
378
379 }
380 }
381
382}
383
384template<>
385inline
386void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
387 const OpAddAssign& op,
388 const QDPExpr<BinaryNode<OpMultiplyAdj,
389 Reference< QDPType<
390 PScalar<PColorMatrix<RComplexFloat, 3> >,
391 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
392 >
393 >,
394
395 UnaryNode< OpIdentity,
396 Reference< QDPType<
397 PScalar<PColorMatrix<RComplexFloat, 3> >,
398 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
399 >
400 > > >,
401 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > >& rhs,
402 const Subset& s) {
403
404 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
405 const C& l = static_cast<const C&>(rhs.expression().left());
406 const C& r = static_cast<const C&>(rhs.expression().right().child());
407
408#ifdef DEBUG_BAGELQDP_LINALG
409 QDPIO::cout << "evaluate(M*A) subset = s " << endl;
410#endif
411 BAGELQDPFloat plus_one[2] QDP_ALIGN16;
412 plus_one[0] = (BAGELQDPFloat)1;
413 plus_one[1] = (BAGELQDPFloat)0;
414
415 if( s.hasOrderedRep() ) {
416 // Do whole subset
417 unsigned int start = s.start();
418 unsigned int end = s.end();
419
420 unsigned long num_sites = end - start + 1;
421
422 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
423 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
424 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
425
426
427 qdp_su3_ma_peq(resptr, plus_one, lptr, rptr, num_sites, (unsigned long)0);
428
429 }
430 else {
431 // Do site by site
432 const int* tab = s.siteTable().slice();
433 unsigned int num_sites = s.numSiteTable();
434 for(unsigned int j=0; j < num_sites; j++) {
435 int i = tab[j];
436 unsigned long one_site = 1;
437 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
438
439 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
440 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
441
442 qdp_su3_ma_peq(resptr, plus_one, lptr, rptr, one_site, (unsigned long)0);
443
444 }
445 }
446
447}
448
449
450template<>
451inline
452void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
453 const OpAddAssign& op,
454 const QDPExpr<BinaryNode<OpAdjMultiplyAdj,
455 UnaryNode<OpIdentity, Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
456 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
457 UnaryNode<OpIdentity, Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
458 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > > >,
459 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > >& rhs,
460 const Subset& s) {
461
462 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
463 const C& l = static_cast<const C&>(rhs.expression().left().child());
464 const C& r = static_cast<const C&>(rhs.expression().right().child());
465
466#ifdef DEBUG_BAGELQDP_LINALG
467 QDPIO::cout << "evaluate(A*A) subset = s " << endl;
468#endif
469 BAGELQDPFloat one_minus_i[2] QDP_ALIGN16;
470 one_minus_i[0] = (BAGELQDPFloat)1;
471 one_minus_i[1] = (BAGELQDPFloat)(-1);
472
473 BAGELQDPFloat plus_one[2] QDP_ALIGN16;
474 plus_one[0] = (BAGELQDPFloat)1;
475 plus_one[1] = (BAGELQDPFloat)0;
476
477 if( s.hasOrderedRep() ) {
478 // Do whole subset
479 unsigned int start = s.start();
480 unsigned int end = s.end();
481
482 unsigned long num_sites = end - start + 1;
483
484 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
485 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
486 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
487
488
489 qdp_su3_aa_peq(resptr, plus_one, lptr, rptr, num_sites, (unsigned long)one_minus_i);
490
491 }
492 else {
493 // Do site by site
494 const int* tab = s.siteTable().slice();
495 unsigned int num_sites = s.numSiteTable();
496 for(unsigned int j=0; j < num_sites; j++) {
497 int i = tab[j];
498 unsigned long one_site = 1;
499 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
500
501 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
502 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
503
504 qdp_su3_aa_peq(resptr, plus_one, lptr, rptr, one_site, (unsigned long)one_minus_i);
505
506 }
507 }
508
509}
510
511
512
513 // -=
514
515template<>
516inline
517void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
518 const OpSubtractAssign& op,
519 const QDPExpr<BinaryNode<OpMultiply,
520 Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
521 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > >,
522 Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
523 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
524 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > >& rhs,
525 const Subset& s) {
526
527 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
528 const C& l = static_cast<const C&>(rhs.expression().left());
529 const C& r = static_cast<const C&>(rhs.expression().right());
530
531#ifdef DEBUG_BAGELQDP_LINALG
532 QDPIO::cout << "evaluate(M*M) subset = s " << endl;
533#endif
534
535 BAGELQDPFloat minus_one[2] QDP_ALIGN16;
536 minus_one[0] = (BAGELQDPFloat)-1;
537 minus_one[1] = (BAGELQDPFloat)0;
538
539 if( s.hasOrderedRep() ) {
540 // Do whole subset
541 unsigned int start = s.start();
542 unsigned int end = s.end();
543
544 unsigned long num_sites = end - start + 1;
545
546 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
547 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
548 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
549
550 qdp_su3_mm_peq(resptr, minus_one, lptr, rptr, num_sites, (unsigned long)0);
551
552 }
553 else {
554 // Do site by site
555 const int* tab = s.siteTable().slice();
556 unsigned int num_sites = s.numSiteTable();
557 for(unsigned int j=0; j < num_sites; j++) {
558 int i = tab[j];
559 unsigned long one_site = 1;
560 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
561
562 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
563 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
564
565 qdp_su3_mm_peq(resptr, minus_one, lptr, rptr, one_site, (unsigned long)0);
566
567 }
568 }
569
570}
571
572template<>
573inline
574void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
575 const OpSubtractAssign& op,
576 const QDPExpr<BinaryNode<OpAdjMultiply,
577 UnaryNode<OpIdentity, Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
578 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
579 Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
580 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
581 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > >& rhs,
582 const Subset& s) {
583
584 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
585 const C& l = static_cast<const C&>(rhs.expression().left().child());
586 const C& r = static_cast<const C&>(rhs.expression().right());
587
588#ifdef DEBUG_BAGELQDP_LINALG
589 QDPIO::cout << "evaluate(A*M) subset = s " << endl;
590#endif
591
592 BAGELQDPFloat minus_one[2] QDP_ALIGN16;
593 minus_one[0] = (BAGELQDPFloat)-1;
594 minus_one[1] = (BAGELQDPFloat)0;
595
596 if( s.hasOrderedRep() ) {
597 // Do whole subset
598 unsigned int start = s.start();
599 unsigned int end = s.end();
600
601 unsigned long num_sites = end - start + 1;
602
603 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
604 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
605 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
606
607
608 qdp_su3_am_peq(resptr, minus_one, lptr, rptr, num_sites, (unsigned long)0);
609
610 }
611 else {
612 // Do site by site
613 const int* tab = s.siteTable().slice();
614 unsigned int num_sites = s.numSiteTable();
615 for(unsigned int j=0; j < num_sites; j++) {
616 int i = tab[j];
617 unsigned long one_site = 1;
618 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
619
620 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
621 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
622
623 qdp_su3_am_peq(resptr, minus_one, lptr, rptr, one_site, (unsigned long)0);
624
625 }
626 }
627
628}
629
630template<>
631inline
632void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
633 const OpSubtractAssign& op,
634 const QDPExpr<BinaryNode<OpMultiplyAdj,
635 Reference<
636 QDPType<
637 PScalar<PColorMatrix<RComplexFloat, 3> >,
638 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
639 >
640 >,
641
642 UnaryNode< OpIdentity,
643 Reference<
644 QDPType<
645 PScalar<PColorMatrix<RComplexFloat, 3> >,
646 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
647 >
648 >
649 >
650
651 >,
652
653 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
654 >& rhs,
655 const Subset& s) {
656
657 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
658 const C& l = static_cast<const C&>(rhs.expression().left());
659 const C& r = static_cast<const C&>(rhs.expression().right().child());
660
661#ifdef DEBUG_BAGELQDP_LINALG
662 QDPIO::cout << "evaluate(M*A) subset = s " << endl;
663#endif
664 BAGELQDPFloat minus_one[2] QDP_ALIGN16;
665 minus_one[0] = (BAGELQDPFloat)-1;
666 minus_one[1] = (BAGELQDPFloat)0;
667
668 if( s.hasOrderedRep() ) {
669 // Do whole subset
670 unsigned int start = s.start();
671 unsigned int end = s.end();
672
673 unsigned long num_sites = end - start + 1;
674
675 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
676 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
677 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
678
679
680 qdp_su3_ma_peq(resptr, minus_one, lptr, rptr, num_sites, (unsigned long)0);
681
682 }
683 else {
684 // Do site by site
685 const int* tab = s.siteTable().slice();
686 unsigned int num_sites = s.numSiteTable();
687 for(unsigned int j=0; j < num_sites; j++) {
688 int i = tab[j];
689 unsigned long one_site = 1;
690 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
691
692 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
693 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
694
695 qdp_su3_ma_peq(resptr, minus_one, lptr, rptr, one_site, (unsigned long)0);
696
697 }
698 }
699
700}
701
702
703template<>
704inline
705void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
706 const OpSubtractAssign& op,
707 const QDPExpr<BinaryNode<OpAdjMultiplyAdj,
708 UnaryNode<OpIdentity, Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
709 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > >,
710 UnaryNode<OpIdentity, Reference<QDPType<PScalar<PColorMatrix<RComplexFloat, 3> >,
711 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > > > > >,
712 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > >& rhs,
713 const Subset& s) {
714
715 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
716 const C& l = static_cast<const C&>(rhs.expression().left().child());
717 const C& r = static_cast<const C&>(rhs.expression().right().child());
718
719#ifdef DEBUG_BAGELQDP_LINALG
720 QDPIO::cout << "evaluate(A*A) subset = s " << endl;
721#endif
722 BAGELQDPFloat one_minus_i[2] QDP_ALIGN16;
723 one_minus_i[0] = (BAGELQDPFloat)1;
724 one_minus_i[1] = (BAGELQDPFloat)(-1);
725
726 BAGELQDPFloat minus_one[2] QDP_ALIGN16;
727 minus_one[0] = (BAGELQDPFloat)-1;
728 minus_one[1] = (BAGELQDPFloat)0;
729
730 if( s.hasOrderedRep() ) {
731 // Do whole subset
732 unsigned int start = s.start();
733 unsigned int end = s.end();
734
735 unsigned long num_sites = end - start + 1;
736
737 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
738 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
739 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
740
741
742 qdp_su3_aa_peq(resptr, minus_one, lptr, rptr, num_sites, (unsigned long)one_minus_i);
743
744 }
745 else {
746 // Do site by site
747 const int* tab = s.siteTable().slice();
748 unsigned int num_sites = s.numSiteTable();
749 for(unsigned int j=0; j < num_sites; j++) {
750 int i = tab[j];
751 unsigned long one_site = 1;
752 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
753
754 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
755 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
756
757 qdp_su3_aa_peq(resptr, minus_one, lptr, rptr, one_site, (unsigned long)one_minus_i);
758
759 }
760 }
761
762}
763
764
765
766
767 // += *
768template<>
769inline
770void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
771 const OpAddAssign& op,
772
773 const QDPExpr<
774 BinaryNode< OpMultiply,
775 BinaryNode<OpMultiply,
776 Reference<
777 QDPType<
778 PScalar<PScalar<RScalar<BAGELQDPFloat> > >,
779 OScalar<PScalar<PScalar<RScalar<BAGELQDPFloat> > > >
780 >
781 >,
782 Reference<
783 QDPType<
784 PScalar<PColorMatrix<RComplexFloat, 3> >,
785 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
786 >
787 >
788 >,
789 Reference<
790 QDPType<
791 PScalar<PColorMatrix<RComplexFloat, 3> >,
792 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
793 >
794 >
795 >,
796 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
797 >& rhs,
798
799 const Subset& s) {
800
801 typedef OScalar<PScalar<PScalar<RScalar<BAGELQDPFloat> > > > F;
802
803 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
804
805 typedef BinaryNode<OpMultiply,
806 Reference<
807 QDPType<
808 PScalar<PScalar<RScalar<BAGELQDPFloat> > >,
809 OScalar<PScalar<PScalar<RScalar<BAGELQDPFloat> > > >
810 >
811 >,
812 Reference<
813 QDPType<
814 PScalar<PColorMatrix<RComplexFloat, 3> >,
815 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
816 >
817 >
818 > BN;
819
820
821 const BN& node = static_cast<const BN&>(rhs.expression().left());
822 const F& scal = static_cast<const F&>(node.left());
823 const C& l = static_cast<const C&>(node.right());
824 const C& r = static_cast<const C&>(rhs.expression().right());
825
826#if DEBUG_BAGELQDP_LINALG
827 QDPIO::cout << "evaluate(M += alpha*M*M 2 ) subset = s " << endl;
828#endif
829
830 BAGELQDPFloat scalar[2] QDP_ALIGN16;
831 scalar[0] = (BAGELQDPFloat)(scal.elem().elem().elem().elem());
832 scalar[1] = (BAGELQDPFloat)0;
833
834 if( s.hasOrderedRep() ) {
835 // Do whole subset
836 unsigned int start = s.start();
837 unsigned int end = s.end();
838
839 unsigned long num_sites = end - start + 1;
840
841 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
842 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
843 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
844
845 qdp_su3_mm_peq(resptr, scalar, lptr, rptr, num_sites, (unsigned long)0);
846
847 }
848 else {
849 // Do site by site
850 const int* tab = s.siteTable().slice();
851 unsigned int num_sites = s.numSiteTable();
852 for(unsigned int j=0; j < num_sites; j++) {
853 int i = tab[j];
854 unsigned long one_site = 1;
855 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
856
857 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
858 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
859
860 qdp_su3_mm_peq(resptr, scalar, lptr, rptr, one_site, (unsigned long)0);
861
862 }
863 }
864
865}
866
867 // +-= *
868template<>
869inline
870void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
871 const OpSubtractAssign& op,
872
873 const QDPExpr<
874 BinaryNode< OpMultiply,
875 BinaryNode<OpMultiply,
876 Reference<
877 QDPType<
878 PScalar<PScalar<RScalar<BAGELQDPFloat> > >,
879 OScalar<PScalar<PScalar<RScalar<BAGELQDPFloat> > > >
880 >
881 >,
882 Reference<
883 QDPType<
884 PScalar<PColorMatrix<RComplexFloat, 3> >,
885 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
886 >
887 >
888 >,
889 Reference<
890 QDPType<
891 PScalar<PColorMatrix<RComplexFloat, 3> >,
892 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
893 >
894 >
895 >,
896 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
897 >& rhs,
898
899 const Subset& s) {
900
901 typedef OScalar<PScalar<PScalar<RScalar<BAGELQDPFloat> > > > F;
902
903 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
904
905 typedef BinaryNode<OpMultiply,
906 Reference<
907 QDPType<
908 PScalar<PScalar<RScalar<BAGELQDPFloat> > >,
909 OScalar<PScalar<PScalar<RScalar<BAGELQDPFloat> > > >
910 >
911 >,
912 Reference<
913 QDPType<
914 PScalar<PColorMatrix<RComplexFloat, 3> >,
915 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
916 >
917 >
918 > BN;
919
920
921 const BN& node = static_cast<const BN&>(rhs.expression().left());
922 const F& scal = static_cast<const F&>(node.left());
923 const C& l = static_cast<const C&>(node.right());
924 const C& r = static_cast<const C&>(rhs.expression().right());
925
926#if DEBUG_BAGELQDP_LINALG
927 QDPIO::cout << "evaluate(M -= alpha*M*M 2 ) subset = s " << endl;
928#endif
929
930 BAGELQDPFloat scalar[2] QDP_ALIGN16;
931 scalar[0] = -(BAGELQDPFloat)(scal.elem().elem().elem().elem());
932 scalar[1] = (BAGELQDPFloat)0;
933
934 if( s.hasOrderedRep() ) {
935 // Do whole subset
936 unsigned int start = s.start();
937 unsigned int end = s.end();
938
939 unsigned long num_sites = end - start + 1;
940
941 BAGELQDPFloat *resptr = &(d.elem(start).elem().elem(0,0).real());
942 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(start).elem().elem(0,0).real()));
943 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
944
945 qdp_su3_mm_peq(resptr, scalar, lptr, rptr, num_sites, (unsigned long)0);
946
947 }
948 else {
949 // Do site by site
950 const int* tab = s.siteTable().slice();
951 unsigned int num_sites = s.numSiteTable();
952 for(unsigned int j=0; j < num_sites; j++) {
953 int i = tab[j];
954 unsigned long one_site = 1;
955 BAGELQDPFloat *resptr = &(d.elem(i).elem().elem(0,0).real());
956
957 BAGELQDPFloat *lptr = const_cast<BAGELQDPFloat*>(&(l.elem(i).elem().elem(0,0).real()));
958 BAGELQDPFloat *rptr = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
959
960 qdp_su3_mm_peq(resptr, scalar, lptr, rptr, one_site, (unsigned long)0);
961
962 }
963 }
964
965}
966
967
968 // +=a * M
969template<>
970inline
971void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
972 const OpAddAssign& op,
973
974 const QDPExpr<
975 BinaryNode< OpMultiply,
976 Reference<
977 QDPType<
978 PScalar<PScalar<RScalar<BAGELQDPFloat> > >,
979 OScalar<PScalar<PScalar<RScalar<BAGELQDPFloat> > > >
980 >
981 >,
982 Reference<
983 QDPType<
984 PScalar<PColorMatrix<RComplexFloat, 3> >,
985 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
986 >
987 >
988 >,
989 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
990 >& rhs,
991
992 const Subset& s) {
993
994 typedef OScalar<PScalar<PScalar<RScalar<BAGELQDPFloat> > > > F;
995
996 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
997
998
999
1000 const F& scal = static_cast<const F&>(rhs.expression().left());
1001 const C& r = static_cast<const C&>(rhs.expression().right());
1002
1003#if DEBUG_BAGELQDP_LINALG
1004 QDPIO::cout << "evaluate(M += alpha*M ) subset = s " << endl;
1005#endif
1006
1007 BAGELQDPFloat *scalar = const_cast<BAGELQDPFloat*>(&(scal.elem().elem().elem().elem()));
1008 if( s.hasOrderedRep() ) {
1009 // Do whole subset
1010 unsigned int start = s.start();
1011 unsigned int end = s.end();
1012
1013 unsigned long num_sites = end - start + 1;
1014
1015 BAGELQDPFloat *y = &(d.elem(start).elem().elem(0,0).real());
1016
1017 BAGELQDPFloat *x = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
1018
1019 qdp_vaxpy3(y, scalar, x, y, 3*num_sites);
1020
1021 }
1022 else {
1023 // Do site by site
1024 const int* tab = s.siteTable().slice();
1025 unsigned int num_sites = s.numSiteTable();
1026 for(unsigned int j=0; j < num_sites; j++) {
1027 int i = tab[j];
1028 BAGELQDPFloat *y = &(d.elem(i).elem().elem(0,0).real());
1029 BAGELQDPFloat *x = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
1030
1031 qdp_vaxpy3(y, scalar, x,y, 3);
1032
1033 }
1034 }
1035
1036}
1037
1038 // -=a * M
1039template<>
1040inline
1041void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
1042 const OpSubtractAssign& op,
1043
1044 const QDPExpr<
1045 BinaryNode< OpMultiply,
1046 Reference<
1047 QDPType<
1048 PScalar<PScalar<RScalar<BAGELQDPFloat> > >,
1049 OScalar<PScalar<PScalar<RScalar<BAGELQDPFloat> > > >
1050 >
1051 >,
1052 Reference<
1053 QDPType<
1054 PScalar<PColorMatrix<RComplexFloat, 3> >,
1055 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
1056 >
1057 >
1058 >,
1059 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
1060 >& rhs,
1061
1062 const Subset& s) {
1063
1064 typedef OScalar<PScalar<PScalar<RScalar<BAGELQDPFloat> > > > F;
1065
1066 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
1067
1068
1069
1070 const F& scal = static_cast<const F&>(rhs.expression().left());
1071 const C& r = static_cast<const C&>(rhs.expression().right());
1072
1073#if DEBUG_BAGELQDP_LINALG
1074 QDPIO::cout << "evaluate(M -= alpha*M ) subset = s " << endl;
1075#endif
1076
1077 BAGELQDPFloat scalar = -scal.elem().elem().elem().elem();
1078 if( s.hasOrderedRep() ) {
1079 // Do whole subset
1080 unsigned int start = s.start();
1081 unsigned int end = s.end();
1082
1083 unsigned long num_sites = end - start + 1;
1084
1085 BAGELQDPFloat *y = &(d.elem(start).elem().elem(0,0).real());
1086
1087 BAGELQDPFloat *x = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
1088
1089 qdp_vaxpy3(y, &scalar, x, y, 3*num_sites);
1090
1091 }
1092 else {
1093 // Do site by site
1094 const int* tab = s.siteTable().slice();
1095 unsigned int num_sites = s.numSiteTable();
1096 for(unsigned int j=0; j < num_sites; j++) {
1097 int i = tab[j];
1098 BAGELQDPFloat *y = &(d.elem(i).elem().elem(0,0).real());
1099 BAGELQDPFloat *x = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
1100
1101 qdp_vaxpy3(y, &scalar, x,y, 3);
1102
1103 }
1104 }
1105
1106}
1107
1108 // += M
1109template<>
1110inline
1111void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
1112 const OpAddAssign& op,
1113
1114 const QDPExpr<
1115 UnaryNode<OpIdentity,
1116 Reference<
1117 QDPType<
1118 PScalar<PColorMatrix<RComplexFloat, 3> >,
1119 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
1120 >
1121 >
1122 >,
1123 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
1124 >& rhs,
1125
1126 const Subset& s) {
1127
1128 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
1129 const C& r = static_cast<const C&>(rhs.expression().child());
1130
1131#if DEBUG_BAGELQDP_LINALG
1132 QDPIO::cout << "evaluate(M += M ) subset = s " << endl;
1133#endif
1134
1135 if( s.hasOrderedRep() ) {
1136 // Do whole subset
1137 unsigned int start = s.start();
1138 unsigned int end = s.end();
1139
1140 unsigned long num_sites = end - start + 1;
1141
1142 BAGELQDPFloat *y = &(d.elem(start).elem().elem(0,0).real());
1143 BAGELQDPFloat *x = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
1144
1145 qdp_vadd3(y, x, y, 3*num_sites);
1146
1147 }
1148 else {
1149 // Do site by site
1150 const int* tab = s.siteTable().slice();
1151 unsigned int num_sites = s.numSiteTable();
1152 for(unsigned int j=0; j < num_sites; j++) {
1153 int i = tab[j];
1154 BAGELQDPFloat *y = &(d.elem(i).elem().elem(0,0).real());
1155 BAGELQDPFloat *x = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
1156
1157 qdp_vadd3(y,x,y, 3);
1158
1159 }
1160 }
1161
1162}
1163
1164 // -=a * M
1165template<>
1166inline
1167void evaluate(OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >& d,
1168 const OpSubtractAssign& op,
1169
1170 const QDPExpr<
1171 UnaryNode< OpIdentity,
1172 Reference<
1173 QDPType<
1174 PScalar<PColorMatrix<RComplexFloat, 3> >,
1175 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
1176 >
1177 >
1178 >,
1179 OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > >
1180 >& rhs,
1181
1182 const Subset& s) {
1183
1184 typedef OLattice<PScalar<PColorMatrix<RComplexFloat, 3> > > C;
1185
1186 const C& r = static_cast<const C&>(rhs.expression().child());
1187
1188#if DEBUG_BAGELQDP_LINALG
1189 QDPIO::cout << "evaluate(M -= M ) subset = s " << endl;
1190#endif
1191
1192 if( s.hasOrderedRep() ) {
1193 // Do whole subset
1194 unsigned int start = s.start();
1195 unsigned int end = s.end();
1196
1197 unsigned long num_sites = end - start + 1;
1198
1199 BAGELQDPFloat *y = &(d.elem(start).elem().elem(0,0).real());
1200
1201 BAGELQDPFloat *x = const_cast<BAGELQDPFloat*>(&(r.elem(start).elem().elem(0,0).real()));
1202
1203 qdp_vsub3(y, y,x, 3*num_sites);
1204
1205 }
1206 else {
1207 // Do site by site
1208 const int* tab = s.siteTable().slice();
1209 unsigned int num_sites = s.numSiteTable();
1210 for(unsigned int j=0; j < num_sites; j++) {
1211 int i = tab[j];
1212 BAGELQDPFloat *y = &(d.elem(i).elem().elem(0,0).real());
1213 BAGELQDPFloat *x = const_cast<BAGELQDPFloat*>(&(r.elem(i).elem().elem(0,0).real()));
1214
1215 qdp_vsub3(y,y,x, 3);
1216
1217 }
1218 }
1219
1220}
1221
1222
1223#endif
1224
1225
1226#if defined(DEBUG_BAGELQDP_LINALG)
1227#undef DEBUG_BAGELQDP_LINALG
1228#endif
1229
1230} // namespace QDP;
1231
1232#endif
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::PColorMatrix
Primitive color Matrix class.
Definition qdp_primcolormat.h:26
QDP::PScalar
Primitive Scalar.
Definition qdp_primscalar.h:25
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::RComplex
Reality complex.
Definition qdp_reality.h:253
QDP::RScalar
Scalar reality (not complex).
Definition qdp_reality.h:28
QDP::Subset
Subsets - controls how lattices are looped.
Definition qdp_subset.h:39
QDP::Subset::end
int end() const
Definition qdp_subset.h:81
QDP::Subset::siteTable
const multi1d< int > & siteTable() const
Definition qdp_subset.h:83
QDP::Subset::numSiteTable
int numSiteTable() const
Definition qdp_subset.h:84
QDP::Subset::hasOrderedRep
bool hasOrderedRep() const
Definition qdp_subset.h:79
QDP::Subset::start
int start() const
Definition qdp_subset.h:80
QDP::multi1d::slice
const T * slice() const
Return ref to a column slice.
Definition qdp_multi.h:225
QDP::evaluate
void evaluate(OLattice< DCol > &d, const OpAssign &op, const QDPExpr< BinaryNode< OpMultiply, Reference< QDPType< DCol, OLattice< DCol > > >, Reference< QDPType< DCol, OLattice< DCol > > > >, OLattice< DCol > > &rhs, const Subset &s)
Definition qdp_scalarsite_linalg_double.cc:78
QDP::C
OLattice< PScalar< PColorMatrix< RComplexFloat, 3 > > > C
Definition qdp_scalarsite_generic_linalg.h:12
QDP::RComplexFloat
RComplex< BAGELQDPFloat > RComplexFloat
Definition qdp_scalarsite_bagel_qdp_linalg.h:29
QDP::QDPIO::cout
StandardOutputStream cout
Definition qdp_stdio.cc:21
QDP
Yet another random number generator.
Definition qdp_bgq_threadbind.cc:16
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_ALIGN16
#define QDP_ALIGN16
Definition qdp.h:61
Generated by doxygen 1.16.1