QDP++
sse_blas_vaxpbyz4_double.cc
Go to the documentation of this file.
1
6
7#include "scalarsite_sse/sse_blas_vaxpbyz4_double.h"
8#include <xmmintrin.h>
9#include "scalarsite_sse/sse_prefetch.h"
10
11namespace QDP {
12
13
14
15#ifndef L2BY2
16#define L2BY2 1365 /* L2 / 2 in SPINORS */
17#endif
18
19
20void vaxpbyz4(REAL64 *z, REAL64 *a, REAL64 *x, REAL64 *b, REAL64 *y, int n_4vec)
21{
22 __m128d a_sse;
23 __m128d b_sse;
24
25 __m128d tmp1;
26 __m128d tmp2;
27 __m128d tmp3;
28 __m128d tmp4;
29
30 __m128d x1;
31 __m128d y1;
32 __m128d z1;
33
34 __m128d x2;
35 __m128d y2;
36 __m128d z2;
37
38 __m128d x3;
39 __m128d y3;
40 __m128d z3;
41
42 // Load the scalar into low bytes of scalar
43 a_sse = _mm_load_sd(a);
44 b_sse = _mm_load_sd(b);
45
46 // cross components into tmp
47 // Zero tmp
48 tmp1 = _mm_setzero_pd();
49
50 tmp1 = _mm_shuffle_pd(a_sse, a_sse, 0x1);
51 a_sse = _mm_add_pd(a_sse, tmp1);
52
53 tmp2 = _mm_setzero_pd();
54 tmp2 = _mm_shuffle_pd(b_sse, b_sse, 0x1);
55 b_sse = _mm_add_pd(b_sse, tmp2);
56
57
58 // Do n_3vec 3vectors.
59 double *x_p=x;
60 double *y_p=y;
61 double *z_p=z;
62
63 if( n_4vec < L2BY2) {
64
65 for(int i=0; i < 3*n_4vec; i++) {
66 PREFETCHW(((const char *)z_p)+16);
67
68 y1 = _mm_load_pd(y_p);
69 x1 = _mm_load_pd(x_p);
70 z1 = _mm_mul_pd(a_sse,x1);
71 tmp1 = _mm_mul_pd(b_sse,y1);
72 z1 = _mm_add_pd(z1,tmp1);
73 _mm_store_pd(z_p, z1);
74
75 y2 = _mm_load_pd(y_p+2);
76 x2 = _mm_load_pd(x_p+2);
77 z2 = _mm_mul_pd(a_sse,x2);
78 tmp2 = _mm_mul_pd(b_sse,y2);
79 z2 = _mm_add_pd(z2,tmp2);
80 _mm_store_pd(z_p+2, z2);
81
82 y3 = _mm_load_pd(y_p+4);
83 x3 = _mm_load_pd(x_p+4);
84 z3 = _mm_mul_pd(a_sse,x3);
85 tmp3 = _mm_mul_pd(b_sse,y3);
86 z3 = _mm_add_pd(z3,tmp3);
87 _mm_store_pd(z_p+4, z3);
88
89
90 y1 = _mm_load_pd(y_p+6);
91 x1 = _mm_load_pd(x_p+6);
92 z1 = _mm_mul_pd(a_sse,x1);
93 tmp1 = _mm_mul_pd(b_sse,y1);
94 z1 = _mm_add_pd(z1,tmp1);
95 _mm_store_pd(z_p+6, z1);
96
97 x_p+=8; y_p+=8; z_p+=8;
98 }
99
100 }
101 else {
102 for(int i=0; i < 3*n_4vec; i++) {
103 PREFETCHNTA(((const char *)x_p)+56);
104 PREFETCHNTA(((const char *)y_p)+56);
105
106 y1 = _mm_load_pd(y_p);
107 x1 = _mm_load_pd(x_p);
108 tmp1 = _mm_mul_pd(b_sse,y1);
109 z1 = _mm_mul_pd(a_sse,x1);
110 z1 = _mm_add_pd(z1,tmp1);
111 _mm_stream_pd(z_p, z1);
112
113 y2 = _mm_load_pd(y_p+2);
114 x2 = _mm_load_pd(x_p+2);
115 tmp2 = _mm_mul_pd(b_sse,y2);
116 z2 = _mm_mul_pd(a_sse,x2);
117 z2 = _mm_add_pd(z2,tmp2);
118 _mm_stream_pd(z_p+2, z2);
119
120 y3 = _mm_load_pd(y_p+4);
121 x3 = _mm_load_pd(x_p+4);
122 tmp3 = _mm_mul_pd(b_sse,y3);
123 z3 = _mm_mul_pd(a_sse,x3);
124 z3 = _mm_add_pd(z3,tmp3);
125 _mm_stream_pd(z_p+4, z3);
126
127
128 y1 = _mm_load_pd(y_p+6);
129 x1 = _mm_load_pd(x_p+6);
130 tmp1 = _mm_mul_pd(b_sse,y1);
131 z1 = _mm_mul_pd(a_sse,x1);
132 z1 = _mm_add_pd(z1,tmp1);
133 _mm_stream_pd(z_p+6, z1);
134
135 x_p +=8; y_p+=8; z_p+=8;
136 }
137 }
138}
139
140
141
142void vaxpby4(REAL64 *y, REAL64 *a, REAL64 *x, REAL64 *b, int n_4vec)
143{
144 __m128d a_sse;
145 __m128d b_sse;
146
147 __m128d tmp1;
148 __m128d tmp2;
149 __m128d tmp3;
150 __m128d tmp4;
151
152 __m128d x1;
153 __m128d y1;
154 __m128d z1;
155
156 __m128d x2;
157 __m128d y2;
158 __m128d z2;
159
160 __m128d x3;
161 __m128d y3;
162 __m128d z3;
163
164 // Load the scalar into low bytes of scalar
165 a_sse = _mm_load_sd(a);
166 b_sse = _mm_load_sd(b);
167
168 // cross components into tmp
169 // Zero tmp
170 tmp1 = _mm_setzero_pd();
171 tmp1 = _mm_shuffle_pd(a_sse, a_sse, 0x1);
172 a_sse = _mm_add_pd(a_sse, tmp1);
173
174 tmp2 = _mm_setzero_pd();
175 tmp2 = _mm_shuffle_pd(b_sse, b_sse, 0x1);
176 b_sse = _mm_add_pd(b_sse, tmp2);
177
178
179 // Do n_3vec 3vectors.
180 double *x_p=x;
181 double *y_p=y;
182
183 if (n_4vec < L2BY2 ) {
184
185 for(int i=0; i < 3*n_4vec; i++) {
186 PREFETCHW(((const char *)y_p)+16);
187 y1 = _mm_load_pd(y_p);
188 x1 = _mm_load_pd(x_p);
189 tmp1 = _mm_mul_pd(b_sse,y1);
190 z1 = _mm_mul_pd(a_sse,x1);
191 z1 = _mm_add_pd(z1,tmp1);
192 _mm_store_pd(y_p, z1);
193
194 y2 = _mm_load_pd(y_p+2);
195 x2 = _mm_load_pd(x_p+2);
196 tmp2 = _mm_mul_pd(b_sse,y2);
197 z2 = _mm_mul_pd(a_sse,x2);
198 z2 = _mm_add_pd(z2,tmp2);
199 _mm_store_pd(y_p+2, z2);
200
201 y3 = _mm_load_pd(y_p+4);
202 x3 = _mm_load_pd(x_p+4);
203 tmp3 = _mm_mul_pd(b_sse,y3);
204 z3 = _mm_mul_pd(a_sse,x3);
205 z3 = _mm_add_pd(z3,tmp3);
206 _mm_store_pd(y_p+4, z3);
207
208 y1 = _mm_load_pd(y_p+6);
209 x1 = _mm_load_pd(x_p+6);
210 z1 = _mm_mul_pd(a_sse,x1);
211 tmp1 = _mm_mul_pd(b_sse,y1);
212 z1 = _mm_add_pd(z1,tmp1);
213 _mm_store_pd(y_p+6, z1);
214
215 x_p += 8;
216 y_p += 8;
217
218 }
219 }
220 else {
221 for(int i=0; i < 3*n_4vec; i++) {
222
223 PREFETCHNTA(((const char *)x_p)+56); // Assume
224 PREFETCHNTA(((const char *)y_p)+56);
225
226 y1 = _mm_load_pd(y_p);
227 x1 = _mm_load_pd(x_p);
228 z1 = _mm_mul_pd(a_sse,x1);
229 tmp1 = _mm_mul_pd(b_sse,y1);
230 z1 = _mm_add_pd(z1,tmp1);
231 _mm_store_pd(y_p, z1);
232
233 y2 = _mm_load_pd(y_p+2);
234 x2 = _mm_load_pd(x_p+2);
235 z2 = _mm_mul_pd(a_sse,x2);
236 tmp2 = _mm_mul_pd(b_sse,y2);
237 z2 = _mm_add_pd(z2,tmp2);
238 _mm_store_pd(y_p+2, z2);
239
240 y3 = _mm_load_pd(y_p+4);
241 x3 = _mm_load_pd(x_p+4);
242 z3 = _mm_mul_pd(a_sse,x3);
243 tmp3 = _mm_mul_pd(b_sse,y3);
244 z3 = _mm_add_pd(z3,tmp3);
245 _mm_store_pd(y_p+4, z3);
246
247 y1 = _mm_load_pd(y_p+6);
248 x1 = _mm_load_pd(x_p+6);
249 z1 = _mm_mul_pd(a_sse,x1);
250 tmp1 = _mm_mul_pd(b_sse,y1);
251 z1 = _mm_add_pd(z1,tmp1);
252 _mm_store_pd(y_p+6, z1);
253
254 x_p+=8;
255 y_p+=8;
256
257#if 0
258 PREFETCHNTA(((const char*)x_p)+80);
259 PREFETCHNTA(((const char*)y_p)+80);
260
261 y2 = _mm_load_pd(y_p+8);
262 x2 = _mm_load_pd(x_p+8);
263 z2 = _mm_mul_pd(a_sse,x2);
264 tmp2 = _mm_mul_pd(b_sse,y2);
265 z2 = _mm_add_pd(z2,tmp2);
266 _mm_store_pd(y_p+8, z2);
267
268 y3 = _mm_load_pd(y_p+10);
269 x3 = _mm_load_pd(x_p+10);
270 z3 = _mm_mul_pd(a_sse,x3);
271 tmp3 = _mm_mul_pd(b_sse,y3);
272 z3 = _mm_add_pd(z3,tmp3);
273 _mm_store_pd(y_p+10, z3);
274
275 y1 = _mm_load_pd(y_p+12);
276 x1 = _mm_load_pd(x_p+12);
277 z1 = _mm_mul_pd(a_sse,x1);
278 tmp1 = _mm_mul_pd(b_sse,y1);
279 z1 = _mm_add_pd(z1,tmp1);
280 _mm_store_pd(y_p+12, z1);
281
282 y2 = _mm_load_pd(y_p+14);
283 x2 = _mm_load_pd(x_p+14);
284 z2 = _mm_mul_pd(a_sse,x2);
285 tmp2 = _mm_mul_pd(b_sse,y2);
286 z2 = _mm_add_pd(z2,tmp2);
287 _mm_store_pd(y_p+14, z2);
288
289 PREFETCHNTA(((const char *)x_p)+88);
290 PREFETCHNTA(((const char *)y_p)+88);
291
292 y3 = _mm_load_pd(y_p+16);
293 x3 = _mm_load_pd(x_p+16);
294 z3 = _mm_mul_pd(a_sse,x3);
295 tmp3 = _mm_mul_pd(b_sse,y3);
296 z3 = _mm_add_pd(z3,tmp3);
297 _mm_store_pd(y_p+16, z3);
298
299
300 y1 = _mm_load_pd(y_p+18);
301 x1 = _mm_load_pd(x_p+18);
302 z1 = _mm_mul_pd(a_sse,x1);
303 tmp1 = _mm_mul_pd(b_sse,y1);
304 z1 = _mm_add_pd(z1,tmp1);
305 _mm_store_pd(y_p+18, z1);
306
307 y2 = _mm_load_pd(y_p+20);
308 x2 = _mm_load_pd(x_p+20);
309 z2 = _mm_mul_pd(a_sse,x2);
310 tmp2 = _mm_mul_pd(b_sse,y2);
311 z2 = _mm_add_pd(z2,tmp2);
312 _mm_store_pd(y_p+20, z2);
313
314
315 y3 = _mm_load_pd(y_p+22);
316 x3 = _mm_load_pd(x_p+22);
317 z3 = _mm_mul_pd(a_sse,x3);
318 tmp3 = _mm_mul_pd(b_sse,y3);
319 z3 = _mm_add_pd(z3,tmp3);
320 _mm_store_pd(y_p+22, z3);
321
322
323 x_p+=24; y_p+=24;
324#endif
325
326 }
327
328 }
329}
330
331
332} // namespace QDP;
QDP::REAL64
double REAL64
Definition qdp_scalarsite_defs.h:22
QDP
Yet another random number generator.
Definition qdp_bgq_threadbind.cc:16
QDP::vaxpby4
void vaxpby4(REAL64 *y, REAL64 *a, REAL64 *x, REAL64 *b, int n_4vec)
Definition sse_blas_vaxpbyz4_double.cc:142
QDP::vaxpbyz4
void vaxpbyz4(REAL64 *z, REAL64 *a, REAL64 *x, REAL64 *b, REAL64 *y, int n_4vec)
Definition sse_blas_vaxpbyz4_double.cc:20
L2BY2
#define L2BY2
Definition sse_blas_vaxmbyz4_double.cc:15
sse_blas_vaxpbyz4_double.h
Generic Scalar VAXPY routine.
sse_prefetch.h
PREFETCHW
#define PREFETCHW(a)
Definition sse_prefetch.h:21
PREFETCHNTA
#define PREFETCHNTA(a)
Definition sse_prefetch.h:9
Generated by doxygen 1.16.1