Line data Source code
1 : // This file is part of ELPA.
2 : //
3 : // The ELPA library was originally created by the ELPA consortium,
4 : // consisting of the following organizations:
5 : //
6 : // - Max Planck Computing and Data Facility (MPCDF), formerly known as
7 : // Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
8 : // - Bergische Universität Wuppertal, Lehrstuhl für angewandte
9 : // Informatik,
10 : // - Technische Universität München, Lehrstuhl für Informatik mit
11 : // Schwerpunkt Wissenschaftliches Rechnen ,
12 : // - Fritz-Haber-Institut, Berlin, Abt. Theorie,
13 : // - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
14 : // Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
15 : // and
16 : // - IBM Deutschland GmbH
17 : //
18 : // This particular source code file contains additions, changes and
19 : // enhancements authored by Intel Corporation which is not part of
20 : // the ELPA consortium.
21 : //
22 : // More information can be found here:
23 : // http://elpa.mpcdf.mpg.de/
24 : //
25 : // ELPA is free software: you can redistribute it and/or modify
26 : // it under the terms of the version 3 of the license of the
27 : // GNU Lesser General Public License as published by the Free
28 : // Software Foundation.
29 : //
30 : // ELPA is distributed in the hope that it will be useful,
31 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
32 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33 : // GNU Lesser General Public License for more details.
34 : //
35 : // You should have received a copy of the GNU Lesser General Public License
36 : // along with ELPA. If not, see <http://www.gnu.org/licenses/>
37 : //
38 : // ELPA reflects a substantial effort on the part of the original
39 : // ELPA consortium, and we ask you to respect the spirit of the
40 : // license that we chose: i.e., please contribute any changes you
41 : // may have back to the original ELPA library distribution, and keep
42 : // any derivatives of ELPA under the same license that we chose for
43 : // the original distribution, the GNU Lesser General Public License.
44 : //
45 : // Author: Andreas Marek, MPCDF, based on the double precision case of A. Heinecke
46 : //
47 : #include "config-f90.h"
48 :
49 : #ifdef HAVE_SSE_INTRINSICS
50 : #include <x86intrin.h>
51 : #endif
52 : #ifdef HAVE_SPARC64_SSE
53 : #include <fjmfunc.h>
54 : #include <emmintrin.h>
55 : #endif
56 : #include <stdio.h>
57 : #include <stdlib.h>
58 :
59 :
60 : #define __forceinline __attribute__((always_inline)) static
61 : #ifdef DOUBLE_PRECISION_REAL
62 : #define __SSE_DATATYPE __m128d
63 : #define _SSE_LOAD _mm_load_pd
64 : #define _SSE_ADD _mm_add_pd
65 : #define _SSE_MUL _mm_mul_pd
66 : #define _SSE_XOR _mm_xor_pd
67 : #define _SSE_STORE _mm_store_pd
68 : #define offset 2
69 : #endif
70 :
71 : #ifdef SINGLE_PRECISION_REAL
72 : #define __SSE_DATATYPE __m128
73 : #define _SSE_LOAD _mm_load_ps
74 : #define _SSE_ADD _mm_add_ps
75 : #define _SSE_MUL _mm_mul_ps
76 : #define _SSE_XOR _mm_xor_ps
77 : #define _SSE_STORE _mm_store_ps
78 : #define offset 4
79 : #endif
80 :
81 :
82 : #ifdef HAVE_SSE_INTRINSICS
83 : #undef __AVX__
84 : #endif
85 :
86 : #ifdef HAVE_SSE_INTRINSICS
87 : //Forward declaration
88 : #ifdef DOUBLE_PRECISION_REAL
89 : __forceinline void hh_trafo_kernel_2_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
90 : __forceinline void hh_trafo_kernel_4_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
91 : __forceinline void hh_trafo_kernel_6_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
92 : __forceinline void hh_trafo_kernel_8_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
93 : __forceinline void hh_trafo_kernel_10_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
94 : __forceinline void hh_trafo_kernel_12_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
95 : #endif
96 : #ifdef SINGLE_PRECISION_REAL
97 : __forceinline void hh_trafo_kernel_4_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
98 : __forceinline void hh_trafo_kernel_8_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
99 : __forceinline void hh_trafo_kernel_12_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
100 : __forceinline void hh_trafo_kernel_16_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
101 : __forceinline void hh_trafo_kernel_20_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
102 : __forceinline void hh_trafo_kernel_24_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
103 : #endif
104 : #endif
105 :
106 : #ifdef HAVE_SPARC64_SSE
107 : //Forward declaration
108 : #ifdef DOUBLE_PRECISION_REAL
109 : __forceinline void hh_trafo_kernel_2_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
110 : __forceinline void hh_trafo_kernel_4_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
111 : __forceinline void hh_trafo_kernel_6_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
112 : __forceinline void hh_trafo_kernel_8_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
113 : __forceinline void hh_trafo_kernel_10_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
114 : __forceinline void hh_trafo_kernel_12_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s);
115 : #endif
116 : #ifdef SINGLE_PRECISION_REAL
117 : __forceinline void hh_trafo_kernel_4_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
118 : __forceinline void hh_trafo_kernel_8_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
119 : __forceinline void hh_trafo_kernel_12_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
120 : __forceinline void hh_trafo_kernel_16_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
121 : __forceinline void hh_trafo_kernel_20_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
122 : __forceinline void hh_trafo_kernel_24_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s);
123 : #endif
124 : #endif
125 :
126 :
127 : #ifdef HAVE_SSE_INTRINSICS
128 : #ifdef DOUBLE_PRECISION_REAL
129 : void double_hh_trafo_real_sse_2hv_double(double* q, double* hh, int* pnb, int* pnq, int* pldq, int* pldh);
130 : #endif
131 : #ifdef SINGLE_PRECISION_REAL
132 : void double_hh_trafo_real_sse_2hv_single_(float* q, float* hh, int* pnb, int* pnq, int* pldq, int* pldh);
133 : #endif
134 : #endif
135 :
136 : #ifdef HAVE_SPARC64_SSE
137 : #ifdef DOUBLE_PRECISION_REAL
138 : void double_hh_trafo_real_sparc64_2hv_double(double* q, double* hh, int* pnb, int* pnq, int* pldq, int* pldh);
139 : #endif
140 : #ifdef SINGLE_PRECISION_REAL
141 : void double_hh_trafo_real_sparc64_2hv_single_(float* q, float* hh, int* pnb, int* pnq, int* pldq, int* pldh);
142 : #endif
143 : #endif
144 :
145 : /*
146 : !f>#ifdef HAVE_SPARC64_SSE
147 : !f> interface
148 : !f> subroutine double_hh_trafo_real_sparc64_2hv_double(q, hh, pnb, pnq, pldq, pldh) &
149 : !f> bind(C, name="double_hh_trafo_real_sparc64_2hv_double")
150 : !f> use, intrinsic :: iso_c_binding
151 : !f> integer(kind=c_int) :: pnb, pnq, pldq, pldh
152 : !f> type(c_ptr), value :: q
153 : !f> real(kind=c_double) :: hh(pnb,6)
154 : !f> end subroutine
155 : !f> end interface
156 : !f>#endif
157 : */
158 :
159 : /*
160 : !f>#ifdef HAVE_SPARC64_SSE
161 : !f> interface
162 : !f> subroutine double_hh_trafo_real_sparc64_2hv_single(q, hh, pnb, pnq, pldq, pldh) &
163 : !f> bind(C, name="double_hh_trafo_real_sparc64_2hv_single")
164 : !f> use, intrinsic :: iso_c_binding
165 : !f> integer(kind=c_int) :: pnb, pnq, pldq, pldh
166 : !f> type(c_ptr), value :: q
167 : !f> real(kind=c_float) :: hh(pnb,6)
168 : !f> end subroutine
169 : !f> end interface
170 : !f>#endif
171 : */
172 : /*
173 : !f>#ifdef HAVE_SSE_INTRINSICS
174 : !f> interface
175 : !f> subroutine double_hh_trafo_real_sse_2hv_double(q, hh, pnb, pnq, pldq, pldh) &
176 : !f> bind(C, name="double_hh_trafo_real_sse_2hv_double")
177 : !f> use, intrinsic :: iso_c_binding
178 : !f> integer(kind=c_int) :: pnb, pnq, pldq, pldh
179 : !f> type(c_ptr), value :: q
180 : !f> real(kind=c_double) :: hh(pnb,6)
181 : !f> end subroutine
182 : !f> end interface
183 : !f>#endif
184 : */
185 :
186 : /*
187 : !f>#ifdef HAVE_SSE_INTRINSICS
188 : !f> interface
189 : !f> subroutine double_hh_trafo_real_sse_2hv_single(q, hh, pnb, pnq, pldq, pldh) &
190 : !f> bind(C, name="double_hh_trafo_real_sse_2hv_single")
191 : !f> use, intrinsic :: iso_c_binding
192 : !f> integer(kind=c_int) :: pnb, pnq, pldq, pldh
193 : !f> type(c_ptr), value :: q
194 : !f> real(kind=c_float) :: hh(pnb,6)
195 : !f> end subroutine
196 : !f> end interface
197 : !f>#endif
198 : */
199 :
200 : #ifdef HAVE_SSE_INTRINSICS
201 : #ifdef DOUBLE_PRECISION_REAL
202 544768 : void double_hh_trafo_real_sse_2hv_double(double* q, double* hh, int* pnb, int* pnq, int* pldq, int* pldh)
203 : #endif
204 : #ifdef SINGLE_PRECISION_REAL
205 102144 : void double_hh_trafo_real_sse_2hv_single(float* q, float* hh, int* pnb, int* pnq, int* pldq, int* pldh)
206 : #endif
207 : #endif
208 : #ifdef HAVE_SPARC64_SSE
209 : #ifdef DOUBLE_PRECISION_REAL
210 : void double_hh_trafo_real_sparc64_2hv_double(double* q, double* hh, int* pnb, int* pnq, int* pldq, int* pldh)
211 : #endif
212 : #ifdef SINGLE_PRECISION_REAL
213 : void double_hh_trafo_real_sparc64_2hv_single(float* q, float* hh, int* pnb, int* pnq, int* pldq, int* pldh)
214 : #endif
215 : #endif
216 : {
217 : int i;
218 646912 : int nb = *pnb;
219 646912 : int nq = *pldq;
220 646912 : int ldq = *pldq;
221 646912 : int ldh = *pldh;
222 : int worked_on;
223 :
224 646912 : worked_on = 0;
225 :
226 : // calculating scalar product to compute
227 : // 2 householder vectors simultaneously
228 : #ifdef DOUBLE_PRECISION_REAL
229 544768 : double s = hh[(ldh)+1]*1.0;
230 : #endif
231 : #ifdef SINGLE_PRECISION_REAL
232 102144 : float s = hh[(ldh)+1]*1.0;
233 : #endif
234 :
235 : #ifdef HAVE_SSE_INTRINSICS
236 : #pragma ivdep
237 : #endif
238 40755456 : for (i = 2; i < nb; i++)
239 : {
240 40108544 : s += hh[i-1] * hh[(i+ldh)];
241 : }
242 :
243 : // Production level kernel calls with padding
244 : #ifdef DOUBLE_PRECISION_REAL
245 2179072 : for (i = 0; i < nq-10; i+=12)
246 : {
247 : #ifdef HAVE_SSE_INTRINSICS
248 1634304 : hh_trafo_kernel_12_SSE_2hv_double(&q[i], hh, nb, ldq, ldh, s);
249 : #endif
250 : #ifdef HAVE_SPARC64_SSE
251 : hh_trafo_kernel_12_SPARC64_2hv_double(&q[i], hh, nb, ldq, ldh, s);
252 : #endif
253 1634304 : worked_on += 12;
254 : }
255 : #endif
256 : #ifdef SINGLE_PRECISION_REAL
257 408576 : for (i = 0; i < nq-20; i+=24)
258 : {
259 : #ifdef HAVE_SSE_INTRINSICS
260 306432 : hh_trafo_kernel_24_SSE_2hv_single(&q[i], hh, nb, ldq, ldh, s);
261 : #endif
262 : #ifdef HAVE_SPARC64_SSE
263 : hh_trafo_kernel_24_SPARC64_2hv_single(&q[i], hh, nb, ldq, ldh, s);
264 : #endif
265 306432 : worked_on += 24;
266 : }
267 : #endif
268 :
269 646912 : if (nq == i)
270 : {
271 0 : return;
272 : }
273 :
274 : #ifdef DOUBLE_PRECISION_REAL
275 544768 : if (nq-i == 10)
276 : {
277 : #ifdef HAVE_SSE_INTRINSICS
278 0 : hh_trafo_kernel_10_SSE_2hv_double(&q[i], hh, nb, ldq, ldh, s);
279 : #endif
280 : #ifdef HAVE_SPARC64_SSE
281 : hh_trafo_kernel_10_SPARC64_2hv_double(&q[i], hh, nb, ldq, ldh, s);
282 : #endif
283 0 : worked_on += 10;
284 : }
285 : #endif
286 :
287 : #ifdef SINGLE_PRECISION_REAL
288 102144 : if (nq-i == 20)
289 : {
290 : #ifdef HAVE_SSE_INTRINSICS
291 0 : hh_trafo_kernel_20_SSE_2hv_single(&q[i], hh, nb, ldq, ldh, s);
292 : #endif
293 : #ifdef HAVE_SPARC64_SSE
294 : hh_trafo_kernel_20_SPARC64_2hv_single(&q[i], hh, nb, ldq, ldh, s);
295 : #endif
296 0 : worked_on += 20;
297 : }
298 : #endif
299 :
300 : #ifdef DOUBLE_PRECISION_REAL
301 544768 : if (nq-i == 8)
302 : {
303 : #ifdef HAVE_SSE_INTRINSICS
304 0 : hh_trafo_kernel_8_SSE_2hv_double(&q[i], hh, nb, ldq, ldh, s);
305 : #endif
306 : #ifdef HAVE_SPARC64_SSE
307 : hh_trafo_kernel_8_SPARC64_2hv_double(&q[i], hh, nb, ldq, ldh, s);
308 : #endif
309 0 : worked_on += 8;
310 : }
311 : #endif
312 :
313 : #ifdef SINGLE_PRECISION_REAL
314 102144 : if (nq-i == 16)
315 : {
316 : #ifdef HAVE_SSE_INTRINSICS
317 0 : hh_trafo_kernel_16_SSE_2hv_single(&q[i], hh, nb, ldq, ldh, s);
318 : #endif
319 : #ifdef HAVE_SPARC64_SSE
320 : hh_trafo_kernel_16_SPARC64_2hv_single(&q[i], hh, nb, ldq, ldh, s);
321 : #endif
322 0 : worked_on += 16;
323 : }
324 : #endif
325 :
326 :
327 : #ifdef DOUBLE_PRECISION_REAL
328 544768 : if (nq-i == 6)
329 : {
330 : #ifdef HAVE_SSE_INTRINSICS
331 0 : hh_trafo_kernel_6_SSE_2hv_double(&q[i], hh, nb, ldq, ldh, s);
332 : #endif
333 : #ifdef HAVE_SPARC64_SSE
334 : hh_trafo_kernel_6_SPARC64_2hv_double(&q[i], hh, nb, ldq, ldh, s);
335 : #endif
336 0 : worked_on += 6;
337 : }
338 : #endif
339 :
340 : #ifdef SINGLE_PRECISION_REAL
341 102144 : if (nq-i == 12)
342 : {
343 : #ifdef HAVE_SSE_INTRINSICS
344 0 : hh_trafo_kernel_12_SSE_2hv_single(&q[i], hh, nb, ldq, ldh, s);
345 : #endif
346 : #ifdef HAVE_SPARC64_SSE
347 : hh_trafo_kernel_12_SPARC64_2hv_single(&q[i], hh, nb, ldq, ldh, s);
348 : #endif
349 :
350 0 : worked_on += 12;
351 : }
352 : #endif
353 :
354 : #ifdef DOUBLE_PRECISION_REAL
355 544768 : if (nq-i == 4)
356 : {
357 : #ifdef HAVE_SSE_INTRINSICS
358 544768 : hh_trafo_kernel_4_SSE_2hv_double(&q[i], hh, nb, ldq, ldh, s);
359 : #endif
360 : #ifdef HAVE_SPARC64_SSE
361 : hh_trafo_kernel_4_SPARC64_2hv_double(&q[i], hh, nb, ldq, ldh, s);
362 : #endif
363 544768 : worked_on += 4;
364 : }
365 : #endif
366 :
367 : #ifdef SINGLE_PRECISION_REAL
368 102144 : if (nq-i == 8)
369 : {
370 : #ifdef HAVE_SSE_INTRINSICS
371 102144 : hh_trafo_kernel_8_SSE_2hv_single(&q[i], hh, nb, ldq, ldh, s);
372 : #endif
373 : #ifdef HAVE_SPARC64_SSE
374 : hh_trafo_kernel_8_SPARC64_2hv_single(&q[i], hh, nb, ldq, ldh, s);
375 : #endif
376 102144 : worked_on += 8;
377 : }
378 : #endif
379 :
380 : #ifdef DOUBLE_PRECISION_REAL
381 544768 : if (nq-i == 2)
382 : {
383 : #ifdef HAVE_SSE_INTRINSICS
384 0 : hh_trafo_kernel_2_SSE_2hv_double(&q[i], hh, nb, ldq, ldh, s);
385 : #endif
386 : #ifdef HAVE_SPARC64_SSE
387 : hh_trafo_kernel_2_SPARC64_2hv_double(&q[i], hh, nb, ldq, ldh, s);
388 : #endif
389 0 : worked_on += 2;
390 : }
391 : #endif
392 :
393 : #ifdef SINGLE_PRECISION_REAL
394 102144 : if (nq-i == 4)
395 : {
396 : #ifdef HAVE_SSE_INTRINSICS
397 0 : hh_trafo_kernel_4_SSE_2hv_single(&q[i], hh, nb, ldq, ldh, s);
398 : #endif
399 : #ifdef HAVE_SPARC64_SSE
400 : hh_trafo_kernel_4_SPARC64_2hv_single(&q[i], hh, nb, ldq, ldh, s);
401 : #endif
402 :
403 0 : worked_on += 4;
404 : }
405 : #endif
406 : #ifdef WITH_DEBUG
407 : if (worked_on != nq)
408 : {
409 : #ifdef HAVE_SSE_INTRINSICS
410 : printf("Error in real SSE BLOCK2 kernel %d %d\n", worked_on, nq);
411 : #endif
412 : #ifdef HAVE_SPARC64_SSE
413 : printf("Error in real SPARC64 BLOCK2 kernel %d %d\n", worked_on, nq);
414 : #endif
415 :
416 : abort();
417 : }
418 : #endif
419 : }
420 :
421 : /**
422 : * Unrolled kernel that computes
423 : #ifdef DOUBLE_PRECISION_REAL
424 : * 12 rows of Q simultaneously, a
425 : #endif
426 : #ifdef SINGLE_PRECISION_REAL
427 : * 24 rows of Q simultaneously, a
428 : #endif
429 : * matrix Vector product with two householder
430 : * vectors + a rank 2 update is performed
431 : */
432 : #ifdef HAVE_SSE_INTRINSICS
433 : #ifdef DOUBLE_PRECISION_REAL
434 : __forceinline void hh_trafo_kernel_12_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
435 : #endif
436 : #ifdef SINGLE_PRECISION_REAL
437 : __forceinline void hh_trafo_kernel_24_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
438 : #endif
439 : #endif
440 : #ifdef HAVE_SPARC64_SSE
441 : #ifdef DOUBLE_PRECISION_REAL
442 : __forceinline void hh_trafo_kernel_12_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
443 : #endif
444 : #ifdef SINGLE_PRECISION_REAL
445 : __forceinline void hh_trafo_kernel_24_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
446 : #endif
447 : #endif
448 : {
449 : /////////////////////////////////////////////////////
450 : // Matrix Vector Multiplication, Q [12 x nb+1] * hh
451 : // hh contains two householder vectors, with offset 1
452 : /////////////////////////////////////////////////////
453 : int i;
454 : #ifdef HAVE_SSE_INTRINSICS
455 : // Needed bit mask for floating point sign flip
456 : #ifdef DOUBLE_PRECISION_REAL
457 1634304 : __SSE_DATATYPE sign = (__SSE_DATATYPE)_mm_set1_epi64x(0x8000000000000000LL);
458 : #endif
459 : #ifdef SINGLE_PRECISION_REAL
460 612864 : __SSE_DATATYPE sign = _mm_castsi128_ps(_mm_set_epi32(0x80000000, 0x80000000, 0x80000000, 0x80000000));
461 : #endif
462 : #endif
463 :
464 3881472 : __SSE_DATATYPE x1 = _SSE_LOAD(&q[ldq]);
465 3881472 : __SSE_DATATYPE x2 = _SSE_LOAD(&q[ldq+offset]);
466 3881472 : __SSE_DATATYPE x3 = _SSE_LOAD(&q[ldq+2*offset]);
467 3881472 : __SSE_DATATYPE x4 = _SSE_LOAD(&q[ldq+3*offset]);
468 3881472 : __SSE_DATATYPE x5 = _SSE_LOAD(&q[ldq+4*offset]);
469 3881472 : __SSE_DATATYPE x6 = _SSE_LOAD(&q[ldq+5*offset]);
470 :
471 : #ifdef HAVE_SSE_INTRINSICS
472 : #ifdef DOUBLE_PRECISION_REAL
473 3268608 : __SSE_DATATYPE h1 = _mm_set1_pd(hh[ldh+1]);
474 : #endif
475 : #ifdef SINGLE_PRECISION_REAL
476 612864 : __SSE_DATATYPE h1 = _mm_set1_ps(hh[ldh+1]);
477 : #endif
478 : #endif
479 : #ifdef HAVE_SPARC64_SSE
480 : #ifdef DOUBLE_PRECISION_REAL
481 : __SSE_DATATYPE h1 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
482 : #endif
483 : #ifdef SINGLE_PRECISION_REAL
484 : __SSE_DATATYPE h1 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
485 : #endif
486 : #endif
487 :
488 : __SSE_DATATYPE h2;
489 :
490 1940736 : __SSE_DATATYPE q1 = _SSE_LOAD(q);
491 3881472 : __SSE_DATATYPE y1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
492 3881472 : __SSE_DATATYPE q2 = _SSE_LOAD(&q[offset]);
493 3881472 : __SSE_DATATYPE y2 = _SSE_ADD(q2, _SSE_MUL(x2, h1));
494 3881472 : __SSE_DATATYPE q3 = _SSE_LOAD(&q[2*offset]);
495 3881472 : __SSE_DATATYPE y3 = _SSE_ADD(q3, _SSE_MUL(x3, h1));
496 3881472 : __SSE_DATATYPE q4 = _SSE_LOAD(&q[3*offset]);
497 3881472 : __SSE_DATATYPE y4 = _SSE_ADD(q4, _SSE_MUL(x4, h1));
498 3881472 : __SSE_DATATYPE q5 = _SSE_LOAD(&q[4*offset]);
499 3881472 : __SSE_DATATYPE y5 = _SSE_ADD(q5, _SSE_MUL(x5, h1));
500 3881472 : __SSE_DATATYPE q6 = _SSE_LOAD(&q[5*offset]);
501 3881472 : __SSE_DATATYPE y6 = _SSE_ADD(q6, _SSE_MUL(x6, h1));
502 122266368 : for(i = 2; i < nb; i++)
503 : {
504 : #ifdef HAVE_SSE_INTRINSICS
505 : #ifdef DOUBLE_PRECISION_REAL
506 202653696 : h1 = _mm_set1_pd(hh[i-1]);
507 202653696 : h2 = _mm_set1_pd(hh[ldh+i]);
508 : #endif
509 : #ifdef SINGLE_PRECISION_REAL
510 37997568 : h1 = _mm_set1_ps(hh[i-1]);
511 37997568 : h2 = _mm_set1_ps(hh[ldh+i]);
512 : #endif
513 : #endif
514 : #ifdef HAVE_SPARC64_SSE
515 : #ifdef DOUBLE_PRECISION_REAL
516 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
517 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
518 : #endif
519 : #ifdef SINGLE_PRECISION_REAL
520 : h1 = _mm_set_ps(hh[i-1], hh[i-1]);
521 : h2 = _mm_set_ps(hh[ldh+i], hh[ldh+i]);
522 : #endif
523 : #endif
524 :
525 240651264 : q1 = _SSE_LOAD(&q[i*ldq]);
526 240651264 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
527 240651264 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
528 240651264 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
529 240651264 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
530 240651264 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
531 240651264 : q3 = _SSE_LOAD(&q[(i*ldq)+2*offset]);
532 240651264 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
533 240651264 : y3 = _SSE_ADD(y3, _SSE_MUL(q3,h2));
534 240651264 : q4 = _SSE_LOAD(&q[(i*ldq)+3*offset]);
535 240651264 : x4 = _SSE_ADD(x4, _SSE_MUL(q4,h1));
536 240651264 : y4 = _SSE_ADD(y4, _SSE_MUL(q4,h2));
537 240651264 : q5 = _SSE_LOAD(&q[(i*ldq)+4*offset]);
538 240651264 : x5 = _SSE_ADD(x5, _SSE_MUL(q5,h1));
539 240651264 : y5 = _SSE_ADD(y5, _SSE_MUL(q5,h2));
540 240651264 : q6 = _SSE_LOAD(&q[(i*ldq)+5*offset]);
541 240651264 : x6 = _SSE_ADD(x6, _SSE_MUL(q6,h1));
542 240651264 : y6 = _SSE_ADD(y6, _SSE_MUL(q6,h2));
543 : }
544 : #ifdef HAVE_SSE_INTRINSICS
545 : #ifdef DOUBLE_PRECISION_REAL
546 3268608 : h1 = _mm_set1_pd(hh[nb-1]);
547 : #endif
548 : #ifdef SINGLE_PRECISION_REAL
549 612864 : h1 = _mm_set1_ps(hh[nb-1]);
550 : #endif
551 : #endif
552 : #ifdef HAVE_SPARC64_SSE
553 : #ifdef DOUBLE_PRECISION_REAL
554 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
555 : #endif
556 : #ifdef SINGLE_PRECISION_REAL
557 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
558 : #endif
559 : #endif
560 :
561 3881472 : q1 = _SSE_LOAD(&q[nb*ldq]);
562 3881472 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
563 3881472 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
564 3881472 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
565 3881472 : q3 = _SSE_LOAD(&q[(nb*ldq)+2*offset]);
566 3881472 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
567 3881472 : q4 = _SSE_LOAD(&q[(nb*ldq)+3*offset]);
568 3881472 : x4 = _SSE_ADD(x4, _SSE_MUL(q4,h1));
569 3881472 : q5 = _SSE_LOAD(&q[(nb*ldq)+4*offset]);
570 3881472 : x5 = _SSE_ADD(x5, _SSE_MUL(q5,h1));
571 3881472 : q6 = _SSE_LOAD(&q[(nb*ldq)+5*offset]);
572 3881472 : x6 = _SSE_ADD(x6, _SSE_MUL(q6,h1));
573 : /////////////////////////////////////////////////////
574 : // Rank-2 update of Q [12 x nb+1]
575 : /////////////////////////////////////////////////////
576 : #ifdef HAVE_SSE_INTRINSICS
577 : #ifdef DOUBLE_PRECISION_REAL
578 3268608 : __SSE_DATATYPE tau1 = _mm_set1_pd(hh[0]);
579 3268608 : __SSE_DATATYPE tau2 = _mm_set1_pd(hh[ldh]);
580 1634304 : __SSE_DATATYPE vs = _mm_set1_pd(s);
581 : #endif
582 : #ifdef SINGLE_PRECISION_REAL
583 612864 : __SSE_DATATYPE tau1 = _mm_set1_ps(hh[0]);
584 612864 : __SSE_DATATYPE tau2 = _mm_set1_ps(hh[ldh]);
585 306432 : __SSE_DATATYPE vs = _mm_set1_ps(s);
586 : #endif
587 : #endif
588 : #ifdef HAVE_SPARC64_SSE
589 : #ifdef DOUBLE_PRECISION_REAL
590 : __SSE_DATATYPE tau1 = _mm_set_pd(hh[0], hh[0]);
591 : __SSE_DATATYPE tau2 = _mm_set_pd(hh[ldh], hh[ldh]);
592 : __SSE_DATATYPE vs = _mm_set_pd(s, s);
593 : #endif
594 : #ifdef SINGLE_PRECISION_REAL
595 : __SSE_DATATYPE tau1 = _mm_set_ps(hh[0], hh[0]);
596 : __SSE_DATATYPE tau2 = _mm_set_ps(hh[ldh], hh[ldh]);
597 : __SSE_DATATYPE vs = _mm_set_ps(s, s);
598 : #endif
599 : #endif
600 :
601 : #ifdef HAVE_SSE_INTRINSICS
602 1940736 : h1 = _SSE_XOR(tau1, sign);
603 : #endif
604 : #ifdef HAVE_SPARC64_SSE
605 : h1 = _fjsp_neg_v2r8(tau1);
606 : #endif
607 1940736 : x1 = _SSE_MUL(x1, h1);
608 1940736 : x2 = _SSE_MUL(x2, h1);
609 1940736 : x3 = _SSE_MUL(x3, h1);
610 1940736 : x4 = _SSE_MUL(x4, h1);
611 1940736 : x5 = _SSE_MUL(x5, h1);
612 1940736 : x6 = _SSE_MUL(x6, h1);
613 : #ifdef HAVE_SSE_INTRINSICS
614 1940736 : h1 = _SSE_XOR(tau2, sign);
615 : #endif
616 : #ifdef HAVE_SPARC64_SSE
617 : h1 = _fjsp_neg_v2r8(tau2);
618 : #endif
619 1940736 : h2 = _SSE_MUL(h1, vs);
620 :
621 5822208 : y1 = _SSE_ADD(_SSE_MUL(y1,h1), _SSE_MUL(x1,h2));
622 5822208 : y2 = _SSE_ADD(_SSE_MUL(y2,h1), _SSE_MUL(x2,h2));
623 5822208 : y3 = _SSE_ADD(_SSE_MUL(y3,h1), _SSE_MUL(x3,h2));
624 5822208 : y4 = _SSE_ADD(_SSE_MUL(y4,h1), _SSE_MUL(x4,h2));
625 5822208 : y5 = _SSE_ADD(_SSE_MUL(y5,h1), _SSE_MUL(x5,h2));
626 5822208 : y6 = _SSE_ADD(_SSE_MUL(y6,h1), _SSE_MUL(x6,h2));
627 1940736 : q1 = _SSE_LOAD(q);
628 1940736 : q1 = _SSE_ADD(q1, y1);
629 : _SSE_STORE(q,q1);
630 3881472 : q2 = _SSE_LOAD(&q[offset]);
631 1940736 : q2 = _SSE_ADD(q2, y2);
632 1940736 : _SSE_STORE(&q[offset],q2);
633 3881472 : q3 = _SSE_LOAD(&q[2*offset]);
634 1940736 : q3 = _SSE_ADD(q3, y3);
635 1940736 : _SSE_STORE(&q[2*offset],q3);
636 3881472 : q4 = _SSE_LOAD(&q[3*offset]);
637 1940736 : q4 = _SSE_ADD(q4, y4);
638 1940736 : _SSE_STORE(&q[3*offset],q4);
639 3881472 : q5 = _SSE_LOAD(&q[4*offset]);
640 1940736 : q5 = _SSE_ADD(q5, y5);
641 1940736 : _SSE_STORE(&q[4*offset],q5);
642 3881472 : q6 = _SSE_LOAD(&q[5*offset]);
643 1940736 : q6 = _SSE_ADD(q6, y6);
644 1940736 : _SSE_STORE(&q[5*offset],q6);
645 :
646 : #ifdef HAVE_SSE_INTRINSICS
647 : #ifdef DOUBLE_PRECISION_REAL
648 3268608 : h2 = _mm_set1_pd(hh[ldh+1]);
649 : #endif
650 : #ifdef SINGLE_PRECISION_REAL
651 612864 : h2 = _mm_set1_ps(hh[ldh+1]);
652 : #endif
653 : #endif
654 :
655 : #ifdef HAVE_SPARC64_SSE
656 : #ifdef DOUBLE_PRECISION_REAL
657 : h2 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
658 : #endif
659 : #ifdef SINGLE_PRECISION_REAL
660 : h2 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
661 : #endif
662 : #endif
663 :
664 3881472 : q1 = _SSE_LOAD(&q[ldq]);
665 5822208 : q1 = _SSE_ADD(q1, _SSE_ADD(x1, _SSE_MUL(y1, h2)));
666 1940736 : _SSE_STORE(&q[ldq],q1);
667 3881472 : q2 = _SSE_LOAD(&q[ldq+offset]);
668 5822208 : q2 = _SSE_ADD(q2, _SSE_ADD(x2, _SSE_MUL(y2, h2)));
669 1940736 : _SSE_STORE(&q[ldq+offset],q2);
670 3881472 : q3 = _SSE_LOAD(&q[ldq+2*offset]);
671 5822208 : q3 = _SSE_ADD(q3, _SSE_ADD(x3, _SSE_MUL(y3, h2)));
672 1940736 : _SSE_STORE(&q[ldq+2*offset],q3);
673 3881472 : q4 = _SSE_LOAD(&q[ldq+3*offset]);
674 5822208 : q4 = _SSE_ADD(q4, _SSE_ADD(x4, _SSE_MUL(y4, h2)));
675 1940736 : _SSE_STORE(&q[ldq+3*offset],q4);
676 3881472 : q5 = _SSE_LOAD(&q[ldq+4*offset]);
677 5822208 : q5 = _SSE_ADD(q5, _SSE_ADD(x5, _SSE_MUL(y5, h2)));
678 1940736 : _SSE_STORE(&q[ldq+4*offset],q5);
679 3881472 : q6 = _SSE_LOAD(&q[ldq+5*offset]);
680 5822208 : q6 = _SSE_ADD(q6, _SSE_ADD(x6, _SSE_MUL(y6, h2)));
681 1940736 : _SSE_STORE(&q[ldq+5*offset],q6);
682 :
683 122266368 : for (i = 2; i < nb; i++)
684 : {
685 : #ifdef HAVE_SSE_INTRINSICS
686 : #ifdef DOUBLE_PRECISION_REAL
687 202653696 : h1 = _mm_set1_pd(hh[i-1]);
688 202653696 : h2 = _mm_set1_pd(hh[ldh+i]);
689 : #endif
690 : #ifdef SINGLE_PRECISION_REAL
691 37997568 : h1 = _mm_set1_ps(hh[i-1]);
692 37997568 : h2 = _mm_set1_ps(hh[ldh+i]);
693 : #endif
694 : #endif
695 : #ifdef HAVE_SPARC64_SSE
696 : #ifdef DOUBLE_PRECISION_REAL
697 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
698 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
699 : #endif
700 : #ifdef SINGLE_PRECISION_REAL
701 : h1 = _mm_set_ps(hh[i-1], hh[i-1]);
702 : h2 = _mm_set_ps(hh[ldh+i], hh[ldh+i]);
703 : #endif
704 : #endif
705 :
706 240651264 : q1 = _SSE_LOAD(&q[i*ldq]);
707 481302528 : q1 = _SSE_ADD(q1, _SSE_ADD(_SSE_MUL(x1,h1), _SSE_MUL(y1, h2)));
708 120325632 : _SSE_STORE(&q[i*ldq],q1);
709 240651264 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
710 481302528 : q2 = _SSE_ADD(q2, _SSE_ADD(_SSE_MUL(x2,h1), _SSE_MUL(y2, h2)));
711 120325632 : _SSE_STORE(&q[(i*ldq)+offset],q2);
712 240651264 : q3 = _SSE_LOAD(&q[(i*ldq)+2*offset]);
713 481302528 : q3 = _SSE_ADD(q3, _SSE_ADD(_SSE_MUL(x3,h1), _SSE_MUL(y3, h2)));
714 120325632 : _SSE_STORE(&q[(i*ldq)+2*offset],q3);
715 240651264 : q4 = _SSE_LOAD(&q[(i*ldq)+3*offset]);
716 481302528 : q4 = _SSE_ADD(q4, _SSE_ADD(_SSE_MUL(x4,h1), _SSE_MUL(y4, h2)));
717 120325632 : _SSE_STORE(&q[(i*ldq)+3*offset],q4);
718 240651264 : q5 = _SSE_LOAD(&q[(i*ldq)+4*offset]);
719 481302528 : q5 = _SSE_ADD(q5, _SSE_ADD(_SSE_MUL(x5,h1), _SSE_MUL(y5, h2)));
720 120325632 : _SSE_STORE(&q[(i*ldq)+4*offset],q5);
721 240651264 : q6 = _SSE_LOAD(&q[(i*ldq)+5*offset]);
722 481302528 : q6 = _SSE_ADD(q6, _SSE_ADD(_SSE_MUL(x6,h1), _SSE_MUL(y6, h2)));
723 120325632 : _SSE_STORE(&q[(i*ldq)+5*offset],q6);
724 : }
725 : #ifdef HAVE_SSE_INTRINSICS
726 : #ifdef DOUBLE_PRECISION_REAL
727 3268608 : h1 = _mm_set1_pd(hh[nb-1]);
728 : #endif
729 : #ifdef SINGLE_PRECISION_REAL
730 612864 : h1 = _mm_set1_ps(hh[nb-1]);
731 : #endif
732 : #endif
733 : #ifdef HAVE_SPARC64_SSE
734 : #ifdef DOUBLE_PRECISION_REAL
735 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
736 : #endif
737 : #ifdef SINGLE_PRECISION_REAL
738 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
739 : #endif
740 : #endif
741 :
742 :
743 3881472 : q1 = _SSE_LOAD(&q[nb*ldq]);
744 3881472 : q1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
745 1940736 : _SSE_STORE(&q[nb*ldq],q1);
746 3881472 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
747 3881472 : q2 = _SSE_ADD(q2, _SSE_MUL(x2, h1));
748 1940736 : _SSE_STORE(&q[(nb*ldq)+offset],q2);
749 3881472 : q3 = _SSE_LOAD(&q[(nb*ldq)+2*offset]);
750 3881472 : q3 = _SSE_ADD(q3, _SSE_MUL(x3, h1));
751 1940736 : _SSE_STORE(&q[(nb*ldq)+2*offset],q3);
752 3881472 : q4 = _SSE_LOAD(&q[(nb*ldq)+3*offset]);
753 3881472 : q4 = _SSE_ADD(q4, _SSE_MUL(x4, h1));
754 1940736 : _SSE_STORE(&q[(nb*ldq)+3*offset],q4);
755 3881472 : q5 = _SSE_LOAD(&q[(nb*ldq)+4*offset]);
756 3881472 : q5 = _SSE_ADD(q5, _SSE_MUL(x5, h1));
757 1940736 : _SSE_STORE(&q[(nb*ldq)+4*offset],q5);
758 3881472 : q6 = _SSE_LOAD(&q[(nb*ldq)+5*offset]);
759 3881472 : q6 = _SSE_ADD(q6, _SSE_MUL(x6, h1));
760 1940736 : _SSE_STORE(&q[(nb*ldq)+5*offset],q6);
761 : }
762 :
763 :
764 :
765 : /**
766 : * Unrolled kernel that computes
767 : #ifdef DOUBLE_PRECISION_REAL
768 : * 10 rows of Q simultaneously, a
769 : #endif
770 : #ifdef SINGLE_PRECISION_REAL
771 : * 20 rows of Q simultaneously, a
772 : #endif
773 : * matrix Vector product with two householder
774 : * vectors + a rank 2 update is performed
775 : */
776 : #ifdef HAVE_SSE_INTRINSICS
777 : #ifdef DOUBLE_PRECISION_REAL
778 : __forceinline void hh_trafo_kernel_10_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
779 : #endif
780 : #ifdef SINGLE_PRECISION_REAL
781 : __forceinline void hh_trafo_kernel_20_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
782 : #endif
783 : #endif
784 : #ifdef HAVE_SPARC64_SSE
785 : #ifdef DOUBLE_PRECISION_REAL
786 : __forceinline void hh_trafo_kernel_10_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
787 : #endif
788 : #ifdef SINGLE_PRECISION_REAL
789 : __forceinline void hh_trafo_kernel_20_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
790 : #endif
791 :
792 : #endif
793 : {
794 : /////////////////////////////////////////////////////
795 : // Matrix Vector Multiplication, Q [12 x nb+1] * hh
796 : // hh contains two householder vectors, with offset 1
797 : /////////////////////////////////////////////////////
798 : int i;
799 : #ifdef HAVE_SSE_INTRINSICS
800 : // Needed bit mask for floating point sign flip
801 : #ifdef DOUBLE_PRECISION_REAL
802 0 : __SSE_DATATYPE sign = (__SSE_DATATYPE)_mm_set1_epi64x(0x8000000000000000LL);
803 : #endif
804 : #ifdef SINGLE_PRECISION_REAL
805 0 : __SSE_DATATYPE sign = _mm_castsi128_ps(_mm_set_epi32(0x80000000, 0x80000000, 0x80000000, 0x80000000));
806 : #endif
807 : #endif
808 :
809 0 : __SSE_DATATYPE x1 = _SSE_LOAD(&q[ldq]);
810 0 : __SSE_DATATYPE x2 = _SSE_LOAD(&q[ldq+offset]);
811 0 : __SSE_DATATYPE x3 = _SSE_LOAD(&q[ldq+2*offset]);
812 0 : __SSE_DATATYPE x4 = _SSE_LOAD(&q[ldq+3*offset]);
813 0 : __SSE_DATATYPE x5 = _SSE_LOAD(&q[ldq+4*offset]);
814 :
815 : #ifdef HAVE_SSE_INTRINSICS
816 : #ifdef DOUBLE_PRECISION_REAL
817 0 : __SSE_DATATYPE h1 = _mm_set1_pd(hh[ldh+1]);
818 : #endif
819 : #ifdef SINGLE_PRECISION_REAL
820 0 : __SSE_DATATYPE h1 = _mm_set1_ps(hh[ldh+1]);
821 : #endif
822 : #endif
823 : #ifdef HAVE_SPARC64_SSE
824 : #ifdef DOUBLE_PRECISION_REAL
825 : __SSE_DATATYPE h1 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
826 : #endif
827 : #ifdef SINGLE_PRECISION_REAL
828 : __SSE_DATATYPE h1 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
829 : #endif
830 : #endif
831 :
832 : __SSE_DATATYPE h2;
833 :
834 0 : __SSE_DATATYPE q1 = _SSE_LOAD(q);
835 0 : __SSE_DATATYPE y1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
836 0 : __SSE_DATATYPE q2 = _SSE_LOAD(&q[offset]);
837 0 : __SSE_DATATYPE y2 = _SSE_ADD(q2, _SSE_MUL(x2, h1));
838 0 : __SSE_DATATYPE q3 = _SSE_LOAD(&q[2*offset]);
839 0 : __SSE_DATATYPE y3 = _SSE_ADD(q3, _SSE_MUL(x3, h1));
840 0 : __SSE_DATATYPE q4 = _SSE_LOAD(&q[3*offset]);
841 0 : __SSE_DATATYPE y4 = _SSE_ADD(q4, _SSE_MUL(x4, h1));
842 0 : __SSE_DATATYPE q5 = _SSE_LOAD(&q[4*offset]);
843 0 : __SSE_DATATYPE y5 = _SSE_ADD(q5, _SSE_MUL(x5, h1));
844 0 : for(i = 2; i < nb; i++)
845 : {
846 : #ifdef HAVE_SSE_INTRINSICS
847 : #ifdef DOUBLE_PRECISION_REAL
848 0 : h1 = _mm_set1_pd(hh[i-1]);
849 0 : h2 = _mm_set1_pd(hh[ldh+i]);
850 : #endif
851 : #ifdef SINGLE_PRECISION_REAL
852 0 : h1 = _mm_set1_ps(hh[i-1]);
853 0 : h2 = _mm_set1_ps(hh[ldh+i]);
854 : #endif
855 : #endif
856 : #ifdef HAVE_SPARC64_SSE
857 : #ifdef DOUBLE_PRECISION_REAL
858 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
859 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
860 : #endif
861 : #ifdef SINGLE_PRECISION_REAL
862 : h1 = _mm_set_ps(hh[i-1], hh[i-1]);
863 : h2 = _mm_set_ps(hh[ldh+i], hh[ldh+i]);
864 : #endif
865 : #endif
866 :
867 :
868 0 : q1 = _SSE_LOAD(&q[i*ldq]);
869 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
870 0 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
871 0 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
872 0 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
873 0 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
874 0 : q3 = _SSE_LOAD(&q[(i*ldq)+2*offset]);
875 0 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
876 0 : y3 = _SSE_ADD(y3, _SSE_MUL(q3,h2));
877 0 : q4 = _SSE_LOAD(&q[(i*ldq)+3*offset]);
878 0 : x4 = _SSE_ADD(x4, _SSE_MUL(q4,h1));
879 0 : y4 = _SSE_ADD(y4, _SSE_MUL(q4,h2));
880 0 : q5 = _SSE_LOAD(&q[(i*ldq)+4*offset]);
881 0 : x5 = _SSE_ADD(x5, _SSE_MUL(q5,h1));
882 0 : y5 = _SSE_ADD(y5, _SSE_MUL(q5,h2));
883 : }
884 :
885 : #ifdef HAVE_SSE_INTRINSICS
886 : #ifdef DOUBLE_PRECISION_REAL
887 0 : h1 = _mm_set1_pd(hh[nb-1]);
888 : #endif
889 : #ifdef SINGLE_PRECISION_REAL
890 0 : h1 = _mm_set1_ps(hh[nb-1]);
891 : #endif
892 : #endif
893 : #ifdef HAVE_SPARC64_SSE
894 : #ifdef DOUBLE_PRECISION_REAL
895 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
896 : #endif
897 : #ifdef SINGLE_PRECISION_REAL
898 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
899 : #endif
900 : #endif
901 :
902 0 : q1 = _SSE_LOAD(&q[nb*ldq]);
903 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
904 0 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
905 0 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
906 0 : q3 = _SSE_LOAD(&q[(nb*ldq)+2*offset]);
907 0 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
908 0 : q4 = _SSE_LOAD(&q[(nb*ldq)+3*offset]);
909 0 : x4 = _SSE_ADD(x4, _SSE_MUL(q4,h1));
910 0 : q5 = _SSE_LOAD(&q[(nb*ldq)+4*offset]);
911 0 : x5 = _SSE_ADD(x5, _SSE_MUL(q5,h1));
912 : /////////////////////////////////////////////////////
913 : // Rank-2 update of Q [12 x nb+1]
914 : /////////////////////////////////////////////////////
915 : #ifdef HAVE_SSE_INTRINSICS
916 : #ifdef DOUBLE_PRECISION_REAL
917 0 : __SSE_DATATYPE tau1 = _mm_set1_pd(hh[0]);
918 0 : __SSE_DATATYPE tau2 = _mm_set1_pd(hh[ldh]);
919 0 : __SSE_DATATYPE vs = _mm_set1_pd(s);
920 : #endif
921 : #ifdef SINGLE_PRECISION_REAL
922 0 : __SSE_DATATYPE tau1 = _mm_set1_ps(hh[0]);
923 0 : __SSE_DATATYPE tau2 = _mm_set1_ps(hh[ldh]);
924 0 : __SSE_DATATYPE vs = _mm_set1_ps(s);
925 :
926 : #endif
927 : #endif
928 : #ifdef HAVE_SPARC64_SSE
929 : #ifdef DOUBLE_PRECISION_REAL
930 : __SSE_DATATYPE tau1 = _mm_set_pd(hh[0], hh[0]);
931 : __SSE_DATATYPE tau2 = _mm_set_pd(hh[ldh], hh[ldh]);
932 : __SSE_DATATYPE vs = _mm_set_pd(s, s);
933 : #endif
934 : #ifdef SINGLE_PRECISION_REAL
935 : __SSE_DATATYPE tau1 = _mm_set_ps(hh[0], hh[0]);
936 : __SSE_DATATYPE tau2 = _mm_set_ps(hh[ldh], hh[ldh]);
937 : __SSE_DATATYPE vs = _mm_set_ps(s, s);
938 :
939 : #endif
940 : #endif
941 :
942 : #ifdef HAVE_SSE_INTRINSICS
943 0 : h1 = _SSE_XOR(tau1, sign);
944 : #endif
945 : #ifdef HAVE_SPARC64_SSE
946 : h1 = _fjsp_neg_v2r8(tau1);
947 : #endif
948 0 : x1 = _SSE_MUL(x1, h1);
949 0 : x2 = _SSE_MUL(x2, h1);
950 0 : x3 = _SSE_MUL(x3, h1);
951 0 : x4 = _SSE_MUL(x4, h1);
952 0 : x5 = _SSE_MUL(x5, h1);
953 : #ifdef HAVE_SSE_INTRINSICS
954 0 : h1 = _SSE_XOR(tau2, sign);
955 : #endif
956 : #ifdef HAVE_SPARC64_SSE
957 : h1 = _fjsp_neg_v2r8(tau2);
958 : #endif
959 0 : h2 = _SSE_MUL(h1, vs);
960 :
961 0 : y1 = _SSE_ADD(_SSE_MUL(y1,h1), _SSE_MUL(x1,h2));
962 0 : y2 = _SSE_ADD(_SSE_MUL(y2,h1), _SSE_MUL(x2,h2));
963 0 : y3 = _SSE_ADD(_SSE_MUL(y3,h1), _SSE_MUL(x3,h2));
964 0 : y4 = _SSE_ADD(_SSE_MUL(y4,h1), _SSE_MUL(x4,h2));
965 0 : y5 = _SSE_ADD(_SSE_MUL(y5,h1), _SSE_MUL(x5,h2));
966 0 : q1 = _SSE_LOAD(q);
967 0 : q1 = _SSE_ADD(q1, y1);
968 : _SSE_STORE(q,q1);
969 0 : q2 = _SSE_LOAD(&q[offset]);
970 0 : q2 = _SSE_ADD(q2, y2);
971 0 : _SSE_STORE(&q[offset],q2);
972 0 : q3 = _SSE_LOAD(&q[2*offset]);
973 0 : q3 = _SSE_ADD(q3, y3);
974 0 : _SSE_STORE(&q[2*offset],q3);
975 0 : q4 = _SSE_LOAD(&q[3*offset]);
976 0 : q4 = _SSE_ADD(q4, y4);
977 0 : _SSE_STORE(&q[3*offset],q4);
978 0 : q5 = _SSE_LOAD(&q[4*offset]);
979 0 : q5 = _SSE_ADD(q5, y5);
980 0 : _SSE_STORE(&q[4*offset],q5);
981 :
982 : #ifdef HAVE_SSE_INTRINSICS
983 : #ifdef DOUBLE_PRECISION_REAL
984 0 : h2 = _mm_set1_pd(hh[ldh+1]);
985 : #endif
986 : #ifdef SINGLE_PRECISION_REAL
987 0 : h2 = _mm_set1_ps(hh[ldh+1]);
988 : #endif
989 : #endif
990 : #ifdef HAVE_SPARC64_SSE
991 : #ifdef DOUBLE_PRECISION_REAL
992 : h2 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
993 : #endif
994 : #ifdef SINGLE_PRECISION_REAL
995 : h2 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
996 : #endif
997 : #endif
998 :
999 0 : q1 = _SSE_LOAD(&q[ldq]);
1000 0 : q1 = _SSE_ADD(q1, _SSE_ADD(x1, _SSE_MUL(y1, h2)));
1001 0 : _SSE_STORE(&q[ldq],q1);
1002 0 : q2 = _SSE_LOAD(&q[ldq+offset]);
1003 0 : q2 = _SSE_ADD(q2, _SSE_ADD(x2, _SSE_MUL(y2, h2)));
1004 0 : _SSE_STORE(&q[ldq+offset],q2);
1005 0 : q3 = _SSE_LOAD(&q[ldq+2*offset]);
1006 0 : q3 = _SSE_ADD(q3, _SSE_ADD(x3, _SSE_MUL(y3, h2)));
1007 0 : _SSE_STORE(&q[ldq+2*offset],q3);
1008 0 : q4 = _SSE_LOAD(&q[ldq+3*offset]);
1009 0 : q4 = _SSE_ADD(q4, _SSE_ADD(x4, _SSE_MUL(y4, h2)));
1010 0 : _SSE_STORE(&q[ldq+3*offset],q4);
1011 0 : q5 = _SSE_LOAD(&q[ldq+4*offset]);
1012 0 : q5 = _SSE_ADD(q5, _SSE_ADD(x5, _SSE_MUL(y5, h2)));
1013 0 : _SSE_STORE(&q[ldq+4*offset],q5);
1014 :
1015 0 : for (i = 2; i < nb; i++)
1016 : {
1017 : #ifdef HAVE_SSE_INTRINSICS
1018 : #ifdef DOUBLE_PRECISION_REAL
1019 0 : h1 = _mm_set1_pd(hh[i-1]);
1020 0 : h2 = _mm_set1_pd(hh[ldh+i]);
1021 : #endif
1022 : #ifdef SINGLE_PRECISION_REAL
1023 0 : h1 = _mm_set1_ps(hh[i-1]);
1024 0 : h2 = _mm_set1_ps(hh[ldh+i]);
1025 : #endif
1026 : #endif
1027 : #ifdef HAVE_SPARC64_SSE
1028 : #ifdef DOUBLE_PRECISION_REAL
1029 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
1030 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
1031 : #endif
1032 : #ifdef SINGLE_PRECISION_REAL
1033 : h1 = _mm_set_ps(hh[i-1], hh[i-1]);
1034 : h2 = _mm_set_ps(hh[ldh+i], hh[ldh+i]);
1035 : #endif
1036 : #endif
1037 :
1038 0 : q1 = _SSE_LOAD(&q[i*ldq]);
1039 0 : q1 = _SSE_ADD(q1, _SSE_ADD(_SSE_MUL(x1,h1), _SSE_MUL(y1, h2)));
1040 0 : _SSE_STORE(&q[i*ldq],q1);
1041 0 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
1042 0 : q2 = _SSE_ADD(q2, _SSE_ADD(_SSE_MUL(x2,h1), _SSE_MUL(y2, h2)));
1043 0 : _SSE_STORE(&q[(i*ldq)+offset],q2);
1044 0 : q3 = _SSE_LOAD(&q[(i*ldq)+2*offset]);
1045 0 : q3 = _SSE_ADD(q3, _SSE_ADD(_SSE_MUL(x3,h1), _SSE_MUL(y3, h2)));
1046 0 : _SSE_STORE(&q[(i*ldq)+2*offset],q3);
1047 0 : q4 = _SSE_LOAD(&q[(i*ldq)+3*offset]);
1048 0 : q4 = _SSE_ADD(q4, _SSE_ADD(_SSE_MUL(x4,h1), _SSE_MUL(y4, h2)));
1049 0 : _SSE_STORE(&q[(i*ldq)+3*offset],q4);
1050 0 : q5 = _SSE_LOAD(&q[(i*ldq)+4*offset]);
1051 0 : q5 = _SSE_ADD(q5, _SSE_ADD(_SSE_MUL(x5,h1), _SSE_MUL(y5, h2)));
1052 0 : _SSE_STORE(&q[(i*ldq)+4*offset],q5);
1053 : }
1054 : #ifdef HAVE_SSE_INTRINSICS
1055 : #ifdef DOUBLE_PRECISION_REAL
1056 0 : h1 = _mm_set1_pd(hh[nb-1]);
1057 : #endif
1058 : #ifdef SINGLE_PRECISION_REAL
1059 0 : h1 = _mm_set1_ps(hh[nb-1]);
1060 : #endif
1061 : #endif
1062 : #ifdef HAVE_SPARC64_SSE
1063 : #ifdef DOUBLE_PRECISION_REAL
1064 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
1065 : #endif
1066 : #ifdef SINGLE_PRECISION_REAL
1067 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
1068 : #endif
1069 : #endif
1070 :
1071 0 : q1 = _SSE_LOAD(&q[nb*ldq]);
1072 0 : q1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
1073 0 : _SSE_STORE(&q[nb*ldq],q1);
1074 0 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
1075 0 : q2 = _SSE_ADD(q2, _SSE_MUL(x2, h1));
1076 0 : _SSE_STORE(&q[(nb*ldq)+offset],q2);
1077 0 : q3 = _SSE_LOAD(&q[(nb*ldq)+2*offset]);
1078 0 : q3 = _SSE_ADD(q3, _SSE_MUL(x3, h1));
1079 0 : _SSE_STORE(&q[(nb*ldq)+2*offset],q3);
1080 0 : q4 = _SSE_LOAD(&q[(nb*ldq)+3*offset]);
1081 0 : q4 = _SSE_ADD(q4, _SSE_MUL(x4, h1));
1082 0 : _SSE_STORE(&q[(nb*ldq)+3*offset],q4);
1083 0 : q5 = _SSE_LOAD(&q[(nb*ldq)+4*offset]);
1084 0 : q5 = _SSE_ADD(q5, _SSE_MUL(x5, h1));
1085 0 : _SSE_STORE(&q[(nb*ldq)+4*offset],q5);
1086 : }
1087 :
1088 : /**
1089 : * Unrolled kernel that computes
1090 : #ifdef DOUBLE_PRECISION_REAL
1091 : * 8 rows of Q simultaneously, a
1092 : #endif
1093 : #ifdef SINGLE_PRECISION_REAL
1094 : * 16 rows of Q simultaneously, a
1095 : #endif
1096 : * matrix Vector product with two householder
1097 : * vectors + a rank 2 update is performed
1098 : */
1099 : #ifdef HAVE_SSE_INTRINSICS
1100 : #ifdef DOUBLE_PRECISION_REAL
1101 : __forceinline void hh_trafo_kernel_8_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
1102 : #endif
1103 : #ifdef SINGLE_PRECISION_REAL
1104 : __forceinline void hh_trafo_kernel_16_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
1105 : #endif
1106 : #endif
1107 : #ifdef HAVE_SPARC64_SSE
1108 : #ifdef DOUBLE_PRECISION_REAL
1109 : __forceinline void hh_trafo_kernel_8_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
1110 : #endif
1111 : #ifdef SINGLE_PRECISION_REAL
1112 : __forceinline void hh_trafo_kernel_16_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
1113 : #endif
1114 :
1115 : #endif
1116 : {
1117 : /////////////////////////////////////////////////////
1118 : // Matrix Vector Multiplication, Q [12 x nb+1] * hh
1119 : // hh contains two householder vectors, with offset 1
1120 : /////////////////////////////////////////////////////
1121 : int i;
1122 :
1123 : #ifdef HAVE_SSE_INTRINSICS
1124 : // Needed bit mask for floating point sign flip
1125 : #ifdef DOUBLE_PRECISION_REAL
1126 0 : __SSE_DATATYPE sign = (__SSE_DATATYPE)_mm_set1_epi64x(0x8000000000000000LL);
1127 : #endif
1128 : #ifdef SINGLE_PRECISION_REAL
1129 0 : __SSE_DATATYPE sign = _mm_castsi128_ps(_mm_set_epi32(0x80000000, 0x80000000, 0x80000000, 0x80000000));
1130 : #endif
1131 : #endif
1132 0 : __SSE_DATATYPE x1 = _SSE_LOAD(&q[ldq]);
1133 0 : __SSE_DATATYPE x2 = _SSE_LOAD(&q[ldq+offset]);
1134 0 : __SSE_DATATYPE x3 = _SSE_LOAD(&q[ldq+2*offset]);
1135 0 : __SSE_DATATYPE x4 = _SSE_LOAD(&q[ldq+3*offset]);
1136 :
1137 : #ifdef HAVE_SSE_INTRINSICS
1138 : #ifdef DOUBLE_PRECISION_REAL
1139 0 : __SSE_DATATYPE h1 = _mm_set1_pd(hh[ldh+1]);
1140 : #endif
1141 : #ifdef SINGLE_PRECISION_REAL
1142 0 : __SSE_DATATYPE h1 = _mm_set1_ps(hh[ldh+1]);
1143 : #endif
1144 : #endif
1145 : #ifdef HAVE_SPARC64_SSE
1146 : #ifdef DOUBLE_PRECISION_REAL
1147 : __SSE_DATATYPE h1 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1148 : #endif
1149 : #ifdef SINGLE_PRECISION_REAL
1150 : __SSE_DATATYPE h1 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
1151 : #endif
1152 : #endif
1153 :
1154 : __SSE_DATATYPE h2;
1155 :
1156 0 : __SSE_DATATYPE q1 = _SSE_LOAD(q);
1157 0 : __SSE_DATATYPE y1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
1158 0 : __SSE_DATATYPE q2 = _SSE_LOAD(&q[offset]);
1159 0 : __SSE_DATATYPE y2 = _SSE_ADD(q2, _SSE_MUL(x2, h1));
1160 0 : __SSE_DATATYPE q3 = _SSE_LOAD(&q[2*offset]);
1161 0 : __SSE_DATATYPE y3 = _SSE_ADD(q3, _SSE_MUL(x3, h1));
1162 0 : __SSE_DATATYPE q4 = _SSE_LOAD(&q[3*offset]);
1163 0 : __SSE_DATATYPE y4 = _SSE_ADD(q4, _SSE_MUL(x4, h1));
1164 0 : for(i = 2; i < nb; i++)
1165 : {
1166 : #ifdef HAVE_SSE_INTRINSICS
1167 : #ifdef DOUBLE_PRECISION_REAL
1168 0 : h1 = _mm_set1_pd(hh[i-1]);
1169 0 : h2 = _mm_set1_pd(hh[ldh+i]);
1170 : #endif
1171 : #ifdef SINGLE_PRECISION_REAL
1172 0 : h1 = _mm_set1_ps(hh[i-1]);
1173 0 : h2 = _mm_set1_ps(hh[ldh+i]);
1174 : #endif
1175 : #endif
1176 : #ifdef HAVE_SPARC64_SSE
1177 : #ifdef DOUBLE_PRECISION_REAL
1178 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
1179 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
1180 : #endif
1181 : #ifdef SINGLE_PRECISION_REAL
1182 : h1 = _mm_set_ps(hh[i-1], hh[i-1]);
1183 : h2 = _mm_set_ps(hh[ldh+i], hh[ldh+i]);
1184 : #endif
1185 : #endif
1186 :
1187 :
1188 0 : q1 = _SSE_LOAD(&q[i*ldq]);
1189 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
1190 0 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
1191 0 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
1192 0 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
1193 0 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
1194 0 : q3 = _SSE_LOAD(&q[(i*ldq)+2*offset]);
1195 0 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
1196 0 : y3 = _SSE_ADD(y3, _SSE_MUL(q3,h2));
1197 0 : q4 = _SSE_LOAD(&q[(i*ldq)+3*offset]);
1198 0 : x4 = _SSE_ADD(x4, _SSE_MUL(q4,h1));
1199 0 : y4 = _SSE_ADD(y4, _SSE_MUL(q4,h2));
1200 : }
1201 : #ifdef HAVE_SSE_INTRINSICS
1202 : #ifdef DOUBLE_PRECISION_REAL
1203 0 : h1 = _mm_set1_pd(hh[nb-1]);
1204 : #endif
1205 : #ifdef SINGLE_PRECISION_REAL
1206 0 : h1 = _mm_set1_ps(hh[nb-1]);
1207 : #endif
1208 : #endif
1209 : #ifdef HAVE_SPARC64_SSE
1210 : #ifdef DOUBLE_PRECISION_REAL
1211 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
1212 : #endif
1213 : #ifdef SINGLE_PRECISION_REAL
1214 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
1215 : #endif
1216 : #endif
1217 0 : q1 = _SSE_LOAD(&q[nb*ldq]);
1218 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
1219 0 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
1220 0 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
1221 0 : q3 = _SSE_LOAD(&q[(nb*ldq)+2*offset]);
1222 0 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
1223 0 : q4 = _SSE_LOAD(&q[(nb*ldq)+3*offset]);
1224 0 : x4 = _SSE_ADD(x4, _SSE_MUL(q4,h1));
1225 : /////////////////////////////////////////////////////
1226 : // Rank-2 update of Q [12 x nb+1]
1227 : /////////////////////////////////////////////////////
1228 : #ifdef HAVE_SSE_INTRINSICS
1229 : #ifdef DOUBLE_PRECISION_REAL
1230 0 : __SSE_DATATYPE tau1 = _mm_set1_pd(hh[0]);
1231 0 : __SSE_DATATYPE tau2 = _mm_set1_pd(hh[ldh]);
1232 0 : __SSE_DATATYPE vs = _mm_set1_pd(s);
1233 : #endif
1234 : #ifdef SINGLE_PRECISION_REAL
1235 0 : __SSE_DATATYPE tau1 = _mm_set1_ps(hh[0]);
1236 0 : __SSE_DATATYPE tau2 = _mm_set1_ps(hh[ldh]);
1237 0 : __SSE_DATATYPE vs = _mm_set1_ps(s);
1238 :
1239 : #endif
1240 : #endif
1241 : #ifdef HAVE_SPARC64_SSE
1242 : #ifdef DOUBLE_PRECISION_REAL
1243 : __SSE_DATATYPE tau1 = _mm_set_pd(hh[0], hh[0]);
1244 : __SSE_DATATYPE tau2 = _mm_set_pd(hh[ldh], hh[ldh]);
1245 : __SSE_DATATYPE vs = _mm_set_pd(s, s);
1246 : #endif
1247 : #ifdef SINGLE_PRECISION_REAL
1248 : __SSE_DATATYPE tau1 = _mm_set_ps(hh[0], hh[0]);
1249 : __SSE_DATATYPE tau2 = _mm_set_ps(hh[ldh], hh[ldh]);
1250 : __SSE_DATATYPE vs = _mm_set_ps(s, s);
1251 :
1252 : #endif
1253 : #endif
1254 :
1255 : #ifdef HAVE_SSE_INTRINSICS
1256 0 : h1 = _SSE_XOR(tau1, sign);
1257 : #endif
1258 : #ifdef HAVE_SPARC64_SSE
1259 : h1 = _fjsp_neg_v2r8(tau1);
1260 : #endif
1261 0 : x1 = _SSE_MUL(x1, h1);
1262 0 : x2 = _SSE_MUL(x2, h1);
1263 0 : x3 = _SSE_MUL(x3, h1);
1264 0 : x4 = _SSE_MUL(x4, h1);
1265 : #ifdef HAVE_SSE_INTRINSICS
1266 0 : h1 = _SSE_XOR(tau2, sign);
1267 : #endif
1268 : #ifdef HAVE_SPARC64_SSE
1269 : h1 = _fjsp_neg_v2r8(tau2);
1270 : #endif
1271 0 : h2 = _SSE_MUL(h1, vs);
1272 :
1273 0 : y1 = _SSE_ADD(_SSE_MUL(y1,h1), _SSE_MUL(x1,h2));
1274 0 : y2 = _SSE_ADD(_SSE_MUL(y2,h1), _SSE_MUL(x2,h2));
1275 0 : y3 = _SSE_ADD(_SSE_MUL(y3,h1), _SSE_MUL(x3,h2));
1276 0 : y4 = _SSE_ADD(_SSE_MUL(y4,h1), _SSE_MUL(x4,h2));
1277 0 : q1 = _SSE_LOAD(q);
1278 0 : q1 = _SSE_ADD(q1, y1);
1279 : _SSE_STORE(q,q1);
1280 0 : q2 = _SSE_LOAD(&q[offset]);
1281 0 : q2 = _SSE_ADD(q2, y2);
1282 0 : _SSE_STORE(&q[offset],q2);
1283 0 : q3 = _SSE_LOAD(&q[2*offset]);
1284 0 : q3 = _SSE_ADD(q3, y3);
1285 0 : _SSE_STORE(&q[2*offset],q3);
1286 0 : q4 = _SSE_LOAD(&q[3*offset]);
1287 0 : q4 = _SSE_ADD(q4, y4);
1288 0 : _SSE_STORE(&q[3*offset],q4);
1289 :
1290 : #ifdef HAVE_SSE_INTRINSICS
1291 : #ifdef DOUBLE_PRECISION_REAL
1292 0 : h2 = _mm_set1_pd(hh[ldh+1]);
1293 : #endif
1294 : #ifdef SINGLE_PRECISION_REAL
1295 0 : h2 = _mm_set1_ps(hh[ldh+1]);
1296 : #endif
1297 : #endif
1298 : #ifdef HAVE_SPARC64_SSE
1299 : #ifdef DOUBLE_PRECISION_REAL
1300 : h2 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1301 : #endif
1302 : #ifdef SINGLE_PRECISION_REAL
1303 : h2 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
1304 : #endif
1305 : #endif
1306 :
1307 0 : q1 = _SSE_LOAD(&q[ldq]);
1308 0 : q1 = _SSE_ADD(q1, _SSE_ADD(x1, _SSE_MUL(y1, h2)));
1309 0 : _SSE_STORE(&q[ldq],q1);
1310 0 : q2 = _SSE_LOAD(&q[ldq+offset]);
1311 0 : q2 = _SSE_ADD(q2, _SSE_ADD(x2, _SSE_MUL(y2, h2)));
1312 0 : _SSE_STORE(&q[ldq+offset],q2);
1313 0 : q3 = _SSE_LOAD(&q[ldq+2*offset]);
1314 0 : q3 = _SSE_ADD(q3, _SSE_ADD(x3, _SSE_MUL(y3, h2)));
1315 0 : _SSE_STORE(&q[ldq+2*offset],q3);
1316 0 : q4 = _SSE_LOAD(&q[ldq+3*offset]);
1317 0 : q4 = _SSE_ADD(q4, _SSE_ADD(x4, _SSE_MUL(y4, h2)));
1318 0 : _SSE_STORE(&q[ldq+3*offset],q4);
1319 :
1320 0 : for (i = 2; i < nb; i++)
1321 : {
1322 : #ifdef HAVE_SSE_INTRINSICS
1323 : #ifdef DOUBLE_PRECISION_REAL
1324 0 : h1 = _mm_set1_pd(hh[i-1]);
1325 0 : h2 = _mm_set1_pd(hh[ldh+i]);
1326 : #endif
1327 : #ifdef SINGLE_PRECISION_REAL
1328 0 : h1 = _mm_set1_ps(hh[i-1]);
1329 0 : h2 = _mm_set1_ps(hh[ldh+i]);
1330 : #endif
1331 : #endif
1332 : #ifdef HAVE_SPARC64_SSE
1333 : #ifdef DOUBLE_PRECISION_REAL
1334 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
1335 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
1336 : #endif
1337 : #ifdef SINGLE_PRECISION_REAL
1338 : h1 = _mm_set1_ps(hh[i-1], hh[i-1]);
1339 : h2 = _mm_set1_ps(hh[ldh+i], hh[ldh+i]);
1340 : #endif
1341 : #endif
1342 :
1343 :
1344 0 : q1 = _SSE_LOAD(&q[i*ldq]);
1345 0 : q1 = _SSE_ADD(q1, _SSE_ADD(_SSE_MUL(x1,h1), _SSE_MUL(y1, h2)));
1346 0 : _SSE_STORE(&q[i*ldq],q1);
1347 0 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
1348 0 : q2 = _SSE_ADD(q2, _SSE_ADD(_SSE_MUL(x2,h1), _SSE_MUL(y2, h2)));
1349 0 : _SSE_STORE(&q[(i*ldq)+offset],q2);
1350 0 : q3 = _SSE_LOAD(&q[(i*ldq)+2*offset]);
1351 0 : q3 = _SSE_ADD(q3, _SSE_ADD(_SSE_MUL(x3,h1), _SSE_MUL(y3, h2)));
1352 0 : _SSE_STORE(&q[(i*ldq)+2*offset],q3);
1353 0 : q4 = _SSE_LOAD(&q[(i*ldq)+3*offset]);
1354 0 : q4 = _SSE_ADD(q4, _SSE_ADD(_SSE_MUL(x4,h1), _SSE_MUL(y4, h2)));
1355 0 : _SSE_STORE(&q[(i*ldq)+3*offset],q4);
1356 : }
1357 : #ifdef HAVE_SSE_INTRINSICS
1358 : #ifdef DOUBLE_PRECISION_REAL
1359 0 : h1 = _mm_set1_pd(hh[nb-1]);
1360 : #endif
1361 : #ifdef SINGLE_PRECISION_REAL
1362 0 : h1 = _mm_set1_ps(hh[nb-1]);
1363 : #endif
1364 : #endif
1365 : #ifdef HAVE_SPARC64_SSE
1366 : #ifdef DOUBLE_PRECISION_REAL
1367 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
1368 : #endif
1369 : #ifdef SINGLE_PRECISION_REAL
1370 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
1371 : #endif
1372 : #endif
1373 :
1374 0 : q1 = _SSE_LOAD(&q[nb*ldq]);
1375 0 : q1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
1376 0 : _SSE_STORE(&q[nb*ldq],q1);
1377 0 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
1378 0 : q2 = _SSE_ADD(q2, _SSE_MUL(x2, h1));
1379 0 : _SSE_STORE(&q[(nb*ldq)+offset],q2);
1380 0 : q3 = _SSE_LOAD(&q[(nb*ldq)+2*offset]);
1381 0 : q3 = _SSE_ADD(q3, _SSE_MUL(x3, h1));
1382 0 : _SSE_STORE(&q[(nb*ldq)+2*offset],q3);
1383 0 : q4 = _SSE_LOAD(&q[(nb*ldq)+3*offset]);
1384 0 : q4 = _SSE_ADD(q4, _SSE_MUL(x4, h1));
1385 0 : _SSE_STORE(&q[(nb*ldq)+3*offset],q4);
1386 : }
1387 :
1388 : /**
1389 : * Unrolled kernel that computes
1390 : #ifdef DOUBLE_PRECISION_REAL
1391 : * 6 rows of Q simultaneously, a
1392 : #endif
1393 : #ifdef SINGLE_PRECISION_REAL
1394 : * 12 rows of Q simultaneously, a
1395 : #endif
1396 : * matrix Vector product with two householder
1397 : * vectors + a rank 2 update is performed
1398 : */
1399 : #ifdef HAVE_SSE_INTRINSICS
1400 : #ifdef DOUBLE_PRECISION_REAL
1401 : __forceinline void hh_trafo_kernel_6_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
1402 : #endif
1403 : #ifdef SINGLE_PRECISION_REAL
1404 : __forceinline void hh_trafo_kernel_12_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
1405 : #endif
1406 : #endif
1407 : #ifdef HAVE_SPARC64_SSE
1408 : #ifdef DOUBLE_PRECISION_REAL
1409 : __forceinline void hh_trafo_kernel_6_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
1410 : #endif
1411 : #ifdef SINGLE_PRECISION_REAL
1412 : __forceinline void hh_trafo_kernel_12_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
1413 : #endif
1414 :
1415 : #endif
1416 : {
1417 : /////////////////////////////////////////////////////
1418 : // Matrix Vector Multiplication, Q [12 x nb+1] * hh
1419 : // hh contains two householder vectors, with offset 1
1420 : /////////////////////////////////////////////////////
1421 : int i;
1422 : #ifdef HAVE_SSE_INTRINSICS
1423 : // Needed bit mask for floating point sign flip
1424 : #ifdef DOUBLE_PRECISION_REAL
1425 0 : __SSE_DATATYPE sign = (__SSE_DATATYPE)_mm_set1_epi64x(0x8000000000000000LL);
1426 : #endif
1427 : #ifdef SINGLE_PRECISION_REAL
1428 0 : __SSE_DATATYPE sign = _mm_castsi128_ps(_mm_set_epi32(0x80000000, 0x80000000, 0x80000000, 0x80000000));
1429 : #endif
1430 : #endif
1431 :
1432 0 : __SSE_DATATYPE x1 = _SSE_LOAD(&q[ldq]);
1433 0 : __SSE_DATATYPE x2 = _SSE_LOAD(&q[ldq+offset]);
1434 0 : __SSE_DATATYPE x3 = _SSE_LOAD(&q[ldq+2*offset]);
1435 :
1436 : #ifdef HAVE_SSE_INTRINSICS
1437 : #ifdef DOUBLE_PRECISION_REAL
1438 0 : __SSE_DATATYPE h1 = _mm_set1_pd(hh[ldh+1]);
1439 : #endif
1440 : #ifdef SINGLE_PRECISION_REAL
1441 0 : __SSE_DATATYPE h1 = _mm_set1_ps(hh[ldh+1]);
1442 : #endif
1443 : #endif
1444 : #ifdef HAVE_SPARC64_SSE
1445 : #ifdef DOUBLE_PRECISION_REAL
1446 : __SSE_DATATYPE h1 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1447 : #endif
1448 : #ifdef SINGLE_PRECISION_REAL
1449 : __SSE_DATATYPE h1 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
1450 : #endif
1451 : #endif
1452 : __SSE_DATATYPE h2;
1453 :
1454 0 : __SSE_DATATYPE q1 = _SSE_LOAD(q);
1455 0 : __SSE_DATATYPE y1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
1456 0 : __SSE_DATATYPE q2 = _SSE_LOAD(&q[offset]);
1457 0 : __SSE_DATATYPE y2 = _SSE_ADD(q2, _SSE_MUL(x2, h1));
1458 0 : __SSE_DATATYPE q3 = _SSE_LOAD(&q[2*offset]);
1459 0 : __SSE_DATATYPE y3 = _SSE_ADD(q3, _SSE_MUL(x3, h1));
1460 0 : for(i = 2; i < nb; i++)
1461 : {
1462 : #ifdef HAVE_SSE_INTRINSICS
1463 : #ifdef DOUBLE_PRECISION_REAL
1464 0 : h1 = _mm_set1_pd(hh[i-1]);
1465 0 : h2 = _mm_set1_pd(hh[ldh+i]);
1466 : #endif
1467 : #ifdef SINGLE_PRECISION_REAL
1468 0 : h1 = _mm_set1_ps(hh[i-1]);
1469 0 : h2 = _mm_set1_ps(hh[ldh+i]);
1470 : #endif
1471 : #endif
1472 : #ifdef HAVE_SPARC64_SSE
1473 : #ifdef DOUBLE_PRECISION_REAL
1474 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
1475 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
1476 : #endif
1477 : #ifdef SINGLE_PRECISION_REAL
1478 : h1 = _mm_set1_ps(hh[i-1], hh[i-1]);
1479 : h2 = _mm_set1_ps(hh[ldh+i], hh[ldh+i]);
1480 : #endif
1481 : #endif
1482 :
1483 :
1484 0 : q1 = _SSE_LOAD(&q[i*ldq]);
1485 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
1486 0 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
1487 0 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
1488 0 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
1489 0 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
1490 0 : q3 = _SSE_LOAD(&q[(i*ldq)+2*offset]);
1491 0 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
1492 0 : y3 = _SSE_ADD(y3, _SSE_MUL(q3,h2));
1493 : }
1494 : #ifdef HAVE_SSE_INTRINSICS
1495 : #ifdef DOUBLE_PRECISION_REAL
1496 0 : h1 = _mm_set1_pd(hh[nb-1]);
1497 : #endif
1498 : #ifdef SINGLE_PRECISION_REAL
1499 0 : h1 = _mm_set1_ps(hh[nb-1]);
1500 : #endif
1501 : #endif
1502 : #ifdef HAVE_SPARC64_SSE
1503 : #ifdef DOUBLE_PRECISION_REAL
1504 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
1505 : #endif
1506 : #ifdef SINGLE_PRECISION_REAL
1507 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
1508 : #endif
1509 : #endif
1510 :
1511 0 : q1 = _SSE_LOAD(&q[nb*ldq]);
1512 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
1513 0 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
1514 0 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
1515 0 : q3 = _SSE_LOAD(&q[(nb*ldq)+2*offset]);
1516 0 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
1517 : /////////////////////////////////////////////////////
1518 : // Rank-2 update of Q [12 x nb+1]
1519 : /////////////////////////////////////////////////////
1520 : #ifdef HAVE_SSE_INTRINSICS
1521 : #ifdef DOUBLE_PRECISION_REAL
1522 0 : __SSE_DATATYPE tau1 = _mm_set1_pd(hh[0]);
1523 0 : __SSE_DATATYPE tau2 = _mm_set1_pd(hh[ldh]);
1524 0 : __SSE_DATATYPE vs = _mm_set1_pd(s);
1525 : #endif
1526 : #ifdef SINGLE_PRECISION_REAL
1527 0 : __SSE_DATATYPE tau1 = _mm_set1_ps(hh[0]);
1528 0 : __SSE_DATATYPE tau2 = _mm_set1_ps(hh[ldh]);
1529 0 : __SSE_DATATYPE vs = _mm_set1_ps(s);
1530 :
1531 : #endif
1532 : #endif
1533 : #ifdef HAVE_SPARC64_SSE
1534 : #ifdef DOUBLE_PRECISION_REAL
1535 : __SSE_DATATYPE tau1 = _mm_set_pd(hh[0], hh[0]);
1536 : __SSE_DATATYPE tau2 = _mm_set_pd(hh[ldh], hh[ldh]);
1537 : __SSE_DATATYPE vs = _mm_set_pd(s, s);
1538 : #endif
1539 : #ifdef SINGLE_PRECISION_REAL
1540 : __SSE_DATATYPE tau1 = _mm_set1_ps(hh[0], hh[0]);
1541 : __SSE_DATATYPE tau2 = _mm_set1_ps(hh[ldh], hh[ldh]);
1542 : __SSE_DATATYPE vs = _mm_set1_ps(s, s);
1543 :
1544 : #endif
1545 : #endif
1546 :
1547 :
1548 : #ifdef HAVE_SSE_INTRINSICS
1549 0 : h1 = _SSE_XOR(tau1, sign);
1550 : #endif
1551 : #ifdef HAVE_SPARC64_SSE
1552 : h1 = _fjsp_neg_v2r8(tau1);
1553 : #endif
1554 0 : x1 = _SSE_MUL(x1, h1);
1555 0 : x2 = _SSE_MUL(x2, h1);
1556 0 : x3 = _SSE_MUL(x3, h1);
1557 : #ifdef HAVE_SSE_INTRINSICS
1558 0 : h1 = _SSE_XOR(tau2, sign);
1559 : #endif
1560 : #ifdef HAVE_SPARC64_SSE
1561 : h1 = _fjsp_neg_v2r8(tau2);
1562 : #endif
1563 0 : h2 = _SSE_MUL(h1, vs);
1564 :
1565 0 : y1 = _SSE_ADD(_SSE_MUL(y1,h1), _SSE_MUL(x1,h2));
1566 0 : y2 = _SSE_ADD(_SSE_MUL(y2,h1), _SSE_MUL(x2,h2));
1567 0 : y3 = _SSE_ADD(_SSE_MUL(y3,h1), _SSE_MUL(x3,h2));
1568 0 : q1 = _SSE_LOAD(q);
1569 0 : q1 = _SSE_ADD(q1, y1);
1570 : _SSE_STORE(q,q1);
1571 0 : q2 = _SSE_LOAD(&q[offset]);
1572 0 : q2 = _SSE_ADD(q2, y2);
1573 0 : _SSE_STORE(&q[offset],q2);
1574 0 : q3 = _SSE_LOAD(&q[2*offset]);
1575 0 : q3 = _SSE_ADD(q3, y3);
1576 0 : _SSE_STORE(&q[2*offset],q3);
1577 :
1578 : #ifdef HAVE_SSE_INTRINSICS
1579 : #ifdef DOUBLE_PRECISION_REAL
1580 0 : h2 = _mm_set1_pd(hh[ldh+1]);
1581 : #endif
1582 : #ifdef SINGLE_PRECISION_REAL
1583 0 : h2 = _mm_set1_ps(hh[ldh+1]);
1584 : #endif
1585 : #endif
1586 : #ifdef HAVE_SPARC64_SSE
1587 : #ifdef DOUBLE_PRECISION_REAL
1588 : h2 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1589 : #endif
1590 : #ifdef SINGLE_PRECISION_REAL
1591 : h2 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
1592 : #endif
1593 : #endif
1594 :
1595 :
1596 :
1597 0 : q1 = _SSE_LOAD(&q[ldq]);
1598 0 : q1 = _SSE_ADD(q1, _SSE_ADD(x1, _SSE_MUL(y1, h2)));
1599 0 : _SSE_STORE(&q[ldq],q1);
1600 0 : q2 = _SSE_LOAD(&q[ldq+offset]);
1601 0 : q2 = _SSE_ADD(q2, _SSE_ADD(x2, _SSE_MUL(y2, h2)));
1602 0 : _SSE_STORE(&q[ldq+offset],q2);
1603 0 : q3 = _SSE_LOAD(&q[ldq+2*offset]);
1604 0 : q3 = _SSE_ADD(q3, _SSE_ADD(x3, _SSE_MUL(y3, h2)));
1605 0 : _SSE_STORE(&q[ldq+2*offset],q3);
1606 :
1607 0 : for (i = 2; i < nb; i++)
1608 : {
1609 : #ifdef HAVE_SSE_INTRINSICS
1610 : #ifdef DOUBLE_PRECISION_REAL
1611 0 : h1 = _mm_set1_pd(hh[i-1]);
1612 0 : h2 = _mm_set1_pd(hh[ldh+i]);
1613 : #endif
1614 : #ifdef SINGLE_PRECISION_REAL
1615 0 : h1 = _mm_set1_ps(hh[i-1]);
1616 0 : h2 = _mm_set1_ps(hh[ldh+i]);
1617 : #endif
1618 : #endif
1619 : #ifdef HAVE_SPARC64_SSE
1620 : #ifdef DOUBLE_PRECISION_REAL
1621 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
1622 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
1623 : #endif
1624 : #ifdef SINGLE_PRECISION_REAL
1625 : h1 = _mm_set1_ps(hh[i-1], hh[i-1]);
1626 : h2 = _mm_set1_ps(hh[ldh+i], hh[ldh+i]);
1627 : #endif
1628 : #endif
1629 :
1630 0 : q1 = _SSE_LOAD(&q[i*ldq]);
1631 0 : q1 = _SSE_ADD(q1, _SSE_ADD(_SSE_MUL(x1,h1), _SSE_MUL(y1, h2)));
1632 0 : _SSE_STORE(&q[i*ldq],q1);
1633 0 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
1634 0 : q2 = _SSE_ADD(q2, _SSE_ADD(_SSE_MUL(x2,h1), _SSE_MUL(y2, h2)));
1635 0 : _SSE_STORE(&q[(i*ldq)+offset],q2);
1636 0 : q3 = _SSE_LOAD(&q[(i*ldq)+2*offset]);
1637 0 : q3 = _SSE_ADD(q3, _SSE_ADD(_SSE_MUL(x3,h1), _SSE_MUL(y3, h2)));
1638 0 : _SSE_STORE(&q[(i*ldq)+2*offset],q3);
1639 : }
1640 : #ifdef HAVE_SSE_INTRINSICS
1641 : #ifdef DOUBLE_PRECISION_REAL
1642 0 : h1 = _mm_set1_pd(hh[nb-1]);
1643 : #endif
1644 : #ifdef SINGLE_PRECISION_REAL
1645 0 : h1 = _mm_set1_ps(hh[nb-1]);
1646 : #endif
1647 : #endif
1648 : #ifdef HAVE_SPARC64_SSE
1649 : #ifdef DOUBLE_PRECISION_REAL
1650 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
1651 : #endif
1652 : #ifdef SINGLE_PRECISION_REAL
1653 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
1654 : #endif
1655 : #endif
1656 0 : q1 = _SSE_LOAD(&q[nb*ldq]);
1657 0 : q1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
1658 0 : _SSE_STORE(&q[nb*ldq],q1);
1659 0 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
1660 0 : q2 = _SSE_ADD(q2, _SSE_MUL(x2, h1));
1661 0 : _SSE_STORE(&q[(nb*ldq)+offset],q2);
1662 0 : q3 = _SSE_LOAD(&q[(nb*ldq)+2*offset]);
1663 0 : q3 = _SSE_ADD(q3, _SSE_MUL(x3, h1));
1664 0 : _SSE_STORE(&q[(nb*ldq)+2*offset],q3);
1665 : }
1666 :
1667 :
1668 : /**
1669 : * Unrolled kernel that computes
1670 : #ifdef DOUBLE_PRECISION_REAL
1671 : * 4 rows of Q simultaneously, a
1672 : #endif
1673 : #ifdef SINGLE_PRECISION_REAL
1674 : * 8 rows of Q simultaneously, a
1675 : #endif
1676 : * matrix Vector product with two householder
1677 : * vectors + a rank 2 update is performed
1678 : */
1679 : #ifdef HAVE_SSE_INTRINSICS
1680 : #ifdef DOUBLE_PRECISION_REAL
1681 : __forceinline void hh_trafo_kernel_4_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
1682 : #endif
1683 : #ifdef SINGLE_PRECISION_REAL
1684 : __forceinline void hh_trafo_kernel_8_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
1685 : #endif
1686 : #endif
1687 : #ifdef HAVE_SPARC64_SSE
1688 : #ifdef DOUBLE_PRECISION_REAL
1689 : __forceinline void hh_trafo_kernel_4_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
1690 : #endif
1691 : #ifdef SINGLE_PRECISION_REAL
1692 : __forceinline void hh_trafo_kernel_8_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
1693 : #endif
1694 :
1695 : #endif
1696 : {
1697 : /////////////////////////////////////////////////////
1698 : // Matrix Vector Multiplication, Q [12 x nb+1] * hh
1699 : // hh contains two householder vectors, with offset 1
1700 : /////////////////////////////////////////////////////
1701 : int i;
1702 : #ifdef HAVE_SSE_INTRINSICS
1703 : // Needed bit mask for floating point sign flip
1704 : #ifdef DOUBLE_PRECISION_REAL
1705 544768 : __SSE_DATATYPE sign = (__SSE_DATATYPE)_mm_set1_epi64x(0x8000000000000000LL);
1706 : #endif
1707 : #ifdef SINGLE_PRECISION_REAL
1708 204288 : __SSE_DATATYPE sign = _mm_castsi128_ps(_mm_set_epi32(0x80000000, 0x80000000, 0x80000000, 0x80000000));
1709 : #endif
1710 : #endif
1711 :
1712 1293824 : __SSE_DATATYPE x1 = _SSE_LOAD(&q[ldq]);
1713 1293824 : __SSE_DATATYPE x2 = _SSE_LOAD(&q[ldq+offset]);
1714 :
1715 : #ifdef HAVE_SSE_INTRINSICS
1716 : #ifdef DOUBLE_PRECISION_REAL
1717 1089536 : __SSE_DATATYPE h1 = _mm_set1_pd(hh[ldh+1]);
1718 : #endif
1719 : #ifdef SINGLE_PRECISION_REAL
1720 204288 : __SSE_DATATYPE h1 = _mm_set1_ps(hh[ldh+1]);
1721 : #endif
1722 : #endif
1723 : #ifdef HAVE_SPARC64_SSE
1724 : #ifdef DOUBLE_PRECISION_REAL
1725 : __SSE_DATATYPE h1 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1726 : #endif
1727 : #ifdef SINGLE_PRECISION_REAL
1728 : __SSE_DATATYPE h1 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
1729 : #endif
1730 : #endif
1731 : __SSE_DATATYPE h2;
1732 :
1733 646912 : __SSE_DATATYPE q1 = _SSE_LOAD(q);
1734 1293824 : __SSE_DATATYPE y1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
1735 1293824 : __SSE_DATATYPE q2 = _SSE_LOAD(&q[offset]);
1736 1293824 : __SSE_DATATYPE y2 = _SSE_ADD(q2, _SSE_MUL(x2, h1));
1737 40755456 : for(i = 2; i < nb; i++)
1738 : {
1739 : #ifdef HAVE_SSE_INTRINSICS
1740 : #ifdef DOUBLE_PRECISION_REAL
1741 67551232 : h1 = _mm_set1_pd(hh[i-1]);
1742 67551232 : h2 = _mm_set1_pd(hh[ldh+i]);
1743 : #endif
1744 : #ifdef SINGLE_PRECISION_REAL
1745 12665856 : h1 = _mm_set1_ps(hh[i-1]);
1746 12665856 : h2 = _mm_set1_ps(hh[ldh+i]);
1747 : #endif
1748 : #endif
1749 : #ifdef HAVE_SPARC64_SSE
1750 : #ifdef DOUBLE_PRECISION_REAL
1751 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
1752 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
1753 : #endif
1754 : #ifdef SINGLE_PRECISION_REAL
1755 : h1 = _mm_set_ps(hh[i-1], hh[i-1]);
1756 : h2 = _mm_set_ps(hh[ldh+i], hh[ldh+i]);
1757 : #endif
1758 : #endif
1759 :
1760 80217088 : q1 = _SSE_LOAD(&q[i*ldq]);
1761 80217088 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
1762 80217088 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
1763 80217088 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
1764 80217088 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
1765 80217088 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
1766 : }
1767 :
1768 : #ifdef HAVE_SSE_INTRINSICS
1769 : #ifdef DOUBLE_PRECISION_REAL
1770 1089536 : h1 = _mm_set1_pd(hh[nb-1]);
1771 : #endif
1772 : #ifdef SINGLE_PRECISION_REAL
1773 204288 : h1 = _mm_set1_ps(hh[nb-1]);
1774 : #endif
1775 : #endif
1776 : #ifdef HAVE_SPARC64_SSE
1777 : #ifdef DOUBLE_PRECISION_REAL
1778 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
1779 : #endif
1780 : #ifdef SINGLE_PRECISION_REAL
1781 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
1782 : #endif
1783 : #endif
1784 :
1785 1293824 : q1 = _SSE_LOAD(&q[nb*ldq]);
1786 1293824 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
1787 1293824 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
1788 1293824 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
1789 : /////////////////////////////////////////////////////
1790 : // Rank-2 update of Q [12 x nb+1]
1791 : /////////////////////////////////////////////////////
1792 : #ifdef HAVE_SSE_INTRINSICS
1793 : #ifdef DOUBLE_PRECISION_REAL
1794 1089536 : __SSE_DATATYPE tau1 = _mm_set1_pd(hh[0]);
1795 1089536 : __SSE_DATATYPE tau2 = _mm_set1_pd(hh[ldh]);
1796 544768 : __SSE_DATATYPE vs = _mm_set1_pd(s);
1797 : #endif
1798 : #ifdef SINGLE_PRECISION_REAL
1799 204288 : __SSE_DATATYPE tau1 = _mm_set1_ps(hh[0]);
1800 204288 : __SSE_DATATYPE tau2 = _mm_set1_ps(hh[ldh]);
1801 102144 : __SSE_DATATYPE vs = _mm_set1_ps(s);
1802 :
1803 : #endif
1804 : #endif
1805 : #ifdef HAVE_SPARC64_SSE
1806 : #ifdef DOUBLE_PRECISION_REAL
1807 : __SSE_DATATYPE tau1 = _mm_set_pd(hh[0], hh[0]);
1808 : __SSE_DATATYPE tau2 = _mm_set_pd(hh[ldh], hh[ldh]);
1809 : __SSE_DATATYPE vs = _mm_set_pd(s, s);
1810 : #endif
1811 : #ifdef SINGLE_PRECISION_REAL
1812 : __SSE_DATATYPE tau1 = _mm_set_ps(hh[0], hh[0]);
1813 : __SSE_DATATYPE tau2 = _mm_set_ps(hh[ldh], hh[ldh]);
1814 : __SSE_DATATYPE vs = _mm_set_ps(s, s);
1815 :
1816 : #endif
1817 : #endif
1818 :
1819 : #ifdef HAVE_SSE_INTRINSICS
1820 646912 : h1 = _SSE_XOR(tau1, sign);
1821 : #endif
1822 : #ifdef HAVE_SPARC64_SSE
1823 : h1 = _fjsp_neg_v2r8(tau1);
1824 : #endif
1825 646912 : x1 = _SSE_MUL(x1, h1);
1826 646912 : x2 = _SSE_MUL(x2, h1);
1827 : #ifdef HAVE_SSE_INTRINSICS
1828 646912 : h1 = _SSE_XOR(tau2, sign);
1829 : #endif
1830 : #ifdef HAVE_SPARC64_SSE
1831 : h1 = _fjsp_neg_v2r8(tau2);
1832 : #endif
1833 646912 : h2 = _SSE_MUL(h1, vs);
1834 :
1835 1940736 : y1 = _SSE_ADD(_SSE_MUL(y1,h1), _SSE_MUL(x1,h2));
1836 1940736 : y2 = _SSE_ADD(_SSE_MUL(y2,h1), _SSE_MUL(x2,h2));
1837 646912 : q1 = _SSE_LOAD(q);
1838 646912 : q1 = _SSE_ADD(q1, y1);
1839 : _SSE_STORE(q,q1);
1840 1293824 : q2 = _SSE_LOAD(&q[offset]);
1841 646912 : q2 = _SSE_ADD(q2, y2);
1842 646912 : _SSE_STORE(&q[offset],q2);
1843 :
1844 : #ifdef HAVE_SSE_INTRINSICS
1845 : #ifdef DOUBLE_PRECISION_REAL
1846 1089536 : h2 = _mm_set1_pd(hh[ldh+1]);
1847 : #endif
1848 : #ifdef SINGLE_PRECISION_REAL
1849 204288 : h2 = _mm_set1_ps(hh[ldh+1]);
1850 : #endif
1851 : #endif
1852 : #ifdef HAVE_SPARC64_SSE
1853 : #ifdef DOUBLE_PRECISION_REAL
1854 : h2 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1855 : #endif
1856 : #ifdef SINGLE_PRECISION_REAL
1857 : h2 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
1858 : #endif
1859 : #endif
1860 :
1861 1293824 : q1 = _SSE_LOAD(&q[ldq]);
1862 1940736 : q1 = _SSE_ADD(q1, _SSE_ADD(x1, _SSE_MUL(y1, h2)));
1863 646912 : _SSE_STORE(&q[ldq],q1);
1864 1293824 : q2 = _SSE_LOAD(&q[ldq+offset]);
1865 1940736 : q2 = _SSE_ADD(q2, _SSE_ADD(x2, _SSE_MUL(y2, h2)));
1866 646912 : _SSE_STORE(&q[ldq+offset],q2);
1867 :
1868 40755456 : for (i = 2; i < nb; i++)
1869 : {
1870 : #ifdef HAVE_SSE_INTRINSICS
1871 : #ifdef DOUBLE_PRECISION_REAL
1872 67551232 : h1 = _mm_set1_pd(hh[i-1]);
1873 67551232 : h2 = _mm_set1_pd(hh[ldh+i]);
1874 : #endif
1875 : #ifdef SINGLE_PRECISION_REAL
1876 12665856 : h1 = _mm_set1_ps(hh[i-1]);
1877 12665856 : h2 = _mm_set1_ps(hh[ldh+i]);
1878 : #endif
1879 : #endif
1880 : #ifdef HAVE_SPARC64_SSE
1881 : #ifdef DOUBLE_PRECISION_REAL
1882 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
1883 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
1884 : #endif
1885 : #ifdef SINGLE_PRECISION_REAL
1886 : h1 = _mm_set_ps(hh[i-1], hh[i-1]);
1887 : h2 = _mm_set_ps(hh[ldh+i], hh[ldh+i]);
1888 : #endif
1889 : #endif
1890 :
1891 80217088 : q1 = _SSE_LOAD(&q[i*ldq]);
1892 160434176 : q1 = _SSE_ADD(q1, _SSE_ADD(_SSE_MUL(x1,h1), _SSE_MUL(y1, h2)));
1893 40108544 : _SSE_STORE(&q[i*ldq],q1);
1894 80217088 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
1895 160434176 : q2 = _SSE_ADD(q2, _SSE_ADD(_SSE_MUL(x2,h1), _SSE_MUL(y2, h2)));
1896 40108544 : _SSE_STORE(&q[(i*ldq)+offset],q2);
1897 : }
1898 : #ifdef HAVE_SSE_INTRINSICS
1899 : #ifdef DOUBLE_PRECISION_REAL
1900 1089536 : h1 = _mm_set1_pd(hh[nb-1]);
1901 : #endif
1902 : #ifdef SINGLE_PRECISION_REAL
1903 204288 : h1 = _mm_set1_ps(hh[nb-1]);
1904 : #endif
1905 : #endif
1906 : #ifdef HAVE_SPARC64_SSE
1907 : #ifdef DOUBLE_PRECISION_REAL
1908 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
1909 : #endif
1910 : #ifdef SINGLE_PRECISION_REAL
1911 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
1912 : #endif
1913 : #endif
1914 1293824 : q1 = _SSE_LOAD(&q[nb*ldq]);
1915 1293824 : q1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
1916 646912 : _SSE_STORE(&q[nb*ldq],q1);
1917 1293824 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
1918 1293824 : q2 = _SSE_ADD(q2, _SSE_MUL(x2, h1));
1919 646912 : _SSE_STORE(&q[(nb*ldq)+offset],q2);
1920 : }
1921 :
1922 :
1923 : /**
1924 : * Unrolled kernel that computes
1925 : #ifdef DOUBLE_PRECISION_REAL
1926 : * 2 rows of Q simultaneously, a
1927 : #endif
1928 : #ifdef SINGLE_PRECISION_REAL
1929 : * 4 rows of Q simultaneously, a
1930 : #endif
1931 : * matrix Vector product with two householder
1932 : * vectors + a rank 2 update is performed
1933 : */
1934 : #ifdef HAVE_SSE_INTRINSICS
1935 : #ifdef DOUBLE_PRECISION_REAL
1936 : __forceinline void hh_trafo_kernel_2_SSE_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
1937 : #endif
1938 : #ifdef SINGLE_PRECISION_REAL
1939 : __forceinline void hh_trafo_kernel_4_SSE_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
1940 : #endif
1941 : #endif
1942 : #ifdef HAVE_SPARC64_SSE
1943 : #ifdef DOUBLE_PRECISION_REAL
1944 : __forceinline void hh_trafo_kernel_2_SPARC64_2hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s)
1945 : #endif
1946 : #ifdef SINGLE_PRECISION_REAL
1947 : __forceinline void hh_trafo_kernel_4_SPARC64_2hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s)
1948 : #endif
1949 :
1950 : #endif
1951 : {
1952 : /////////////////////////////////////////////////////
1953 : // Matrix Vector Multiplication, Q [12 x nb+1] * hh
1954 : // hh contains two householder vectors, with offset 1
1955 : /////////////////////////////////////////////////////
1956 : int i;
1957 : #ifdef HAVE_SSE_INTRINSICS
1958 : // Needed bit mask for floating point sign flip
1959 : #ifdef DOUBLE_PRECISION_REAL
1960 0 : __SSE_DATATYPE sign = (__SSE_DATATYPE)_mm_set1_epi64x(0x8000000000000000LL);
1961 : #endif
1962 : #ifdef SINGLE_PRECISION_REAL
1963 0 : __SSE_DATATYPE sign = _mm_castsi128_ps(_mm_set_epi32(0x80000000, 0x80000000, 0x80000000, 0x80000000));
1964 : #endif
1965 : #endif
1966 :
1967 0 : __SSE_DATATYPE x1 = _SSE_LOAD(&q[ldq]);
1968 :
1969 : #ifdef HAVE_SSE_INTRINSICS
1970 : #ifdef DOUBLE_PRECISION_REAL
1971 0 : __SSE_DATATYPE h1 = _mm_set1_pd(hh[ldh+1]);
1972 : #endif
1973 : #ifdef SINGLE_PRECISION_REAL
1974 0 : __SSE_DATATYPE h1 = _mm_set1_ps(hh[ldh+1]);
1975 : #endif
1976 : #endif
1977 : #ifdef HAVE_SPARC64_SSE
1978 : #ifdef DOUBLE_PRECISION_REAL
1979 : __SSE_DATATYPE h1 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1980 : #endif
1981 : #ifdef SINGLE_PRECISION_REAL
1982 : __SSE_DATATYPE h1 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
1983 : #endif
1984 : #endif
1985 : __SSE_DATATYPE h2;
1986 :
1987 0 : __SSE_DATATYPE q1 = _SSE_LOAD(q);
1988 0 : __SSE_DATATYPE y1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
1989 0 : for(i = 2; i < nb; i++)
1990 : {
1991 : #ifdef HAVE_SSE_INTRINSICS
1992 : #ifdef DOUBLE_PRECISION_REAL
1993 0 : h1 = _mm_set1_pd(hh[i-1]);
1994 0 : h2 = _mm_set1_pd(hh[ldh+i]);
1995 : #endif
1996 : #ifdef SINGLE_PRECISION_REAL
1997 0 : h1 = _mm_set1_ps(hh[i-1]);
1998 0 : h2 = _mm_set1_ps(hh[ldh+i]);
1999 : #endif
2000 : #endif
2001 : #ifdef HAVE_SPARC64_SSE
2002 : #ifdef DOUBLE_PRECISION_REAL
2003 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
2004 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
2005 : #endif
2006 : #ifdef SINGLE_PRECISION_REAL
2007 : h1 = _mm_set_ps(hh[i-1], hh[i-1]);
2008 : h2 = _mm_set_ps(hh[ldh+i], hh[ldh+i]);
2009 : #endif
2010 : #endif
2011 :
2012 0 : q1 = _SSE_LOAD(&q[i*ldq]);
2013 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
2014 0 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
2015 : }
2016 :
2017 : #ifdef HAVE_SSE_INTRINSICS
2018 : #ifdef DOUBLE_PRECISION_REAL
2019 0 : h1 = _mm_set1_pd(hh[nb-1]);
2020 : #endif
2021 : #ifdef SINGLE_PRECISION_REAL
2022 0 : h1 = _mm_set1_ps(hh[nb-1]);
2023 : #endif
2024 : #endif
2025 : #ifdef HAVE_SPARC64_SSE
2026 : #ifdef DOUBLE_PRECISION_REAL
2027 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
2028 : #endif
2029 : #ifdef SINGLE_PRECISION_REAL
2030 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
2031 : #endif
2032 : #endif
2033 :
2034 :
2035 0 : q1 = _SSE_LOAD(&q[nb*ldq]);
2036 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
2037 : /////////////////////////////////////////////////////
2038 : // Rank-2 update of Q [12 x nb+1]
2039 : /////////////////////////////////////////////////////
2040 : #ifdef HAVE_SSE_INTRINSICS
2041 : #ifdef DOUBLE_PRECISION_REAL
2042 0 : __SSE_DATATYPE tau1 = _mm_set1_pd(hh[0]);
2043 0 : __SSE_DATATYPE tau2 = _mm_set1_pd(hh[ldh]);
2044 0 : __SSE_DATATYPE vs = _mm_set1_pd(s);
2045 : #endif
2046 : #ifdef SINGLE_PRECISION_REAL
2047 0 : __SSE_DATATYPE tau1 = _mm_set1_ps(hh[0]);
2048 0 : __SSE_DATATYPE tau2 = _mm_set1_ps(hh[ldh]);
2049 0 : __SSE_DATATYPE vs = _mm_set1_ps(s);
2050 :
2051 : #endif
2052 : #endif
2053 : #ifdef HAVE_SPARC64_SSE
2054 : #ifdef DOUBLE_PRECISION_REAL
2055 : __SSE_DATATYPE tau1 = _mm_set_pd(hh[0], hh[0]);
2056 : __SSE_DATATYPE tau2 = _mm_set_pd(hh[ldh], hh[ldh]);
2057 : __SSE_DATATYPE vs = _mm_set_pd(s,s );
2058 : #endif
2059 : #ifdef SINGLE_PRECISION_REAL
2060 : __SSE_DATATYPE tau1 = _mm_set_ps(hh[0], hh[0]);
2061 : __SSE_DATATYPE tau2 = _mm_set_ps(hh[ldh], hh[ldh]);
2062 : __SSE_DATATYPE vs = _mm_set_ps(s, s);
2063 :
2064 : #endif
2065 : #endif
2066 :
2067 : #ifdef HAVE_SSE_INTRINSICS
2068 0 : h1 = _SSE_XOR(tau1, sign);
2069 : #endif
2070 : #ifdef HAVE_SPARC64_SSE
2071 : h1 = _fjsp_neg_v2r8(tau1);
2072 : #endif
2073 0 : x1 = _SSE_MUL(x1, h1);
2074 : #ifdef HAVE_SSE_INTRINSICS
2075 0 : h1 = _SSE_XOR(tau2, sign);
2076 : #endif
2077 : #ifdef HAVE_SPARC64_SSE
2078 : h1 = _fjsp_neg_v2r8(tau2);
2079 : #endif
2080 0 : h2 = _SSE_MUL(h1, vs);
2081 :
2082 0 : y1 = _SSE_ADD(_SSE_MUL(y1,h1), _SSE_MUL(x1,h2));
2083 0 : q1 = _SSE_LOAD(q);
2084 0 : q1 = _SSE_ADD(q1, y1);
2085 : _SSE_STORE(q,q1);
2086 :
2087 : #ifdef HAVE_SSE_INTRINSICS
2088 : #ifdef DOUBLE_PRECISION_REAL
2089 0 : h2 = _mm_set1_pd(hh[ldh+1]);
2090 : #endif
2091 : #ifdef SINGLE_PRECISION_REAL
2092 0 : h2 = _mm_set1_ps(hh[ldh+1]);
2093 : #endif
2094 : #endif
2095 : #ifdef HAVE_SPARC64_SSE
2096 : #ifdef DOUBLE_PRECISION_REAL
2097 : h2 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
2098 : #endif
2099 : #ifdef SINGLE_PRECISION_REAL
2100 : h2 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
2101 : #endif
2102 : #endif
2103 :
2104 0 : q1 = _SSE_LOAD(&q[ldq]);
2105 0 : q1 = _SSE_ADD(q1, _SSE_ADD(x1, _SSE_MUL(y1, h2)));
2106 0 : _SSE_STORE(&q[ldq],q1);
2107 :
2108 0 : for (i = 2; i < nb; i++)
2109 : {
2110 : #ifdef HAVE_SSE_INTRINSICS
2111 : #ifdef DOUBLE_PRECISION_REAL
2112 0 : h1 = _mm_set1_pd(hh[i-1]);
2113 0 : h2 = _mm_set1_pd(hh[ldh+i]);
2114 : #endif
2115 : #ifdef SINGLE_PRECISION_REAL
2116 0 : h1 = _mm_set1_ps(hh[i-1]);
2117 0 : h2 = _mm_set1_ps(hh[ldh+i]);
2118 : #endif
2119 : #endif
2120 : #ifdef HAVE_SPARC64_SSE
2121 : #ifdef DOUBLE_PRECISION_REAL
2122 : h1 = _mm_set_pd(hh[i-1], hh[i-1]);
2123 : h2 = _mm_set_pd(hh[ldh+i], hh[ldh+i]);
2124 : #endif
2125 : #ifdef SINGLE_PRECISION_REAL
2126 : h1 = _mm_set_ps(hh[i-1], hh[i-1]);
2127 : h2 = _mm_set_ps(hh[ldh+i]), hh[ldh+i];
2128 : #endif
2129 : #endif
2130 :
2131 0 : q1 = _SSE_LOAD(&q[i*ldq]);
2132 0 : q1 = _SSE_ADD(q1, _SSE_ADD(_SSE_MUL(x1,h1), _SSE_MUL(y1, h2)));
2133 0 : _SSE_STORE(&q[i*ldq],q1);
2134 : }
2135 : #ifdef HAVE_SSE_INTRINSICS
2136 : #ifdef DOUBLE_PRECISION_REAL
2137 0 : h1 = _mm_set1_pd(hh[nb-1]);
2138 : #endif
2139 : #ifdef SINGLE_PRECISION_REAL
2140 0 : h1 = _mm_set1_ps(hh[nb-1]);
2141 : #endif
2142 : #endif
2143 : #ifdef HAVE_SPARC64_SSE
2144 : #ifdef DOUBLE_PRECISION_REAL
2145 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
2146 : #endif
2147 : #ifdef SINGLE_PRECISION_REAL
2148 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
2149 : #endif
2150 : #endif
2151 0 : q1 = _SSE_LOAD(&q[nb*ldq]);
2152 0 : q1 = _SSE_ADD(q1, _SSE_MUL(x1, h1));
2153 0 : _SSE_STORE(&q[nb*ldq],q1);
2154 : }
|
