QDP++
generic_blas_vcaxpby3.h
Go to the documentation of this file.
1#ifndef QDP_GENERIC_BLAS_VCAXPBY3
2#define QDP_GENERIC_BLAS_VCAXPBY3
3
4namespace QDP {
5inline
6void vcaxpby3(REAL* Out, REAL* ap, REAL* xp, REAL* bp, REAL *yp, int n_3vec)
7{
8 double a_r;
9 double a_i;
10 double b_r;
11 double b_i;
12
13 double x0r;
14 double x0i;
15
16 double x1r;
17 double x1i;
18
19 double x2r;
20 double x2i;
21
22 double y0r;
23 double y0i;
24
25 double y1r;
26 double y1i;
27
28 double y2r;
29 double y2i;
30
31 double z0r;
32 double z0i;
33
34 double z1r;
35 double z1i;
36
37 double z2r;
38 double z2i;
39
40 a_r =(double)(*ap);
41 a_i =(double)*(ap+1);
42
43 b_r =(double)(*bp);
44 b_i =(double)*(bp+1);
45
46 int index_x = 0;
47 int index_y = 0;
48 int index_z = 0;
49
50 int counter;
51
52 if( n_3vec > 0 ) {
53 // Prefetch whole vectors
54 x0r = (double)xp[index_x++];
55 x0i = (double)xp[index_x++];
56 y0r = (double)yp[index_y++];
57 y0i = (double)yp[index_y++];
58
59 x1r = (double)xp[index_x++];
60 x1i = (double)xp[index_x++];
61 y1r = (double)yp[index_y++];
62 y1i = (double)yp[index_y++];
63
64 x2r = (double)xp[index_x++];
65 x2i = (double)xp[index_x++];
66 y2r = (double)yp[index_y++];
67 y2i = (double)yp[index_y++];
68
69 int len = 4*n_3vec;
70 for( counter = 0; counter < len-1; counter++) {
71
72 z0r = a_r * x0r;
73 z0i = a_i * x0r;
74 x0r = (double)xp[index_x++];
75 z0r += b_r * y0r;
76 z0i += b_i * y0r;
77 y0r = (double)yp[index_y++];
78
79 z0r -= a_i * x0i;
80 z0i += a_r * x0i;
81 x0i = (double)xp[index_x++];
82 z0r -= b_i * y0i;
83 z0i += b_r * y0i;
84 y0i = (double)yp[index_y++];
85 Out[index_z++] = (REAL)z0r;
86 Out[index_z++] = (REAL)z0i;
87
88 z1r = a_r * x1r;
89 z1i = a_i * x1r;
90 x1r = (double)xp[index_x++];
91 z1r += b_r * y1r;
92 z1i += b_i * y1r;
93 y1r = (double)yp[index_y++];
94
95 z1r -= a_i * x1i;
96 z1i += a_r * x1i;
97 x1i = (double)xp[index_x++];
98 z1r -= b_i * y1i;
99 z1i += b_r * y1i;
100 y1i = (double)yp[index_y++];
101 Out[index_z++] = (REAL)z1r;
102 Out[index_z++] = (REAL)z1i;
103
104 z2r = a_r * x2r;
105 z2i = a_i * x2r;
106 x2r = (double)xp[index_x++];
107 z2r += b_r * y2r;
108 z2i += b_i * y2r;
109 y2r = (double)yp[index_y++];
110
111 z2r -= a_i * x2i;
112 z2i += a_r * x2i;
113 x2i = (double)xp[index_x++];
114 z2r -= b_i * y2i;
115 z2i += b_r * y2i;
116 y2i = (double)yp[index_y++];
117 Out[index_z++] = (REAL)z2r;
118 Out[index_z++] = (REAL)z2i;
119 }
120
121 z0r = a_r * x0r;
122 z0i = a_i * x0r;
123 z0r += b_r * y0r;
124 z0i += b_i * y0r;
125
126 z0r -= a_i * x0i;
127 z0i += a_r * x0i;
128 z0r -= b_i * y0i;
129 Out[index_z++] = (REAL)z0r;
130 z0i += b_r * y0i;
131 Out[index_z++] = (REAL)z0i;
132
133 z1r = a_r * x1r;
134 z1i = a_i * x1r;
135 z1r += b_r * y1r;
136 z1i += b_i * y1r;
137
138 z1r -= a_i * x1i;
139 z1i += a_r * x1i;
140 z1r -= b_i * y1i;
141 Out[index_z++] = (REAL)z1r;
142 z1i += b_r * y1i;
143 Out[index_z++] = (REAL)z1i;
144
145 z2r = a_r * x2r;
146 z2i = a_i * x2r;
147 z2r += b_r * y2r;
148 z2i += b_i * y2r;
149
150 z2r -= a_i * x2i;
151 z2i += a_r * x2i;
152 z2r -= b_i * y2i;
153 Out[index_z++] = (REAL)z2r;
154 z2i += b_r * y2i;
155 Out[index_z++] = (REAL)z2i;
156 }
157}
158
159
160} // namespace QDP;
161
162#endif // guard
REAL32 REAL
Yet another random number generator.
void vcaxpby3(REAL *Out, REAL *ap, REAL *xp, REAL *bp, REAL *yp, int n_3vec)