QDP++
generic_blas_g5.h
Go to the documentation of this file.
1#ifndef QDP_GENERIC_BLAS_G5
2#define QDP_GENERIC_BLAS_G5
3
4namespace QDP {
5
6// (Vector) out = (Scalar) (*scalep) * (Vector) In
7inline
8void scal_g5(REAL *Out, REAL *scalep, REAL *In, int n_4vec)
9{
10 double a;
11 double x0r;
12 double x0i;
13
14 double x1r;
15 double x1i;
16
17 double x2r;
18 double x2i;
19
20 double z0r;
21 double z0i;
22
23 double z1r;
24 double z1i;
25
26 double z2r;
27 double z2i;
28
29 a = *scalep;
30
31 int index_x = 0;
32 int index_z = 0;
33
34 int counter;
35
36 for( counter = 0; counter < n_4vec; counter++) {
37 // Spin Component 0
38 x0r = (double)In[index_x++];
39 z0r = a*x0r;
40 Out[index_z++] =(REAL) z0r;
41
42 x0i = (double)In[index_x++];
43 z0i = a*x0i;
44 Out[index_z++] =(REAL) z0i;
45
46 x1r = (double)In[index_x++];
47 z1r = a*x1r;
48 Out[index_z++] = (REAL)z1r;
49
50 x1i = (double)In[index_x++];
51 z1i = a*x1i;
52 Out[index_z++] = (REAL)z1i;
53
54 x2r = (double)In[index_x++];
55 z2r = a*x2r;
56 Out[index_z++] = (REAL)z2r;
57
58 x2i = (double)In[index_x++];
59 z2i = a*x2i;
60 Out[index_z++] = (REAL)z2i;
61
62 // Spin Component 1
63 x0r = (double)In[index_x++];
64 z0r = a*x0r;
65 Out[index_z++] =(REAL) z0r;
66
67 x0i = (double)In[index_x++];
68 z0i = a*x0i;
69 Out[index_z++] =(REAL) z0i;
70
71 x1r = (double)In[index_x++];
72 z1r = a*x1r;
73 Out[index_z++] = (REAL)z1r;
74
75 x1i = (double)In[index_x++];
76 z1i = a*x1i;
77 Out[index_z++] = (REAL)z1i;
78
79 x2r = (double)In[index_x++];
80 z2r = a*x2r;
81 Out[index_z++] = (REAL)z2r;
82
83 x2i = (double)In[index_x++];
84 z2i = a*x2i;
85 Out[index_z++] = (REAL)z2i;
86
87 // Spin Component 2
88 x0r = (double)In[index_x++];
89 z0r = a*x0r;
90 Out[index_z++] =-(REAL) z0r;
91
92 x0i = (double)In[index_x++];
93 z0i = a*x0i;
94 Out[index_z++] =-(REAL) z0i;
95
96 x1r = (double)In[index_x++];
97 z1r = a*x1r;
98 Out[index_z++] =-(REAL)z1r;
99
100 x1i = (double)In[index_x++];
101 z1i = a*x1i;
102 Out[index_z++] =-(REAL)z1i;
103
104 x2r = (double)In[index_x++];
105 z2r = a*x2r;
106 Out[index_z++] =-(REAL)z2r;
107
108 x2i = (double)In[index_x++];
109 z2i = a*x2i;
110 Out[index_z++] =-(REAL)z2i;
111
112 // Spin Component 3
113 x0r = (double)In[index_x++];
114 z0r = a*x0r;
115 Out[index_z++] =-(REAL) z0r;
116
117 x0i = (double)In[index_x++];
118 z0i = a*x0i;
119 Out[index_z++] =-(REAL) z0i;
120
121 x1r = (double)In[index_x++];
122 z1r = a*x1r;
123 Out[index_z++] =-(REAL)z1r;
124
125 x1i = (double)In[index_x++];
126 z1i = a*x1i;
127 Out[index_z++] =-(REAL)z1i;
128
129 x2r = (double)In[index_x++];
130 z2r = a*x2r;
131 Out[index_z++] =-(REAL)z2r;
132
133 x2i = (double)In[index_x++];
134 z2i = a*x2i;
135 Out[index_z++] =-(REAL)z2i;
136 }
137}
138
139// (Vector) out = (Scalar) (*scalep) * (Vector) InScale + (scalep2)*g5*vector)Add)
140inline
141void axpbyz_g5(REAL *Out,REAL *scalep,REAL *InScale, REAL *scalep2, REAL *Add,int n_4vec)
142{
143 double a;
144 double b;
145
146 double x0r;
147 double x0i;
148
149 double x1r;
150 double x1i;
151
152 double x2r;
153 double x2i;
154
155 double y0r;
156 double y0i;
157
158 double y1r;
159 double y1i;
160
161 double y2r;
162 double y2i;
163
164 double z0r;
165 double z0i;
166
167 double z1r;
168 double z1i;
169
170 double z2r;
171 double z2i;
172
173 a = *scalep;
174 b = *scalep2;
175
176 int index_x = 0;
177 int index_y = 0;
178 int index_z = 0;
179
180 int counter;
181
182 for( counter = 0; counter < n_4vec; counter++) {
183 // Spin Component 0 (AXPY3)
184 x0r = (double)InScale[index_x++];
185 y0r = (double)Add[index_y++];
186 z0r = a*x0r ;
187 z0r += b*y0r;
188 Out[index_z++] =(REAL) z0r;
189
190 x0i = (double)InScale[index_x++];
191 y0i = (double)Add[index_y++];
192 z0i = a*x0i;
193 z0i += b*y0i;
194 Out[index_z++] =(REAL) z0i;
195
196 x1r = (double)InScale[index_x++];
197 y1r = (double)Add[index_y++];
198 z1r = a*x1r ;
199 z1r += b*y1r;
200 Out[index_z++] = (REAL)z1r;
201
202 x1i = (double)InScale[index_x++];
203 y1i = (double)Add[index_y++];
204 z1i = a*x1i;
205 z1i += b*y1i;
206 Out[index_z++] = (REAL)z1i;
207
208 x2r = (double)InScale[index_x++];
209 y2r = (double)Add[index_y++];
210 z2r = a*x2r ;
211 z2r += b*y2r;
212 Out[index_z++] = (REAL)z2r;
213
214 x2i = (double)InScale[index_x++];
215 y2i = (double)Add[index_y++];
216 z2i = a*x2i ;
217 z2i += b*y2i;
218 Out[index_z++] = (REAL)z2i;
219
220 // Spin Component 1
221 x0r = (double)InScale[index_x++];
222 y0r = (double)Add[index_y++];
223 z0r = a*x0r;
224 z0r += b*y0r;
225 Out[index_z++] =(REAL) z0r;
226
227 x0i = (double)InScale[index_x++];
228 y0i = (double)Add[index_y++];
229 z0i = a*x0i;
230 z0i += b*y0i;
231 Out[index_z++] =(REAL) z0i;
232
233 x1r = (double)InScale[index_x++];
234 y1r = (double)Add[index_y++];
235 z1r = a*x1r ;
236 z1r += b*y1r;
237 Out[index_z++] = (REAL)z1r;
238
239 x1i = (double)InScale[index_x++];
240 y1i = (double)Add[index_y++];
241 z1i = a*x1i;
242 z1i += b*y1i;
243 Out[index_z++] = (REAL)z1i;
244
245 x2r = (double)InScale[index_x++];
246 y2r = (double)Add[index_y++];
247 z2r = a*x2r;
248 z2r += b*y2r;
249 Out[index_z++] = (REAL)z2r;
250
251 x2i = (double)InScale[index_x++];
252 y2i = (double)Add[index_y++];
253 z2i = a*x2i;
254 z2i += b*y2i;
255 Out[index_z++] = (REAL)z2i;
256
257 // Spin Component 2 (AXPY3)
258 x0r = (double)InScale[index_x++];
259 y0r = (double)Add[index_y++];
260 z0r = a*x0r ;
261 z0r -= b*y0r;
262 Out[index_z++] =(REAL) z0r;
263
264 x0i = (double)InScale[index_x++];
265 y0i = (double)Add[index_y++];
266 z0i = a*x0i;
267 z0i -= b*y0i;
268 Out[index_z++] =(REAL) z0i;
269
270 x1r = (double)InScale[index_x++];
271 y1r = (double)Add[index_y++];
272 z1r = a*x1r ;
273 z1r -= b*y1r;
274 Out[index_z++] = (REAL)z1r;
275
276 x1i = (double)InScale[index_x++];
277 y1i = (double)Add[index_y++];
278 z1i = a*x1i;
279 z1i -= b*y1i;
280 Out[index_z++] = (REAL)z1i;
281
282 x2r = (double)InScale[index_x++];
283 y2r = (double)Add[index_y++];
284 z2r = a*x2r ;
285 z2r -= b*y2r;
286 Out[index_z++] = (REAL)z2r;
287
288 x2i = (double)InScale[index_x++];
289 y2i = (double)Add[index_y++];
290 z2i = a*x2i ;
291 z2i -= b*y2i;
292 Out[index_z++] = (REAL)z2i;
293
294 // Spin Component 1
295 x0r = (double)InScale[index_x++];
296 y0r = (double)Add[index_y++];
297 z0r = a*x0r;
298 z0r -= b*y0r;
299 Out[index_z++] =(REAL) z0r;
300
301 x0i = (double)InScale[index_x++];
302 y0i = (double)Add[index_y++];
303 z0i = a*x0i;
304 z0i -= b*y0i;
305 Out[index_z++] =(REAL) z0i;
306
307 x1r = (double)InScale[index_x++];
308 y1r = (double)Add[index_y++];
309 z1r = a*x1r ;
310 z1r -= b*y1r;
311 Out[index_z++] = (REAL)z1r;
312
313 x1i = (double)InScale[index_x++];
314 y1i = (double)Add[index_y++];
315 z1i = a*x1i;
316 z1i -= b*y1i;
317 Out[index_z++] = (REAL)z1i;
318
319 x2r = (double)InScale[index_x++];
320 y2r = (double)Add[index_y++];
321 z2r = a*x2r;
322 z2r -= b*y2r;
323 Out[index_z++] = (REAL)z2r;
324
325 x2i = (double)InScale[index_x++];
326 y2i = (double)Add[index_y++];
327 z2i = a*x2i;
328 z2i -= b*y2i;
329 Out[index_z++] = (REAL)z2i;
330 }
331}
332
333// (Vector) out = (Vector) Add + (Scalar) (*scalep) * (Vector) P{+} InScale
334inline
335void xmayz_g5(REAL *Out,REAL *scalep,REAL *Add, REAL *InScale,int n_4vec)
336{
337 double a;
338 double x0r;
339 double x0i;
340
341 double x1r;
342 double x1i;
343
344 double x2r;
345 double x2i;
346
347 double y0r;
348 double y0i;
349
350 double y1r;
351 double y1i;
352
353 double y2r;
354 double y2i;
355
356 double z0r;
357 double z0i;
358
359 double z1r;
360 double z1i;
361
362 double z2r;
363 double z2i;
364
365 a = *scalep;
366
367 int index_x = 0;
368 int index_y = 0;
369 int index_z = 0;
370
371 int counter;
372
373 for( counter = 0; counter < n_4vec; counter++) {
374 // Spin Component 0 (AYPX)
375 x0r = (double)Add[index_x++];
376 y0r = (double)InScale[index_y++];
377 z0r = x0r - a*y0r;
378 Out[index_z++] =(REAL) z0r;
379
380 x0i = (double)Add[index_x++];
381 y0i = (double)InScale[index_y++];
382 z0i = x0i - a*y0i;
383 Out[index_z++] =(REAL) z0i;
384
385 x1r = (double)Add[index_x++];
386 y1r = (double)InScale[index_y++];
387 z1r = x1r - a*y1r;
388 Out[index_z++] = (REAL)z1r;
389
390 x1i = (double)Add[index_x++];
391 y1i = (double)InScale[index_y++];
392 z1i = x1i - a*y1i;
393 Out[index_z++] = (REAL)z1i;
394
395 x2r = (double)Add[index_x++];
396 y2r = (double)InScale[index_y++];
397 z2r = x2r - a*y2r;
398 Out[index_z++] = (REAL)z2r;
399
400 x2i = (double)Add[index_x++];
401 y2i = (double)InScale[index_y++];
402 z2i = x2i - a*y2i;
403 Out[index_z++] = (REAL)z2i;
404
405 // Spin Component 1 (AYPX)
406 x0r = (double)Add[index_x++];
407 y0r = (double)InScale[index_y++];
408 z0r = x0r - a*y0r;
409 Out[index_z++] =(REAL) z0r;
410
411 x0i = (double)Add[index_x++];
412 y0i = (double)InScale[index_y++];
413 z0i = x0i - a*y0i;
414 Out[index_z++] =(REAL) z0i;
415
416 x1r = (double)Add[index_x++];
417 y1r = (double)InScale[index_y++];
418 z1r = x1r - a*y1r;
419 Out[index_z++] = (REAL)z1r;
420
421 x1i = (double)Add[index_x++];
422 y1i = (double)InScale[index_y++];
423 z1i = x1i - a*y1i;
424 Out[index_z++] = (REAL)z1i;
425
426 x2r = (double)Add[index_x++];
427 y2r = (double)InScale[index_y++];
428 z2r = x2r - a*y2r;
429 Out[index_z++] = (REAL)z2r;
430
431 x2i = (double)Add[index_x++];
432 y2i = (double)InScale[index_y++];
433 z2i = x2i - a*y2i;
434 Out[index_z++] = (REAL)z2i;
435
436 // Spin Component 2 (AYPX)
437 x0r = (double)Add[index_x++];
438 y0r = (double)InScale[index_y++];
439 z0r = x0r + a*y0r;
440 Out[index_z++] =(REAL) z0r;
441
442 x0i = (double)Add[index_x++];
443 y0i = (double)InScale[index_y++];
444 z0i = x0i + a*y0i;
445 Out[index_z++] =(REAL) z0i;
446
447 x1r = (double)Add[index_x++];
448 y1r = (double)InScale[index_y++];
449 z1r = x1r + a*y1r;
450 Out[index_z++] = (REAL)z1r;
451
452 x1i = (double)Add[index_x++];
453 y1i = (double)InScale[index_y++];
454 z1i = x1i + a*y1i;
455 Out[index_z++] = (REAL)z1i;
456
457 x2r = (double)Add[index_x++];
458 y2r = (double)InScale[index_y++];
459 z2r = x2r + a*y2r;
460 Out[index_z++] = (REAL)z2r;
461
462 x2i = (double)Add[index_x++];
463 y2i = (double)InScale[index_y++];
464 z2i = x2i + a*y2i;
465 Out[index_z++] = (REAL)z2i;
466
467 // Spin Component 1 (AYPX)
468 x0r = (double)Add[index_x++];
469 y0r = (double)InScale[index_y++];
470 z0r = x0r + a*y0r;
471 Out[index_z++] =(REAL) z0r;
472
473 x0i = (double)Add[index_x++];
474 y0i = (double)InScale[index_y++];
475 z0i = x0i + a*y0i;
476 Out[index_z++] =(REAL) z0i;
477
478 x1r = (double)Add[index_x++];
479 y1r = (double)InScale[index_y++];
480 z1r = x1r + a*y1r;
481 Out[index_z++] = (REAL)z1r;
482
483 x1i = (double)Add[index_x++];
484 y1i = (double)InScale[index_y++];
485 z1i = x1i + a*y1i;
486 Out[index_z++] = (REAL)z1i;
487
488 x2r = (double)Add[index_x++];
489 y2r = (double)InScale[index_y++];
490 z2r = x2r + a*y2r;
491 Out[index_z++] = (REAL)z2r;
492
493 x2i = (double)Add[index_x++];
494 y2i = (double)InScale[index_y++];
495 z2i = x2i + a*y2i;
496 Out[index_z++] = (REAL)z2i;
497
498
499
500 }
501}
502
503// (Vector) out = (Scalar) (*scalep) * (Vector) InScale + (scalep2)*g5*vector)Add)
504inline
505void g5_axmbyz(REAL *Out,REAL *scalep,REAL *InScale, REAL *scalep2, REAL *Add,int n_4vec)
506{
507 double a;
508 double b;
509
510 double x0r;
511 double x0i;
512
513 double x1r;
514 double x1i;
515
516 double x2r;
517 double x2i;
518
519 double y0r;
520 double y0i;
521
522 double y1r;
523 double y1i;
524
525 double y2r;
526 double y2i;
527
528 double z0r;
529 double z0i;
530
531 double z1r;
532 double z1i;
533
534 double z2r;
535 double z2i;
536
537 a = *scalep;
538 b = *scalep2;
539
540 int index_x = 0;
541 int index_y = 0;
542 int index_z = 0;
543
544 int counter;
545
546 for( counter = 0; counter < n_4vec; counter++) {
547 // Spin Component 0 (AXPY3)
548 x0r = (double)InScale[index_x++];
549 y0r = (double)Add[index_y++];
550 z0r = a*x0r ;
551 z0r -= b*y0r;
552 Out[index_z++] =(REAL) z0r;
553
554 x0i = (double)InScale[index_x++];
555 y0i = (double)Add[index_y++];
556 z0i = a*x0i;
557 z0i -= b*y0i;
558 Out[index_z++] =(REAL) z0i;
559
560 x1r = (double)InScale[index_x++];
561 y1r = (double)Add[index_y++];
562 z1r = a*x1r ;
563 z1r -= b*y1r;
564 Out[index_z++] = (REAL)z1r;
565
566 x1i = (double)InScale[index_x++];
567 y1i = (double)Add[index_y++];
568 z1i = a*x1i;
569 z1i -= b*y1i;
570 Out[index_z++] = (REAL)z1i;
571
572 x2r = (double)InScale[index_x++];
573 y2r = (double)Add[index_y++];
574 z2r = a*x2r ;
575 z2r -= b*y2r;
576 Out[index_z++] = (REAL)z2r;
577
578 x2i = (double)InScale[index_x++];
579 y2i = (double)Add[index_y++];
580 z2i = a*x2i ;
581 z2i -= b*y2i;
582 Out[index_z++] = (REAL)z2i;
583
584 // Spin Component 1
585 x0r = (double)InScale[index_x++];
586 y0r = (double)Add[index_y++];
587 z0r = a*x0r;
588 z0r -= b*y0r;
589 Out[index_z++] =(REAL) z0r;
590
591 x0i = (double)InScale[index_x++];
592 y0i = (double)Add[index_y++];
593 z0i = a*x0i;
594 z0i -= b*y0i;
595 Out[index_z++] =(REAL) z0i;
596
597 x1r = (double)InScale[index_x++];
598 y1r = (double)Add[index_y++];
599 z1r = a*x1r ;
600 z1r -= b*y1r;
601 Out[index_z++] = (REAL)z1r;
602
603 x1i = (double)InScale[index_x++];
604 y1i = (double)Add[index_y++];
605 z1i = a*x1i;
606 z1i -= b*y1i;
607 Out[index_z++] = (REAL)z1i;
608
609 x2r = (double)InScale[index_x++];
610 y2r = (double)Add[index_y++];
611 z2r = a*x2r;
612 z2r -= b*y2r;
613 Out[index_z++] = (REAL)z2r;
614
615 x2i = (double)InScale[index_x++];
616 y2i = (double)Add[index_y++];
617 z2i = a*x2i;
618 z2i -= b*y2i;
619 Out[index_z++] = (REAL)z2i;
620
621 // Spin Component 2 (AXPY3)
622 x0r = (double)InScale[index_x++];
623 y0r = (double)Add[index_y++];
624 z0r = b*y0r;
625 z0r -= a*x0r ;
626
627 Out[index_z++] =(REAL) z0r;
628
629 x0i = (double)InScale[index_x++];
630 y0i = (double)Add[index_y++];
631 z0i = b*y0i;
632 z0i -= a*x0i;
633 Out[index_z++] =(REAL) z0i;
634
635 x1r = (double)InScale[index_x++];
636 y1r = (double)Add[index_y++];
637 z1r = b*y1r;
638 z1r -= a*x1r ;
639 Out[index_z++] = (REAL)z1r;
640
641 x1i = (double)InScale[index_x++];
642 y1i = (double)Add[index_y++];
643 z1i = b*y1i;
644 z1i -= a*x1i;
645 Out[index_z++] = (REAL)z1i;
646
647 x2r = (double)InScale[index_x++];
648 y2r = (double)Add[index_y++];
649 z2r = b*y2r;
650 z2r -= a*x2r ;
651 Out[index_z++] = (REAL)z2r;
652
653 x2i = (double)InScale[index_x++];
654 y2i = (double)Add[index_y++];
655 z2i = b*y2i;
656 z2i -= a*x2i ;
657
658 Out[index_z++] = (REAL)z2i;
659
660 // Spin Component 2 (AXPY3)
661 x0r = (double)InScale[index_x++];
662 y0r = (double)Add[index_y++];
663 z0r = b*y0r;
664 z0r -= a*x0r ;
665
666 Out[index_z++] =(REAL) z0r;
667
668 // Spin Component 3
669 x0i = (double)InScale[index_x++];
670 y0i = (double)Add[index_y++];
671 z0i = b*y0i;
672 z0i -= a*x0i;
673 Out[index_z++] =(REAL) z0i;
674
675 x1r = (double)InScale[index_x++];
676 y1r = (double)Add[index_y++];
677 z1r = b*y1r;
678 z1r -= a*x1r ;
679 Out[index_z++] = (REAL)z1r;
680
681 x1i = (double)InScale[index_x++];
682 y1i = (double)Add[index_y++];
683 z1i = b*y1i;
684 z1i -= a*x1i;
685 Out[index_z++] = (REAL)z1i;
686
687 x2r = (double)InScale[index_x++];
688 y2r = (double)Add[index_y++];
689 z2r = b*y2r;
690 z2r -= a*x2r ;
691 Out[index_z++] = (REAL)z2r;
692
693 x2i = (double)InScale[index_x++];
694 y2i = (double)Add[index_y++];
695 z2i = b*y2i;
696 z2i -= a*x2i ;
697
698 Out[index_z++] = (REAL)z2i;
699 }
700}
701
702// (Vector) out = (Scalar) (*scalep) * (Vector) InScale + (scalep2)*ig5*vector)Add)
703inline
704void axpbyz_ig5(REAL *Out,REAL *scalep,REAL *InScale, REAL *scalep2, REAL *Add,int n_4vec)
705{
706 double a;
707 double b;
708
709 double x0r;
710 double x0i;
711
712 double x1r;
713 double x1i;
714
715 double x2r;
716 double x2i;
717
718 double y0r;
719 double y0i;
720
721 double y1r;
722 double y1i;
723
724 double y2r;
725 double y2i;
726
727 double z0r;
728 double z0i;
729
730 double z1r;
731 double z1i;
732
733 double z2r;
734 double z2i;
735
736 a = *scalep;
737 b = *scalep2;
738
739 int index_x = 0;
740 int index_y = 0;
741 int index_z = 0;
742
743 int counter;
744
745 for( counter = 0; counter < n_4vec; counter++) {
746
747 // Spin Component 0 (AXPY3)
748 x0r = (double)InScale[index_x++];
749 x0i = (double)InScale[index_x++];
750 y0r = (double)Add[index_y++];
751 y0i = (double)Add[index_y++];
752
753 z0r = a*x0r ;
754 z0r -= b*y0i;
755 Out[index_z++] =(REAL) z0r;
756 z0i = a*x0i;
757 z0i += b*y0r;
758 Out[index_z++] =(REAL) z0i;
759
760 x1r = (double)InScale[index_x++];
761 x1i = (double)InScale[index_x++];
762 y1r = (double)Add[index_y++];
763 y1i = (double)Add[index_y++];
764
765 z1r = a*x1r ;
766 z1r -= b*y1i;
767 Out[index_z++] = (REAL)z1r;
768 z1i = a*x1i;
769 z1i += b*y1r;
770 Out[index_z++] = (REAL)z1i;
771
772 x2r = (double)InScale[index_x++];
773 x2i = (double)InScale[index_x++];
774 y2r = (double)Add[index_y++];
775 y2i = (double)Add[index_y++];
776
777 z2r = a*x2r ;
778 z2r -= b*y2i;
779 Out[index_z++] = (REAL)z2r;
780 z2i = a*x2i ;
781 z2i += b*y2r;
782 Out[index_z++] = (REAL)z2i;
783
784 // Spin Component 1
785 x0r = (double)InScale[index_x++];
786 x0i = (double)InScale[index_x++];
787 y0r = (double)Add[index_y++];
788 y0i = (double)Add[index_y++];
789
790 z0r = a*x0r ;
791 z0r -= b*y0i;
792 Out[index_z++] =(REAL) z0r;
793 z0i = a*x0i;
794 z0i += b*y0r;
795 Out[index_z++] =(REAL) z0i;
796
797 x1r = (double)InScale[index_x++];
798 x1i = (double)InScale[index_x++];
799 y1r = (double)Add[index_y++];
800 y1i = (double)Add[index_y++];
801
802 z1r = a*x1r ;
803 z1r -= b*y1i;
804 Out[index_z++] = (REAL)z1r;
805 z1i = a*x1i;
806 z1i += b*y1r;
807 Out[index_z++] = (REAL)z1i;
808
809 x2r = (double)InScale[index_x++];
810 x2i = (double)InScale[index_x++];
811 y2r = (double)Add[index_y++];
812 y2i = (double)Add[index_y++];
813
814 z2r = a*x2r ;
815 z2r -= b*y2i;
816 Out[index_z++] = (REAL)z2r;
817 z2i = a*x2i ;
818 z2i += b*y2r;
819 Out[index_z++] = (REAL)z2i;
820
821 // Spin Component 2
822 x0r = (double)InScale[index_x++];
823 x0i = (double)InScale[index_x++];
824 y0r = (double)Add[index_y++];
825 y0i = (double)Add[index_y++];
826
827 z0r = a*x0r ;
828 z0r += b*y0i;
829 Out[index_z++] =(REAL) z0r;
830 z0i = a*x0i;
831 z0i -= b*y0r;
832 Out[index_z++] =(REAL) z0i;
833
834 x1r = (double)InScale[index_x++];
835 x1i = (double)InScale[index_x++];
836 y1r = (double)Add[index_y++];
837 y1i = (double)Add[index_y++];
838
839 z1r = a*x1r ;
840 z1r += b*y1i;
841 Out[index_z++] = (REAL)z1r;
842 z1i = a*x1i;
843 z1i -= b*y1r;
844 Out[index_z++] = (REAL)z1i;
845
846 x2r = (double)InScale[index_x++];
847 x2i = (double)InScale[index_x++];
848 y2r = (double)Add[index_y++];
849 y2i = (double)Add[index_y++];
850
851 z2r = a*x2r ;
852 z2r += b*y2i;
853 Out[index_z++] = (REAL)z2r;
854 z2i = a*x2i ;
855 z2i -= b*y2r;
856 Out[index_z++] = (REAL)z2i;
857
858 // Spin Component 3
859 x0r = (double)InScale[index_x++];
860 x0i = (double)InScale[index_x++];
861 y0r = (double)Add[index_y++];
862 y0i = (double)Add[index_y++];
863
864 z0r = a*x0r ;
865 z0r += b*y0i;
866 Out[index_z++] =(REAL) z0r;
867 z0i = a*x0i;
868 z0i -= b*y0r;
869 Out[index_z++] =(REAL) z0i;
870
871 x1r = (double)InScale[index_x++];
872 x1i = (double)InScale[index_x++];
873 y1r = (double)Add[index_y++];
874 y1i = (double)Add[index_y++];
875
876 z1r = a*x1r ;
877 z1r += b*y1i;
878 Out[index_z++] = (REAL)z1r;
879 z1i = a*x1i;
880 z1i -= b*y1r;
881 Out[index_z++] = (REAL)z1i;
882
883 x2r = (double)InScale[index_x++];
884 x2i = (double)InScale[index_x++];
885 y2r = (double)Add[index_y++];
886 y2i = (double)Add[index_y++];
887
888
889 z2r = a*x2r ;
890 z2r += b*y2i;
891 Out[index_z++] = (REAL)z2r;
892 z2i = a*x2i ;
893 z2i -= b*y2r;
894 Out[index_z++] = (REAL)z2i;
895
896
897 }
898}
899
900
901// (Vector) out = (Scalar) (*scalep) * (Vector) InScale - (scalep2)*ig5*vector)Add)
902inline
903void axmbyz_ig5(REAL *Out,REAL *scalep,REAL *InScale, REAL *scalep2, REAL *Add,int n_4vec)
904{
905 double a;
906 double b;
907
908 double x0r;
909 double x0i;
910
911 double x1r;
912 double x1i;
913
914 double x2r;
915 double x2i;
916
917 double y0r;
918 double y0i;
919
920 double y1r;
921 double y1i;
922
923 double y2r;
924 double y2i;
925
926 double z0r;
927 double z0i;
928
929 double z1r;
930 double z1i;
931
932 double z2r;
933 double z2i;
934
935 a = *scalep;
936 b = *scalep2;
937
938 int index_x = 0;
939 int index_y = 0;
940 int index_z = 0;
941
942 int counter;
943
944 for( counter = 0; counter < n_4vec; counter++) {
945
946 // Spin Component 0 (AXPY3)
947 x0r = (double)InScale[index_x++];
948 x0i = (double)InScale[index_x++];
949 y0r = (double)Add[index_y++];
950 y0i = (double)Add[index_y++];
951
952 z0r = a*x0r ;
953 z0r += b*y0i;
954 Out[index_z++] =(REAL) z0r;
955 z0i = a*x0i;
956 z0i -= b*y0r;
957 Out[index_z++] =(REAL) z0i;
958
959 x1r = (double)InScale[index_x++];
960 x1i = (double)InScale[index_x++];
961 y1r = (double)Add[index_y++];
962 y1i = (double)Add[index_y++];
963
964 z1r = a*x1r ;
965 z1r += b*y1i;
966 Out[index_z++] = (REAL)z1r;
967 z1i = a*x1i;
968 z1i -= b*y1r;
969 Out[index_z++] = (REAL)z1i;
970
971 x2r = (double)InScale[index_x++];
972 x2i = (double)InScale[index_x++];
973 y2r = (double)Add[index_y++];
974 y2i = (double)Add[index_y++];
975
976 z2r = a*x2r ;
977 z2r += b*y2i;
978 Out[index_z++] = (REAL)z2r;
979 z2i = a*x2i ;
980 z2i -= b*y2r;
981 Out[index_z++] = (REAL)z2i;
982
983 // Spin Component 1
984 x0r = (double)InScale[index_x++];
985 x0i = (double)InScale[index_x++];
986 y0r = (double)Add[index_y++];
987 y0i = (double)Add[index_y++];
988
989 z0r = a*x0r ;
990 z0r += b*y0i;
991 Out[index_z++] =(REAL) z0r;
992 z0i = a*x0i;
993 z0i -= b*y0r;
994 Out[index_z++] =(REAL) z0i;
995
996 x1r = (double)InScale[index_x++];
997 x1i = (double)InScale[index_x++];
998 y1r = (double)Add[index_y++];
999 y1i = (double)Add[index_y++];
1000
1001 z1r = a*x1r ;
1002 z1r += b*y1i;
1003 Out[index_z++] = (REAL)z1r;
1004 z1i = a*x1i;
1005 z1i -= b*y1r;
1006 Out[index_z++] = (REAL)z1i;
1007
1008 x2r = (double)InScale[index_x++];
1009 x2i = (double)InScale[index_x++];
1010 y2r = (double)Add[index_y++];
1011 y2i = (double)Add[index_y++];
1012
1013 z2r = a*x2r ;
1014 z2r += b*y2i;
1015 Out[index_z++] = (REAL)z2r;
1016 z2i = a*x2i ;
1017 z2i -= b*y2r;
1018 Out[index_z++] = (REAL)z2i;
1019
1020 // Spin Component 2
1021 x0r = (double)InScale[index_x++];
1022 x0i = (double)InScale[index_x++];
1023 y0r = (double)Add[index_y++];
1024 y0i = (double)Add[index_y++];
1025
1026 z0r = a*x0r ;
1027 z0r -= b*y0i;
1028 Out[index_z++] =(REAL) z0r;
1029 z0i = a*x0i;
1030 z0i += b*y0r;
1031 Out[index_z++] =(REAL) z0i;
1032
1033 x1r = (double)InScale[index_x++];
1034 x1i = (double)InScale[index_x++];
1035 y1r = (double)Add[index_y++];
1036 y1i = (double)Add[index_y++];
1037
1038 z1r = a*x1r ;
1039 z1r -= b*y1i;
1040 Out[index_z++] = (REAL)z1r;
1041 z1i = a*x1i;
1042 z1i += b*y1r;
1043 Out[index_z++] = (REAL)z1i;
1044
1045 x2r = (double)InScale[index_x++];
1046 x2i = (double)InScale[index_x++];
1047 y2r = (double)Add[index_y++];
1048 y2i = (double)Add[index_y++];
1049
1050 z2r = a*x2r ;
1051 z2r -= b*y2i;
1052 Out[index_z++] = (REAL)z2r;
1053 z2i = a*x2i ;
1054 z2i += b*y2r;
1055 Out[index_z++] = (REAL)z2i;
1056
1057 // Spin Component 3
1058 x0r = (double)InScale[index_x++];
1059 x0i = (double)InScale[index_x++];
1060 y0r = (double)Add[index_y++];
1061 y0i = (double)Add[index_y++];
1062
1063 z0r = a*x0r ;
1064 z0r -= b*y0i;
1065 Out[index_z++] =(REAL) z0r;
1066 z0i = a*x0i;
1067 z0i += b*y0r;
1068 Out[index_z++] =(REAL) z0i;
1069
1070 x1r = (double)InScale[index_x++];
1071 x1i = (double)InScale[index_x++];
1072 y1r = (double)Add[index_y++];
1073 y1i = (double)Add[index_y++];
1074
1075 z1r = a*x1r ;
1076 z1r -= b*y1i;
1077 Out[index_z++] = (REAL)z1r;
1078 z1i = a*x1i;
1079 z1i += b*y1r;
1080 Out[index_z++] = (REAL)z1i;
1081
1082 x2r = (double)InScale[index_x++];
1083 x2i = (double)InScale[index_x++];
1084 y2r = (double)Add[index_y++];
1085 y2i = (double)Add[index_y++];
1086
1087
1088 z2r = a*x2r ;
1089 z2r -= b*y2i;
1090 Out[index_z++] = (REAL)z2r;
1091 z2i = a*x2i ;
1092 z2i += b*y2r;
1093 Out[index_z++] = (REAL)z2i;
1094 }
1095}
1096
1097
1098
1099// (Vector) out = (Vector) InScale + (scalep)*ig5*vector)Add)
1100inline
1101void xpayz_ig5(REAL *Out,REAL *scalep, REAL *InScale, REAL *Add,int n_4vec)
1102{
1103 double a;
1104 // double b;
1105
1106 double x0r;
1107 double x0i;
1108
1109 double x1r;
1110 double x1i;
1111
1112 double x2r;
1113 double x2i;
1114
1115 double y0r;
1116 double y0i;
1117
1118 double y1r;
1119 double y1i;
1120
1121 double y2r;
1122 double y2i;
1123
1124 double z0r;
1125 double z0i;
1126
1127 double z1r;
1128 double z1i;
1129
1130 double z2r;
1131 double z2i;
1132
1133 a = *scalep;
1134
1135 int index_x = 0;
1136 int index_y = 0;
1137 int index_z = 0;
1138
1139 int counter;
1140
1141 for( counter = 0; counter < n_4vec; counter++) {
1142
1143 // Spin Component 0 (AXPY3)
1144 x0r = (double)InScale[index_x++];
1145 x0i = (double)InScale[index_x++];
1146 y0r = (double)Add[index_y++];
1147 y0i = (double)Add[index_y++];
1148
1149 z0r = x0r - a*y0i ;
1150 Out[index_z++] =(REAL) z0r;
1151 z0i = x0i + a*y0r ;
1152 Out[index_z++] =(REAL) z0i;
1153
1154 x1r = (double)InScale[index_x++];
1155 x1i = (double)InScale[index_x++];
1156 y1r = (double)Add[index_y++];
1157 y1i = (double)Add[index_y++];
1158
1159 z1r = x1r - a*y1i ;
1160 Out[index_z++] = (REAL)z1r;
1161 z1i = x1i + a*y1r;
1162 Out[index_z++] = (REAL)z1i;
1163
1164 x2r = (double)InScale[index_x++];
1165 x2i = (double)InScale[index_x++];
1166 y2r = (double)Add[index_y++];
1167 y2i = (double)Add[index_y++];
1168
1169 z2r = x2r - a*y2i;
1170 Out[index_z++] = (REAL)z2r;
1171 z2i = x2i + a*y2r;
1172 Out[index_z++] = (REAL)z2i;
1173
1174 // Spin Component 1
1175 x0r = (double)InScale[index_x++];
1176 x0i = (double)InScale[index_x++];
1177 y0r = (double)Add[index_y++];
1178 y0i = (double)Add[index_y++];
1179
1180 z0r = x0r - a*y0i ;
1181 Out[index_z++] =(REAL) z0r;
1182 z0i = x0i + a*y0r ;
1183 Out[index_z++] =(REAL) z0i;
1184
1185 x1r = (double)InScale[index_x++];
1186 x1i = (double)InScale[index_x++];
1187 y1r = (double)Add[index_y++];
1188 y1i = (double)Add[index_y++];
1189
1190 z1r = x1r - a*y1i ;
1191 Out[index_z++] = (REAL)z1r;
1192 z1i = x1i + a*y1r;
1193 Out[index_z++] = (REAL)z1i;
1194
1195 x2r = (double)InScale[index_x++];
1196 x2i = (double)InScale[index_x++];
1197 y2r = (double)Add[index_y++];
1198 y2i = (double)Add[index_y++];
1199
1200 z2r = x2r - a*y2i;
1201 Out[index_z++] = (REAL)z2r;
1202 z2i = x2i + a*y2r;
1203 Out[index_z++] = (REAL)z2i;
1204
1205 // Spin Component 2
1206 x0r = (double)InScale[index_x++];
1207 x0i = (double)InScale[index_x++];
1208 y0r = (double)Add[index_y++];
1209 y0i = (double)Add[index_y++];
1210
1211 z0r = x0r + a*y0i ;
1212 Out[index_z++] =(REAL) z0r;
1213 z0i = x0i - a*y0r ;
1214 Out[index_z++] =(REAL) z0i;
1215
1216 x1r = (double)InScale[index_x++];
1217 x1i = (double)InScale[index_x++];
1218 y1r = (double)Add[index_y++];
1219 y1i = (double)Add[index_y++];
1220
1221 z1r = x1r + a*y1i ;
1222 Out[index_z++] = (REAL)z1r;
1223 z1i = x1i - a*y1r;
1224 Out[index_z++] = (REAL)z1i;
1225
1226 x2r = (double)InScale[index_x++];
1227 x2i = (double)InScale[index_x++];
1228 y2r = (double)Add[index_y++];
1229 y2i = (double)Add[index_y++];
1230
1231 z2r = x2r + a*y2i;
1232 Out[index_z++] = (REAL)z2r;
1233 z2i = x2i - a*y2r;
1234 Out[index_z++] = (REAL)z2i;
1235
1236 // Spin Component 3
1237 x0r = (double)InScale[index_x++];
1238 x0i = (double)InScale[index_x++];
1239 y0r = (double)Add[index_y++];
1240 y0i = (double)Add[index_y++];
1241
1242 z0r = x0r + a*y0i ;
1243 Out[index_z++] =(REAL) z0r;
1244 z0i = x0i - a*y0r ;
1245 Out[index_z++] =(REAL) z0i;
1246
1247 x1r = (double)InScale[index_x++];
1248 x1i = (double)InScale[index_x++];
1249 y1r = (double)Add[index_y++];
1250 y1i = (double)Add[index_y++];
1251
1252 z1r = x1r + a*y1i ;
1253 Out[index_z++] = (REAL)z1r;
1254 z1i = x1i - a*y1r;
1255 Out[index_z++] = (REAL)z1i;
1256
1257 x2r = (double)InScale[index_x++];
1258 x2i = (double)InScale[index_x++];
1259 y2r = (double)Add[index_y++];
1260 y2i = (double)Add[index_y++];
1261
1262 z2r = x2r + a*y2i;
1263 Out[index_z++] = (REAL)z2r;
1264 z2i = x2i - a*y2r;
1265 Out[index_z++] = (REAL)z2i;
1266 }
1267}
1268
1269// (Vector) out = (Vector) InScale - (scalep)*ig5*vector)Add)
1270inline
1271void xmayz_ig5(REAL *Out,REAL *scalep, REAL *InScale, REAL *Add,int n_4vec)
1272{
1273 double a;
1274 // double b;
1275
1276 double x0r;
1277 double x0i;
1278
1279 double x1r;
1280 double x1i;
1281
1282 double x2r;
1283 double x2i;
1284
1285 double y0r;
1286 double y0i;
1287
1288 double y1r;
1289 double y1i;
1290
1291 double y2r;
1292 double y2i;
1293
1294 double z0r;
1295 double z0i;
1296
1297 double z1r;
1298 double z1i;
1299
1300 double z2r;
1301 double z2i;
1302
1303 a = *scalep;
1304
1305 int index_x = 0;
1306 int index_y = 0;
1307 int index_z = 0;
1308
1309 int counter;
1310
1311 for( counter = 0; counter < n_4vec; counter++) {
1312
1313 // Spin Component 0 (AXPY3)
1314 x0r = (double)InScale[index_x++];
1315 x0i = (double)InScale[index_x++];
1316 y0r = (double)Add[index_y++];
1317 y0i = (double)Add[index_y++];
1318
1319 z0r = x0r + a*y0i ;
1320 Out[index_z++] =(REAL) z0r;
1321 z0i = x0i - a*y0r ;
1322 Out[index_z++] =(REAL) z0i;
1323
1324 x1r = (double)InScale[index_x++];
1325 x1i = (double)InScale[index_x++];
1326 y1r = (double)Add[index_y++];
1327 y1i = (double)Add[index_y++];
1328
1329 z1r = x1r + a*y1i ;
1330 Out[index_z++] = (REAL)z1r;
1331 z1i = x1i - a*y1r;
1332 Out[index_z++] = (REAL)z1i;
1333
1334 x2r = (double)InScale[index_x++];
1335 x2i = (double)InScale[index_x++];
1336 y2r = (double)Add[index_y++];
1337 y2i = (double)Add[index_y++];
1338
1339 z2r = x2r + a*y2i;
1340 Out[index_z++] = (REAL)z2r;
1341 z2i = x2i - a*y2r;
1342 Out[index_z++] = (REAL)z2i;
1343
1344 // Spin Component 1
1345 x0r = (double)InScale[index_x++];
1346 x0i = (double)InScale[index_x++];
1347 y0r = (double)Add[index_y++];
1348 y0i = (double)Add[index_y++];
1349
1350 z0r = x0r + a*y0i ;
1351 Out[index_z++] =(REAL) z0r;
1352 z0i = x0i - a*y0r ;
1353 Out[index_z++] =(REAL) z0i;
1354
1355 x1r = (double)InScale[index_x++];
1356 x1i = (double)InScale[index_x++];
1357 y1r = (double)Add[index_y++];
1358 y1i = (double)Add[index_y++];
1359
1360 z1r = x1r + a*y1i ;
1361 Out[index_z++] = (REAL)z1r;
1362 z1i = x1i - a*y1r;
1363 Out[index_z++] = (REAL)z1i;
1364
1365 x2r = (double)InScale[index_x++];
1366 x2i = (double)InScale[index_x++];
1367 y2r = (double)Add[index_y++];
1368 y2i = (double)Add[index_y++];
1369
1370 z2r = x2r + a*y2i;
1371 Out[index_z++] = (REAL)z2r;
1372 z2i = x2i - a*y2r;
1373 Out[index_z++] = (REAL)z2i;
1374
1375 // Spin Component 2
1376 x0r = (double)InScale[index_x++];
1377 x0i = (double)InScale[index_x++];
1378 y0r = (double)Add[index_y++];
1379 y0i = (double)Add[index_y++];
1380
1381 z0r = x0r - a*y0i ;
1382 Out[index_z++] =(REAL) z0r;
1383 z0i = x0i + a*y0r ;
1384 Out[index_z++] =(REAL) z0i;
1385
1386 x1r = (double)InScale[index_x++];
1387 x1i = (double)InScale[index_x++];
1388 y1r = (double)Add[index_y++];
1389 y1i = (double)Add[index_y++];
1390
1391 z1r = x1r - a*y1i ;
1392 Out[index_z++] = (REAL)z1r;
1393 z1i = x1i + a*y1r;
1394 Out[index_z++] = (REAL)z1i;
1395
1396 x2r = (double)InScale[index_x++];
1397 x2i = (double)InScale[index_x++];
1398 y2r = (double)Add[index_y++];
1399 y2i = (double)Add[index_y++];
1400
1401 z2r = x2r - a*y2i;
1402 Out[index_z++] = (REAL)z2r;
1403 z2i = x2i + a*y2r;
1404 Out[index_z++] = (REAL)z2i;
1405
1406 // Spin Component 3
1407 x0r = (double)InScale[index_x++];
1408 x0i = (double)InScale[index_x++];
1409 y0r = (double)Add[index_y++];
1410 y0i = (double)Add[index_y++];
1411
1412 z0r = x0r - a*y0i ;
1413 Out[index_z++] =(REAL) z0r;
1414 z0i = x0i + a*y0r ;
1415 Out[index_z++] =(REAL) z0i;
1416
1417 x1r = (double)InScale[index_x++];
1418 x1i = (double)InScale[index_x++];
1419 y1r = (double)Add[index_y++];
1420 y1i = (double)Add[index_y++];
1421
1422 z1r = x1r - a*y1i ;
1423 Out[index_z++] = (REAL)z1r;
1424 z1i = x1i + a*y1r;
1425 Out[index_z++] = (REAL)z1i;
1426
1427 x2r = (double)InScale[index_x++];
1428 x2i = (double)InScale[index_x++];
1429 y2r = (double)Add[index_y++];
1430 y2i = (double)Add[index_y++];
1431
1432 z2r = x2r - a*y2i;
1433 Out[index_z++] = (REAL)z2r;
1434 z2i = x2i + a*y2r;
1435 Out[index_z++] = (REAL)z2i;
1436 }
1437}
1438
1439} // namespace QDP;
1440
1441#endif // guard
REAL32 REAL
Yet another random number generator.
void axpbyz_ig5(REAL *Out, REAL *scalep, REAL *InScale, REAL *scalep2, REAL *Add, int n_4vec)
void xpayz_ig5(REAL *Out, REAL *scalep, REAL *InScale, REAL *Add, int n_4vec)
void xmayz_g5(REAL *Out, REAL *scalep, REAL *Add, REAL *InScale, int n_4vec)
void axmbyz_ig5(REAL *Out, REAL *scalep, REAL *InScale, REAL *scalep2, REAL *Add, int n_4vec)
void xmayz_ig5(REAL *Out, REAL *scalep, REAL *InScale, REAL *Add, int n_4vec)
void axpbyz_g5(REAL *Out, REAL *scalep, REAL *InScale, REAL *scalep2, REAL *Add, int n_4vec)
void g5_axmbyz(REAL *Out, REAL *scalep, REAL *InScale, REAL *scalep2, REAL *Add, int n_4vec)
void scal_g5(REAL *Out, REAL *scalep, REAL *In, int n_4vec)