QDP++
sse_linalg_m_peq_m_double.cc
Go to the documentation of this file.
1
6
8#include <xmmintrin.h>
9
10namespace QDP {
11
12
13 typedef union {
14 double c[2];
15 __m128d v;
16 } VD;
17
18
19 /* M2 += M1 */
20 void ssed_m_peq_m(REAL64* m2, REAL64* m1, int n_mat)
21 {
22
23 __m128d tmp1;
24 __m128d tmp2;
25 __m128d tmp3;
26 __m128d tmp4;
27 __m128d tmp5;
28 __m128d tmp6;
29 __m128d tmp7;
30 __m128d tmp8;
31 __m128d tmp9;
32 __m128d tmp10;
33 __m128d tmp11;
34 __m128d tmp12;
35 __m128d tmp13;
36 __m128d tmp14;
37 __m128d tmp15;
38
39
40 REAL64* m1_p=m1;
41 REAL64* m2_p=m2;
42
43 for(int i=0; i < n_mat; i++) {
44 tmp1= _mm_loadu_pd(m1_p);
45 tmp2= _mm_loadu_pd(m2_p);
46 tmp3 = _mm_add_pd(tmp1,tmp2);
47 _mm_storeu_pd(m2_p, tmp3);
48
49 tmp4= _mm_loadu_pd(m1_p+2);
50 tmp5= _mm_loadu_pd(m2_p+2);
51 tmp6 = _mm_add_pd(tmp4,tmp5);
52 _mm_storeu_pd(m2_p+2, tmp6);
53
54 tmp7= _mm_loadu_pd(m1_p+4);
55 tmp8= _mm_loadu_pd(m2_p+4);
56 tmp9 = _mm_add_pd(tmp7,tmp8);
57 _mm_storeu_pd(m2_p+4, tmp9);
58
59 tmp10= _mm_loadu_pd(m1_p+6);
60 tmp11= _mm_loadu_pd(m2_p+6);
61 tmp12 = _mm_add_pd(tmp10,tmp11);
62 _mm_storeu_pd(m2_p+6, tmp12);
63
64 tmp13= _mm_loadu_pd(m1_p+8);
65 tmp14= _mm_loadu_pd(m2_p+8);
66 tmp15 = _mm_add_pd(tmp13,tmp14);
67 _mm_storeu_pd(m2_p+8, tmp15);
68
69 tmp1= _mm_loadu_pd(m1_p+10);
70 tmp2= _mm_loadu_pd(m2_p+10);
71 tmp3 = _mm_add_pd(tmp1,tmp2);
72 _mm_storeu_pd(m2_p+10, tmp3);
73
74 tmp4= _mm_loadu_pd(m1_p+12);
75 tmp5= _mm_loadu_pd(m2_p+12);
76 tmp6 = _mm_add_pd(tmp4,tmp5);
77 _mm_storeu_pd(m2_p+12, tmp6);
78
79 tmp7= _mm_loadu_pd(m1_p+14);
80 tmp8= _mm_loadu_pd(m2_p+14);
81 tmp9 = _mm_add_pd(tmp7,tmp8);
82 _mm_storeu_pd(m2_p+14, tmp9);
83
84 tmp10= _mm_loadu_pd(m1_p+16);
85 tmp11= _mm_loadu_pd(m2_p+16);
86 tmp12 = _mm_add_pd(tmp10,tmp11);
87 _mm_storeu_pd(m2_p+16, tmp12);
88
89 m1_p += 18; m2_p+=18;
90 }
91
92 }
93
94 /* M2 -= M1 */
95 void ssed_m_meq_m(REAL64* m2, REAL64* m1, int n_mat)
96 {
97 __m128d tmp1;
98 __m128d tmp2;
99 __m128d tmp3;
100 __m128d tmp4;
101 __m128d tmp5;
102 __m128d tmp6;
103 __m128d tmp7;
104 __m128d tmp8;
105 __m128d tmp9;
106 __m128d tmp10;
107 __m128d tmp11;
108 __m128d tmp12;
109 __m128d tmp13;
110 __m128d tmp14;
111 __m128d tmp15;
112
113
114 REAL64* m1_p=m1;
115 REAL64* m2_p=m2;
116
117 for(int i=0; i < n_mat; i++) {
118 tmp1= _mm_loadu_pd(m1_p);
119 tmp2= _mm_loadu_pd(m2_p);
120 tmp3 = _mm_sub_pd(tmp2,tmp1);
121 _mm_storeu_pd(m2_p, tmp3);
122
123 tmp4= _mm_loadu_pd(m1_p+2);
124 tmp5= _mm_loadu_pd(m2_p+2);
125 tmp6 = _mm_sub_pd(tmp5,tmp4);
126 _mm_storeu_pd(m2_p+2, tmp6);
127
128 tmp7= _mm_loadu_pd(m1_p+4);
129 tmp8= _mm_loadu_pd(m2_p+4);
130 tmp9 = _mm_sub_pd(tmp8,tmp7);
131 _mm_storeu_pd(m2_p+4, tmp9);
132
133 tmp10= _mm_loadu_pd(m1_p+6);
134 tmp11= _mm_loadu_pd(m2_p+6);
135 tmp12 = _mm_sub_pd(tmp11,tmp10);
136 _mm_storeu_pd(m2_p+6, tmp12);
137
138 tmp13= _mm_loadu_pd(m1_p+8);
139 tmp14= _mm_loadu_pd(m2_p+8);
140 tmp15 = _mm_sub_pd(tmp14,tmp13);
141 _mm_storeu_pd(m2_p+8, tmp15);
142
143 tmp1= _mm_loadu_pd(m1_p+10);
144 tmp2= _mm_loadu_pd(m2_p+10);
145 tmp3 = _mm_sub_pd(tmp2,tmp1);
146 _mm_storeu_pd(m2_p+10, tmp3);
147
148 tmp4= _mm_loadu_pd(m1_p+12);
149 tmp5= _mm_loadu_pd(m2_p+12);
150 tmp6 = _mm_sub_pd(tmp5,tmp4);
151 _mm_storeu_pd(m2_p+12, tmp6);
152
153 tmp7= _mm_loadu_pd(m1_p+14);
154 tmp8= _mm_loadu_pd(m2_p+14);
155 tmp9 = _mm_sub_pd(tmp8,tmp7);
156 _mm_storeu_pd(m2_p+14, tmp9);
157
158 tmp10= _mm_loadu_pd(m1_p+16);
159 tmp11= _mm_loadu_pd(m2_p+16);
160 tmp12 = _mm_sub_pd(tmp11,tmp10);
161 _mm_storeu_pd(m2_p+16, tmp12);
162
163 m1_p += 18; m2_p+=18;
164 }
165 }
166
167
168 /* M2 += adj(M1)*/
169 void ssed_m_peq_h(REAL64* m2, REAL64* m1, int n_mat)
170 {
171 __m128d mfact = _mm_set_pd( (REAL64)(-1), (REAL64)(1) );
172
173 __m128d m1_11;
174 __m128d m1_12;
175 __m128d m1_13;
176 __m128d m1_21;
177 __m128d m1_22;
178 __m128d m1_23;
179 __m128d m1_31;
180 __m128d m1_32;
181 __m128d m1_33;
182
183 __m128d tmp1;
184 __m128d tmp2;
185 __m128d tmp3;
186 __m128d tmp4;
187 __m128d tmp5;
188 __m128d tmp6;
189
190
191 REAL64* m1_p=m1;
192 REAL64* m2_p=m2;
193
194 for(int i=0; i < n_mat; i++) {
195 // Stream in m1
196 m1_11= _mm_loadu_pd(m1_p);
197 m1_12= _mm_loadu_pd(m1_p+2);
198 m1_13= _mm_loadu_pd(m1_p+4);
199 m1_21= _mm_loadu_pd(m1_p+6);
200 m1_22= _mm_loadu_pd(m1_p+8);
201 m1_23= _mm_loadu_pd(m1_p+10);
202 m1_31= _mm_loadu_pd(m1_p+12);
203 m1_32= _mm_loadu_pd(m1_p+14);
204 m1_33= _mm_loadu_pd(m1_p+16);
205
206 tmp1 = _mm_loadu_pd(m2_p);
207 tmp2 = _mm_mul_pd(mfact, m1_11);
208 tmp3 = _mm_add_pd(tmp1, tmp2);
209 _mm_storeu_pd(m2_p, tmp3);
210
211 tmp4 = _mm_loadu_pd(m2_p+2);
212 tmp5 = _mm_mul_pd(mfact, m1_21);
213 tmp6 = _mm_add_pd(tmp4, tmp5);
214 _mm_storeu_pd(m2_p+2,tmp6);
215
216 tmp1 = _mm_loadu_pd(m2_p+4);
217 tmp2 = _mm_mul_pd(mfact, m1_31);
218 tmp3 = _mm_add_pd(tmp1, tmp2);
219 _mm_storeu_pd(m2_p+4, tmp3);
220
221 tmp4 = _mm_loadu_pd(m2_p+6);
222 tmp5 = _mm_mul_pd(mfact, m1_12);
223 tmp6 = _mm_add_pd(tmp4, tmp5);
224 _mm_storeu_pd(m2_p+6, tmp6);
225
226
227 tmp1 = _mm_loadu_pd(m2_p+8);
228 tmp2 = _mm_mul_pd(mfact, m1_22);
229 tmp3 = _mm_add_pd(tmp1, tmp2);
230 _mm_storeu_pd(m2_p+8, tmp3);
231
232 tmp4 = _mm_loadu_pd(m2_p+10);
233 tmp5 = _mm_mul_pd(mfact, m1_32);
234 tmp6 = _mm_add_pd(tmp4, tmp5);
235 _mm_storeu_pd(m2_p+10,tmp6);
236
237 tmp4 = _mm_loadu_pd(m2_p+12);
238 tmp5 = _mm_mul_pd(mfact, m1_13);
239 tmp6 = _mm_add_pd(tmp4, tmp5);
240 _mm_storeu_pd(m2_p+12, tmp6);
241
242
243 tmp1 = _mm_loadu_pd(m2_p+14);
244 tmp2 = _mm_mul_pd(mfact, m1_23);
245 tmp3 = _mm_add_pd(tmp1, tmp2);
246 _mm_storeu_pd(m2_p+14, tmp3);
247
248 tmp4 = _mm_loadu_pd(m2_p+16);
249 tmp5 = _mm_mul_pd(mfact, m1_33);
250 tmp6 = _mm_add_pd(tmp4, tmp5);
251 _mm_storeu_pd(m2_p+16, tmp6);
252
253 m1_p += 18; m2_p+=18;
254 }
255
256 }
257
258 /* M2 -= adj(M1) */
259 void ssed_m_meq_h(REAL64* m2, REAL64* m1, int n_mat)
260 {
261 __m128d m1_11;
262 __m128d m1_12;
263 __m128d m1_13;
264 __m128d m1_21;
265 __m128d m1_22;
266 __m128d m1_23;
267 __m128d m1_31;
268 __m128d m1_32;
269 __m128d m1_33;
270
271 __m128d tmp1;
272 __m128d tmp2;
273 __m128d tmp3;
274 __m128d tmp4;
275 __m128d tmp5;
276 __m128d tmp6;
277
278
279 __m128d mfact = _mm_set_pd( (REAL64)(-1), (REAL64)(1));
280
281 REAL64* m1_p=m1;
282 REAL64* m2_p=m2;
283
284 for(int i=0; i < n_mat; i++) {
285 // Stream in m1
286 m1_11= _mm_loadu_pd(m1_p);
287 m1_12= _mm_loadu_pd(m1_p+2);
288 m1_13= _mm_loadu_pd(m1_p+4);
289 m1_21= _mm_loadu_pd(m1_p+6);
290 m1_22= _mm_loadu_pd(m1_p+8);
291 m1_23= _mm_loadu_pd(m1_p+10);
292 m1_31= _mm_loadu_pd(m1_p+12);
293 m1_32= _mm_loadu_pd(m1_p+14);
294 m1_33= _mm_loadu_pd(m1_p+16);
295
296 tmp1 = _mm_loadu_pd(m2_p);
297 tmp2 = _mm_mul_pd(mfact, m1_11);
298 tmp3 = _mm_sub_pd(tmp1, tmp2);
299 _mm_storeu_pd(m2_p, tmp3);
300
301 tmp4 = _mm_loadu_pd(m2_p+2);
302 tmp5 = _mm_mul_pd(mfact, m1_21);
303 tmp6 = _mm_sub_pd(tmp4, tmp5);
304 _mm_storeu_pd(m2_p+2,tmp6);
305
306 tmp1 = _mm_loadu_pd(m2_p+4);
307 tmp2 = _mm_mul_pd(mfact, m1_31);
308 tmp3 = _mm_sub_pd(tmp1, tmp2);
309 _mm_storeu_pd(m2_p+4, tmp3);
310
311 tmp4 = _mm_loadu_pd(m2_p+6);
312 tmp5 = _mm_mul_pd(mfact, m1_12);
313 tmp6 = _mm_sub_pd(tmp4, tmp5);
314 _mm_storeu_pd(m2_p+6, tmp6);
315
316
317 tmp1 = _mm_loadu_pd(m2_p+8);
318 tmp2 = _mm_mul_pd(mfact, m1_22);
319 tmp3 = _mm_sub_pd(tmp1, tmp2);
320 _mm_storeu_pd(m2_p+8, tmp3);
321
322 tmp4 = _mm_loadu_pd(m2_p+10);
323 tmp5 = _mm_mul_pd(mfact, m1_32);
324 tmp6 = _mm_sub_pd(tmp4, tmp5);
325 _mm_storeu_pd(m2_p+10,tmp6);
326
327 tmp4 = _mm_loadu_pd(m2_p+12);
328 tmp5 = _mm_mul_pd(mfact, m1_13);
329 tmp6 = _mm_sub_pd(tmp4, tmp5);
330 _mm_storeu_pd(m2_p+12, tmp6);
331
332
333 tmp1 = _mm_loadu_pd(m2_p+14);
334 tmp2 = _mm_mul_pd(mfact, m1_23);
335 tmp3 = _mm_sub_pd(tmp1, tmp2);
336 _mm_storeu_pd(m2_p+14, tmp3);
337
338 tmp4 = _mm_loadu_pd(m2_p+16);
339 tmp5 = _mm_mul_pd(mfact, m1_33);
340 tmp6 = _mm_sub_pd(tmp4, tmp5);
341 _mm_storeu_pd(m2_p+16, tmp6);
342
343 m1_p += 18; m2_p+=18;
344 }
345 }
346
347
348
349} // namespace QDP;
350
double REAL64
Yet another random number generator.
void ssed_m_meq_m(REAL64 *m2, REAL64 *m1, int n_mat)
void ssed_m_peq_m(REAL64 *m2, REAL64 *m1, int n_mat)
void ssed_m_peq_h(REAL64 *m2, REAL64 *m1, int n_mat)
void ssed_m_meq_h(REAL64 *m2, REAL64 *m1, int n_mat)