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 : //
46 : // --------------------------------------------------------------------------------------------------
47 : //
48 : // This file contains the compute intensive kernels for the Householder transformations.
49 : // It should be compiled with the highest possible optimization level.
50 : //
51 : // On Intel Nehalem or Intel Westmere or AMD Magny Cours use -O3 -msse3
52 : // On Intel Sandy Bridge use -O3 -mavx
53 : //
54 : // Copyright of the original code rests with the authors inside the ELPA
55 : // consortium. The copyright of any additional modifications shall rest
56 : // with their original authors, but shall adhere to the licensing terms
57 : // distributed along with the original code in the file "COPYING".
58 : //
59 : // Author: Alexander Heinecke (alexander.heinecke@mytum.de)
60 : // Adapted for building a shared-library by Andreas Marek, MPCDF (andreas.marek@mpcdf.mpg.de)
61 : // --------------------------------------------------------------------------------------------------
62 :
63 : #include "config-f90.h"
64 :
65 : #ifdef HAVE_SSE_INTRINSICS
66 : #include <x86intrin.h>
67 : #endif
68 : #ifdef HAVE_SPARC64_SSE
69 : #include <fjmfunc.h>
70 : #include <emmintrin.h>
71 : #endif
72 : #include <stdio.h>
73 : #include <stdlib.h>
74 :
75 : #ifdef DOUBLE_PRECISION_REAL
76 : #define offset 2
77 : #define __SSE_DATATYPE __m128d
78 : #define _SSE_LOAD _mm_load_pd
79 : #define _SSE_ADD _mm_add_pd
80 : #define _SSE_SUB _mm_sub_pd
81 : #define _SSE_MUL _mm_mul_pd
82 : #define _SSE_STORE _mm_store_pd
83 : #endif
84 : #ifdef SINGLE_PRECISION_REAL
85 : #define offset 4
86 : #define __SSE_DATATYPE __m128
87 : #define _SSE_LOAD _mm_load_ps
88 : #define _SSE_ADD _mm_add_ps
89 : #define _SSE_SUB _mm_sub_ps
90 : #define _SSE_MUL _mm_mul_ps
91 : #define _SSE_STORE _mm_store_ps
92 : #endif
93 :
94 : #define __forceinline __attribute__((always_inline)) static
95 :
96 : #ifdef HAVE_SSE_INTRINSICS
97 : #undef __AVX__
98 : #endif
99 :
100 : //Forward declaration
101 : #ifdef HAVE_SSE_INTRINSICS
102 : #ifdef DOUBLE_PRECISION_REAL
103 : __forceinline void hh_trafo_kernel_2_SSE_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4);
104 : __forceinline void hh_trafo_kernel_4_SSE_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4);
105 : __forceinline void hh_trafo_kernel_6_SSE_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4);
106 : #endif
107 : #ifdef SINGLE_PRECISION_REAL
108 : __forceinline void hh_trafo_kernel_4_SSE_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4);
109 : __forceinline void hh_trafo_kernel_8_SSE_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4);
110 : __forceinline void hh_trafo_kernel_12_SSE_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4);
111 : #endif
112 : #endif
113 : #ifdef HAVE_SPARC64_SSE
114 : #ifdef DOUBLE_PRECISION_REAL
115 : __forceinline void hh_trafo_kernel_2_SPARC64_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4);
116 : __forceinline void hh_trafo_kernel_4_SPARC64_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4);
117 : __forceinline void hh_trafo_kernel_6_SPARC64_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4);
118 : #endif
119 : #ifdef SINGLE_PRECISION_REAL
120 : __forceinline void hh_trafo_kernel_4_SPARC64_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4);
121 : __forceinline void hh_trafo_kernel_8_SPARC64_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4);
122 : __forceinline void hh_trafo_kernel_12_SPARC64_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4);
123 : #endif
124 : #endif
125 :
126 : #ifdef HAVE_SSE_INTRINSICS
127 : #ifdef DOUBLE_PRECISION_REAL
128 : void quad_hh_trafo_real_sse_4hv_double(double* q, double* hh, int* pnb, int* pnq, int* pldq, int* pldh);
129 : #endif
130 : #ifdef SINGLE_PRECISION_REAL
131 : void quad_hh_trafo_real_sse_4hv_single(float* q, float* hh, int* pnb, int* pnq, int* pldq, int* pldh);
132 : #endif
133 : #endif
134 : #ifdef HAVE_SPARC64_SSE
135 : #ifdef DOUBLE_PRECISION_REAL
136 : void quad_hh_trafo_real_sparc64_4hv_double(double* q, double* hh, int* pnb, int* pnq, int* pldq, int* pldh);
137 : #endif
138 : #ifdef SINGLE_PRECISION_REAL
139 : void quad_hh_trafo_real_sparc64_4hv_single(float* q, float* hh, int* pnb, int* pnq, int* pldq, int* pldh);
140 : #endif
141 : #endif
142 :
143 : /*
144 : !f>#ifdef HAVE_SSE_INTRINSICS
145 : !f> interface
146 : !f> subroutine quad_hh_trafo_real_sse_4hv_double(q, hh, pnb, pnq, pldq, pldh) &
147 : !f> bind(C, name="quad_hh_trafo_real_sse_4hv_double")
148 : !f> use, intrinsic :: iso_c_binding
149 : !f> integer(kind=c_int) :: pnb, pnq, pldq, pldh
150 : !f> type(c_ptr), value :: q
151 : !f> real(kind=c_double) :: hh(pnb,6)
152 : !f> end subroutine
153 : !f> end interface
154 : !f>#endif
155 : */
156 :
157 : /*
158 : !f>#ifdef HAVE_SPARC64_SSE
159 : !f> interface
160 : !f> subroutine quad_hh_trafo_real_sparc64_4hv_double(q, hh, pnb, pnq, pldq, pldh) &
161 : !f> bind(C, name="quad_hh_trafo_real_sparc64_4hv_double")
162 : !f> use, intrinsic :: iso_c_binding
163 : !f> integer(kind=c_int) :: pnb, pnq, pldq, pldh
164 : !f> type(c_ptr), value :: q
165 : !f> real(kind=c_double) :: hh(pnb,6)
166 : !f> end subroutine
167 : !f> end interface
168 : !f>#endif
169 : */
170 :
171 : /*
172 : !f>#ifdef HAVE_SSE_INTRINSICS
173 : !f> interface
174 : !f> subroutine quad_hh_trafo_real_sse_4hv_single(q, hh, pnb, pnq, pldq, pldh) &
175 : !f> bind(C, name="quad_hh_trafo_real_sse_4hv_single")
176 : !f> use, intrinsic :: iso_c_binding
177 : !f> integer(kind=c_int) :: pnb, pnq, pldq, pldh
178 : !f> type(c_ptr), value :: q
179 : !f> real(kind=c_float) :: hh(pnb,6)
180 : !f> end subroutine
181 : !f> end interface
182 : !f>#endif
183 : */
184 :
185 : /*
186 : !f>#ifdef HAVE_SPARC64_SSE
187 : !f> interface
188 : !f> subroutine quad_hh_trafo_real_sparc64_4hv_single(q, hh, pnb, pnq, pldq, pldh) &
189 : !f> bind(C, name="quad_hh_trafo_real_sparc64_4hv_single")
190 : !f> use, intrinsic :: iso_c_binding
191 : !f> integer(kind=c_int) :: pnb, pnq, pldq, pldh
192 : !f> type(c_ptr), value :: q
193 : !f> real(kind=c_float) :: hh(pnb,6)
194 : !f> end subroutine
195 : !f> end interface
196 : !f>#endif
197 : */
198 :
199 : #ifdef HAVE_SSE_INTRINSICS
200 : #ifdef DOUBLE_PRECISION_REAL
201 268288 : void quad_hh_trafo_real_sse_4hv_double(double* q, double* hh, int* pnb, int* pnq, int* pldq, int* pldh)
202 : #endif
203 : #ifdef SINGLE_PRECISION_REAL
204 50304 : void quad_hh_trafo_real_sse_4hv_single(float* q, float* hh, int* pnb, int* pnq, int* pldq, int* pldh)
205 : #endif
206 : #endif
207 : #ifdef HAVE_SPARC64_SSE
208 : #ifdef DOUBLE_PRECISION_REAL
209 : void quad_hh_trafo_real_sparc64_4hv_double(double* q, double* hh, int* pnb, int* pnq, int* pldq, int* pldh)
210 : #endif
211 : #ifdef SINGLE_PRECISION_REAL
212 : void quad_hh_trafo_real_sparc64_4hv_single(float* q, float* hh, int* pnb, int* pnq, int* pldq, int* pldh)
213 : #endif
214 :
215 : #endif
216 :
217 : {
218 : int i;
219 318592 : int nb = *pnb;
220 318592 : int nq = *pldq;
221 318592 : int ldq = *pldq;
222 318592 : int ldh = *pldh;
223 : int worked_on;
224 :
225 :
226 318592 : worked_on = 0;
227 :
228 : // calculating scalar products to compute
229 : // 4 householder vectors simultaneously
230 : #ifdef DOUBLE_PRECISION_REAL
231 268288 : double s_1_2 = hh[(ldh)+1];
232 268288 : double s_1_3 = hh[(ldh*2)+2];
233 268288 : double s_2_3 = hh[(ldh*2)+1];
234 268288 : double s_1_4 = hh[(ldh*3)+3];
235 268288 : double s_2_4 = hh[(ldh*3)+2];
236 268288 : double s_3_4 = hh[(ldh*3)+1];
237 : #endif
238 : #ifdef SINGLE_PRECISION_REAL
239 50304 : float s_1_2 = hh[(ldh)+1];
240 50304 : float s_1_3 = hh[(ldh*2)+2];
241 50304 : float s_2_3 = hh[(ldh*2)+1];
242 50304 : float s_1_4 = hh[(ldh*3)+3];
243 50304 : float s_2_4 = hh[(ldh*3)+2];
244 50304 : float s_3_4 = hh[(ldh*3)+1];
245 : #endif
246 : // calculate scalar product of first and fourth householder Vector
247 : // loop counter = 2
248 318592 : s_1_2 += hh[2-1] * hh[(2+ldh)];
249 318592 : s_2_3 += hh[(ldh)+2-1] * hh[2+(ldh*2)];
250 318592 : s_3_4 += hh[(ldh*2)+2-1] * hh[2+(ldh*3)];
251 :
252 : // loop counter = 3
253 318592 : s_1_2 += hh[3-1] * hh[(3+ldh)];
254 318592 : s_2_3 += hh[(ldh)+3-1] * hh[3+(ldh*2)];
255 318592 : s_3_4 += hh[(ldh*2)+3-1] * hh[3+(ldh*3)];
256 :
257 318592 : s_1_3 += hh[3-2] * hh[3+(ldh*2)];
258 318592 : s_2_4 += hh[(ldh*1)+3-2] * hh[3+(ldh*3)];
259 :
260 : #pragma ivdep
261 19434112 : for (i = 4; i < nb; i++)
262 : {
263 19115520 : s_1_2 += hh[i-1] * hh[(i+ldh)];
264 19115520 : s_2_3 += hh[(ldh)+i-1] * hh[i+(ldh*2)];
265 19115520 : s_3_4 += hh[(ldh*2)+i-1] * hh[i+(ldh*3)];
266 :
267 19115520 : s_1_3 += hh[i-2] * hh[i+(ldh*2)];
268 19115520 : s_2_4 += hh[(ldh*1)+i-2] * hh[i+(ldh*3)];
269 :
270 19115520 : s_1_4 += hh[i-3] * hh[i+(ldh*3)];
271 : }
272 :
273 : // Production level kernel calls with padding
274 : #ifdef DOUBLE_PRECISION_REAL
275 1878016 : for (i = 0; i < nq-4; i+=6)
276 : {
277 : #ifdef HAVE_SSE_INTRINSICS
278 1609728 : hh_trafo_kernel_6_SSE_4hv_double(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
279 : #endif
280 : #ifdef HAVE_SPARC64_SSE
281 : hh_trafo_kernel_6_SPARC64_4hv_double(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
282 : #endif
283 1609728 : worked_on += 6;
284 : }
285 : #endif
286 : #ifdef SINGLE_PRECISION_REAL
287 352128 : for (i = 0; i < nq-8; i+=12)
288 : {
289 : #ifdef HAVE_SSE_INTRINSICS
290 301824 : hh_trafo_kernel_12_SSE_4hv_single(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
291 : #endif
292 : #ifdef HAVE_SPARC64_SSE
293 : hh_trafo_kernel_12_SPARC64_4hv_single(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
294 : #endif
295 :
296 301824 : worked_on += 12;
297 :
298 : }
299 : #endif
300 318592 : if (nq == i)
301 : {
302 0 : return;
303 : }
304 :
305 : #ifdef DOUBLE_PRECISION_REAL
306 268288 : if (nq-i ==4)
307 : {
308 : #ifdef HAVE_SSE_INTRINSICS
309 268288 : hh_trafo_kernel_4_SSE_4hv_double(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
310 : #endif
311 : #ifdef HAVE_SPARC64_SSE
312 : hh_trafo_kernel_4_SPARC64_4hv_double(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
313 : #endif
314 268288 : worked_on += 4;
315 : }
316 : #endif
317 :
318 : #ifdef SINGLE_PRECISION_REAL
319 50304 : if (nq-i ==8)
320 : {
321 : #ifdef HAVE_SSE_INTRINSICS
322 50304 : hh_trafo_kernel_8_SSE_4hv_single(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
323 : #endif
324 : #ifdef HAVE_SPARC64_SSE
325 : hh_trafo_kernel_8_SPARC64_4hv_single(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
326 : #endif
327 50304 : worked_on += 8;
328 : }
329 : #endif
330 :
331 : #ifdef DOUBLE_PRECISION_REAL
332 268288 : if (nq-i == 2)
333 : {
334 : #ifdef HAVE_SSE_INTRINSICS
335 0 : hh_trafo_kernel_2_SSE_4hv_double(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
336 : #endif
337 : #ifdef HAVE_SPARC64_SSE
338 : hh_trafo_kernel_2_SPARC64_4hv_double(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
339 : #endif
340 :
341 0 : worked_on += 2;
342 : }
343 : #endif
344 :
345 : #ifdef SINGLE_PRECISION_REAL
346 50304 : if (nq-i ==4)
347 : {
348 : #ifdef HAVE_SSE_INTRINSICS
349 0 : hh_trafo_kernel_4_SSE_4hv_single(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
350 : #endif
351 : #ifdef HAVE_SPARC64_SSE
352 : hh_trafo_kernel_4_SPARC64_4hv_single(&q[i], hh, nb, ldq, ldh, s_1_2, s_1_3, s_2_3, s_1_4, s_2_4, s_3_4);
353 : #endif
354 0 : worked_on += 4;
355 : }
356 : #endif
357 : #ifdef WITH_DEBUG
358 : if (worked_on != nq)
359 : {
360 : #ifdef HAVE_SSE_INTRINSICS
361 : printf("Error in real SSE BLOCK4 kernel \n");
362 : #endif
363 : #ifdef HAVE_SPARC64_SSE
364 : printf("Error in real SPARC64 BLOCK4 kernel \n");
365 : #endif
366 :
367 : abort();
368 : }
369 : #endif
370 :
371 : }
372 :
373 : /**
374 : * Unrolled kernel that computes
375 : #ifdef DOUBLE_PRECISION_REAL
376 : * 6 rows of Q simultaneously, a
377 : #endif
378 : #ifdef SINGLE_PRECISION_REAL
379 : * 12 rows of Q simultaneously, a
380 : #endif
381 : * matrix Vector product with two householder
382 : * vectors + a rank 1 update is performed
383 : */
384 : #ifdef HAVE_SSE_INTRINSICS
385 : #ifdef DOUBLE_PRECISION_REAL
386 : __forceinline void hh_trafo_kernel_6_SSE_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4)
387 : #endif
388 : #ifdef SINGLE_PRECISION_REAL
389 : __forceinline void hh_trafo_kernel_12_SSE_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4)
390 : #endif
391 : #endif
392 : #ifdef HAVE_SPARC64_SSE
393 : #ifdef DOUBLE_PRECISION_REAL
394 : __forceinline void hh_trafo_kernel_6_SPARC64_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4)
395 : #endif
396 : #ifdef SINGLE_PRECISION_REAL
397 : __forceinline void hh_trafo_kernel_12_SPARC64_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4)
398 : #endif
399 : #endif
400 : {
401 : /////////////////////////////////////////////////////
402 : // Matrix Vector Multiplication, Q [6 x nb+3] * hh
403 : // hh contains four householder vectors
404 : /////////////////////////////////////////////////////
405 : int i;
406 :
407 3823104 : __SSE_DATATYPE a1_1 = _SSE_LOAD(&q[ldq*3]);
408 3823104 : __SSE_DATATYPE a2_1 = _SSE_LOAD(&q[ldq*2]);
409 3823104 : __SSE_DATATYPE a3_1 = _SSE_LOAD(&q[ldq]);
410 1911552 : __SSE_DATATYPE a4_1 = _SSE_LOAD(&q[0]);
411 :
412 : #ifdef HAVE_SSE_INTRINSICS
413 : #ifdef DOUBLE_PRECISION_REAL
414 3219456 : __SSE_DATATYPE h_2_1 = _mm_set1_pd(hh[ldh+1]);
415 3219456 : __SSE_DATATYPE h_3_2 = _mm_set1_pd(hh[(ldh*2)+1]);
416 3219456 : __SSE_DATATYPE h_3_1 = _mm_set1_pd(hh[(ldh*2)+2]);
417 3219456 : __SSE_DATATYPE h_4_3 = _mm_set1_pd(hh[(ldh*3)+1]);
418 3219456 : __SSE_DATATYPE h_4_2 = _mm_set1_pd(hh[(ldh*3)+2]);
419 3219456 : __SSE_DATATYPE h_4_1 = _mm_set1_pd(hh[(ldh*3)+3]);
420 : #endif
421 :
422 : #ifdef SINGLE_PRECISION_REAL
423 603648 : __m128 h_2_1 = _mm_set1_ps(hh[ldh+1] ); // h_2_1 contains four times hh[ldh+1]
424 603648 : __m128 h_3_2 = _mm_set1_ps(hh[(ldh*2)+1]);
425 603648 : __m128 h_3_1 = _mm_set1_ps(hh[(ldh*2)+2]);
426 603648 : __m128 h_4_3 = _mm_set1_ps(hh[(ldh*3)+1]);
427 603648 : __m128 h_4_2 = _mm_set1_ps(hh[(ldh*3)+2]);
428 603648 : __m128 h_4_1 = _mm_set1_ps(hh[(ldh*3)+3]);
429 : #endif
430 : #endif
431 :
432 : #ifdef HAVE_SPARC64_SSE
433 : #ifdef DOUBLE_PRECISION_REAL
434 : __SSE_DATATYPE h_2_1 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
435 : __SSE_DATATYPE h_3_2 = _mm_set_pd(hh[(ldh*2)+1], hh[(ldh*2)+1]);
436 : __SSE_DATATYPE h_3_1 = _mm_set_pd(hh[(ldh*2)+2], hh[(ldh*2)+2]);
437 : __SSE_DATATYPE h_4_3 = _mm_set_pd(hh[(ldh*3)+1], hh[(ldh*3)+1]);
438 : __SSE_DATATYPE h_4_2 = _mm_set_pd(hh[(ldh*3)+2], hh[(ldh*3)+2]);
439 : __SSE_DATATYPE h_4_1 = _mm_set_pd(hh[(ldh*3)+3], hh[(ldh*3)+3]);
440 : #endif
441 :
442 : #ifdef SINGLE_PRECISION_REAL
443 : __m128 h_2_1 = _mm_set_ps(hh[ldh+1], hh[ldh+1]); // h_2_1 contains four times hh[ldh+1]
444 : __m128 h_3_2 = _mm_set_ps(hh[(ldh*2)+1], hh[(ldh*2)+1]);
445 : __m128 h_3_1 = _mm_set_ps(hh[(ldh*2)+2], hh[(ldh*2)+2]);
446 : __m128 h_4_3 = _mm_set_ps(hh[(ldh*3)+1], hh[(ldh*3)+1]);
447 : __m128 h_4_2 = _mm_set_ps(hh[(ldh*3)+2], hh[(ldh*3)+2]);
448 : __m128 h_4_1 = _mm_set_ps(hh[(ldh*3)+3], hh[(ldh*3)+3]);
449 : #endif
450 : #endif
451 :
452 :
453 :
454 3823104 : register __SSE_DATATYPE w1 = _SSE_ADD(a4_1, _SSE_MUL(a3_1, h_4_3));
455 3823104 : w1 = _SSE_ADD(w1, _SSE_MUL(a2_1, h_4_2));
456 3823104 : w1 = _SSE_ADD(w1, _SSE_MUL(a1_1, h_4_1));
457 3823104 : register __SSE_DATATYPE z1 = _SSE_ADD(a3_1, _SSE_MUL(a2_1, h_3_2));
458 3823104 : z1 = _SSE_ADD(z1, _SSE_MUL(a1_1, h_3_1));
459 3823104 : register __SSE_DATATYPE y1 = _SSE_ADD(a2_1, _SSE_MUL(a1_1, h_2_1));
460 1911552 : register __SSE_DATATYPE x1 = a1_1;
461 :
462 3823104 : __SSE_DATATYPE a1_2 = _SSE_LOAD(&q[(ldq*3)+offset]);
463 3823104 : __SSE_DATATYPE a2_2 = _SSE_LOAD(&q[(ldq*2)+offset]);
464 3823104 : __SSE_DATATYPE a3_2 = _SSE_LOAD(&q[ldq+offset]);
465 3823104 : __SSE_DATATYPE a4_2 = _SSE_LOAD(&q[0+offset]);
466 :
467 3823104 : register __SSE_DATATYPE w2 = _SSE_ADD(a4_2, _SSE_MUL(a3_2, h_4_3));
468 3823104 : w2 = _SSE_ADD(w2, _SSE_MUL(a2_2, h_4_2));
469 3823104 : w2 = _SSE_ADD(w2, _SSE_MUL(a1_2, h_4_1));
470 3823104 : register __SSE_DATATYPE z2 = _SSE_ADD(a3_2, _SSE_MUL(a2_2, h_3_2));
471 3823104 : z2 = _SSE_ADD(z2, _SSE_MUL(a1_2, h_3_1));
472 3823104 : register __SSE_DATATYPE y2 = _SSE_ADD(a2_2, _SSE_MUL(a1_2, h_2_1));
473 1911552 : register __SSE_DATATYPE x2 = a1_2;
474 :
475 3823104 : __SSE_DATATYPE a1_3 = _SSE_LOAD(&q[(ldq*3)+2*offset]);
476 3823104 : __SSE_DATATYPE a2_3 = _SSE_LOAD(&q[(ldq*2)+2*offset]);
477 3823104 : __SSE_DATATYPE a3_3 = _SSE_LOAD(&q[ldq+2*offset]);
478 3823104 : __SSE_DATATYPE a4_3 = _SSE_LOAD(&q[0+2*offset]);
479 :
480 3823104 : register __SSE_DATATYPE w3 = _SSE_ADD(a4_3, _SSE_MUL(a3_3, h_4_3));
481 3823104 : w3 = _SSE_ADD(w3, _SSE_MUL(a2_3, h_4_2));
482 3823104 : w3 = _SSE_ADD(w3, _SSE_MUL(a1_3, h_4_1));
483 3823104 : register __SSE_DATATYPE z3 = _SSE_ADD(a3_3, _SSE_MUL(a2_3, h_3_2));
484 3823104 : z3 = _SSE_ADD(z3, _SSE_MUL(a1_3, h_3_1));
485 3823104 : register __SSE_DATATYPE y3 = _SSE_ADD(a2_3, _SSE_MUL(a1_3, h_2_1));
486 1911552 : register __SSE_DATATYPE x3 = a1_3;
487 :
488 : __SSE_DATATYPE q1;
489 : __SSE_DATATYPE q2;
490 : __SSE_DATATYPE q3;
491 :
492 : __SSE_DATATYPE h1;
493 : __SSE_DATATYPE h2;
494 : __SSE_DATATYPE h3;
495 : __SSE_DATATYPE h4;
496 :
497 116604672 : for(i = 4; i < nb; i++)
498 : {
499 : #ifdef HAVE_SSE_INTRINSICS
500 : #ifdef DOUBLE_PRECISION_REAL
501 193167360 : h1 = _mm_set1_pd(hh[i-3]);
502 : #endif
503 : #ifdef SINGLE_PRECISION_REAL
504 36218880 : h1 = _mm_set1_ps(hh[i-3]);
505 : #endif
506 : #endif
507 :
508 : #ifdef HAVE_SPARC64_SSE
509 : #ifdef DOUBLE_PRECISION_REAL
510 : h1 = _mm_set_pd(hh[i-3], hh[i-3]);
511 : #endif
512 : #ifdef SINGLE_PRECISION_REAL
513 : h1 = _mm_set_ps(hh[i-3], hh[i-3]);
514 : #endif
515 : #endif
516 229386240 : q1 = _SSE_LOAD(&q[i*ldq]);
517 229386240 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
518 229386240 : q3 = _SSE_LOAD(&q[(i*ldq)+2*offset]);
519 :
520 229386240 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
521 229386240 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
522 229386240 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
523 :
524 : #ifdef HAVE_SSE_INTRINSICS
525 : #ifdef DOUBLE_PRECISION_REAL
526 193167360 : h2 = _mm_set1_pd(hh[ldh+i-2]);
527 : #endif
528 : #ifdef SINGLE_PRECISION_REAL
529 36218880 : h2 = _mm_set1_ps(hh[ldh+i-2]);
530 : #endif
531 : #endif
532 :
533 : #ifdef HAVE_SPARC64_SSE
534 : #ifdef DOUBLE_PRECISION_REAL
535 : h2 = _mm_set_pd(hh[ldh+i-2], hh[ldh+i-2]);
536 : #endif
537 : #ifdef SINGLE_PRECISION_REAL
538 : h2 = _mm_set_ps(hh[ldh+i-2], hh[ldh+i-2]);
539 : #endif
540 :
541 : #endif
542 229386240 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
543 229386240 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
544 229386240 : y3 = _SSE_ADD(y3, _SSE_MUL(q3,h2));
545 :
546 : #ifdef HAVE_SSE_INTRINSICS
547 : #ifdef DOUBLE_PRECISION_REAL
548 193167360 : h3 = _mm_set1_pd(hh[(ldh*2)+i-1]);
549 : #endif
550 : #ifdef SINGLE_PRECISION_REAL
551 36218880 : h3 = _mm_set1_ps(hh[(ldh*2)+i-1]);
552 : #endif
553 : #endif
554 :
555 : #ifdef HAVE_SPARC64_SSE
556 : #ifdef DOUBLE_PRECISION_REAL
557 : h3 = _mm_set_pd(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
558 : #endif
559 : #ifdef SINGLE_PRECISION_REAL
560 : h3 = _mm_set_ps(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
561 : #endif
562 : #endif
563 :
564 :
565 229386240 : z1 = _SSE_ADD(z1, _SSE_MUL(q1,h3));
566 229386240 : z2 = _SSE_ADD(z2, _SSE_MUL(q2,h3));
567 229386240 : z3 = _SSE_ADD(z3, _SSE_MUL(q3,h3));
568 : #ifdef HAVE_SSE_INTRINSICS
569 : #ifdef DOUBLE_PRECISION_REAL
570 193167360 : h4 = _mm_set1_pd(hh[(ldh*3)+i]);
571 : #endif
572 : #ifdef SINGLE_PRECISION_REAL
573 36218880 : h4 = _mm_set1_ps(hh[(ldh*3)+i]);
574 : #endif
575 : #endif
576 :
577 : #ifdef HAVE_SPARC64_SSE
578 : #ifdef DOUBLE_PRECISION_REAL
579 : h4 = _mm_set_pd(hh[(ldh*3)+i], hh[(ldh*3)+i]);
580 : #endif
581 : #ifdef SINGLE_PRECISION_REAL
582 : h4 = _mm_set_ps(hh[(ldh*3)+i], hh[(ldh*3)+i]);
583 : #endif
584 : #endif
585 :
586 229386240 : w1 = _SSE_ADD(w1, _SSE_MUL(q1,h4));
587 229386240 : w2 = _SSE_ADD(w2, _SSE_MUL(q2,h4));
588 229386240 : w3 = _SSE_ADD(w3, _SSE_MUL(q3,h4));
589 : }
590 :
591 : #ifdef HAVE_SSE_INTRINSICS
592 : #ifdef DOUBLE_PRECISION_REAL
593 3219456 : h1 = _mm_set1_pd(hh[nb-3]);
594 : #endif
595 : #ifdef SINGLE_PRECISION_REAL
596 603648 : h1 = _mm_set1_ps(hh[nb-3]);
597 : #endif
598 : #endif
599 :
600 : #ifdef HAVE_SPARC64_SSE
601 : #ifdef DOUBLE_PRECISION_REAL
602 : h1 = _mm_set_pd(hh[nb-3], hh[nb-3]);
603 : #endif
604 : #ifdef SINGLE_PRECISION_REAL
605 : h1 = _mm_set_ps(hh[nb-3], hh[nb-3]);
606 : #endif
607 : #endif
608 3823104 : q1 = _SSE_LOAD(&q[nb*ldq]);
609 3823104 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
610 3823104 : q3 = _SSE_LOAD(&q[(nb*ldq)+2*offset]);
611 :
612 3823104 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
613 3823104 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
614 3823104 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
615 :
616 : #ifdef HAVE_SSE_INTRINSICS
617 : #ifdef DOUBLE_PRECISION_REAL
618 3219456 : h2 = _mm_set1_pd(hh[ldh+nb-2]);
619 : #endif
620 : #ifdef SINGLE_PRECISION_REAL
621 603648 : h2 = _mm_set1_ps(hh[ldh+nb-2]);
622 : #endif
623 : #endif
624 :
625 : #ifdef HAVE_SPARC64_SSE
626 : #ifdef DOUBLE_PRECISION_REAL
627 : h2 = _mm_set_pd(hh[ldh+nb-2], hh[ldh+nb-2]);
628 : #endif
629 : #ifdef SINGLE_PRECISION_REAL
630 : h2 = _mm_set_ps(hh[ldh+nb-2], hh[ldh+nb-2]);
631 : #endif
632 : #endif
633 :
634 3823104 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
635 3823104 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
636 3823104 : y3 = _SSE_ADD(y3, _SSE_MUL(q3,h2));
637 :
638 : #ifdef HAVE_SSE_INTRINSICS
639 : #ifdef DOUBLE_PRECISION_REAL
640 3219456 : h3 = _mm_set1_pd(hh[(ldh*2)+nb-1]);
641 : #endif
642 : #ifdef SINGLE_PRECISION_REAL
643 603648 : h3 = _mm_set1_ps(hh[(ldh*2)+nb-1]);
644 : #endif
645 : #endif
646 :
647 : #ifdef HAVE_SPARC64_SSE
648 : #ifdef DOUBLE_PRECISION_REAL
649 : h3 = _mm_set_pd(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
650 : #endif
651 : #ifdef SINGLE_PRECISION_REAL
652 : h3 = _mm_set_ps(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
653 : #endif
654 : #endif
655 :
656 3823104 : z1 = _SSE_ADD(z1, _SSE_MUL(q1,h3));
657 3823104 : z2 = _SSE_ADD(z2, _SSE_MUL(q2,h3));
658 3823104 : z3 = _SSE_ADD(z3, _SSE_MUL(q3,h3));
659 :
660 : #ifdef HAVE_SSE_INTRINSICS
661 : #ifdef DOUBLE_PRECISION_REAL
662 3219456 : h1 = _mm_set1_pd(hh[nb-2]);
663 : #endif
664 : #ifdef SINGLE_PRECISION_REAL
665 603648 : h1 = _mm_set1_ps(hh[nb-2]);
666 : #endif
667 : #endif
668 :
669 : #ifdef HAVE_SPARC64_SSE
670 : #ifdef DOUBLE_PRECISION_REAL
671 : h1 = _mm_set_pd(hh[nb-2], hh[nb-2]);
672 : #endif
673 : #ifdef SINGLE_PRECISION_REAL
674 : h1 = _mm_set_ps(hh[nb-2], hh[nb-2]);
675 : #endif
676 : #endif
677 :
678 3823104 : q1 = _SSE_LOAD(&q[(nb+1)*ldq]);
679 3823104 : q2 = _SSE_LOAD(&q[((nb+1)*ldq)+offset]);
680 3823104 : q3 = _SSE_LOAD(&q[((nb+1)*ldq)+2*offset]);
681 :
682 3823104 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
683 3823104 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
684 3823104 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
685 :
686 : #ifdef HAVE_SSE_INTRINSICS
687 : #ifdef DOUBLE_PRECISION_REAL
688 3219456 : h2 = _mm_set1_pd(hh[(ldh*1)+nb-1]);
689 : #endif
690 : #ifdef SINGLE_PRECISION_REAL
691 603648 : h2 = _mm_set1_ps(hh[ldh+nb-1]);
692 : #endif
693 : #endif
694 :
695 : #ifdef HAVE_SPARC64_SSE
696 : #ifdef DOUBLE_PRECISION_REAL
697 : h2 = _mm_set_pd(hh[(ldh*1)+nb-1], hh[(ldh*1)+nb-1]);
698 : #endif
699 : #ifdef SINGLE_PRECISION_REAL
700 : h2 = _mm_set_ps(hh[ldh+nb-1], hh[(ldh*1)+nb-1]);
701 : #endif
702 : #endif
703 :
704 :
705 3823104 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
706 3823104 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
707 3823104 : y3 = _SSE_ADD(y3, _SSE_MUL(q3,h2));
708 :
709 : #ifdef HAVE_SSE_INTRINSICS
710 : #ifdef DOUBLE_PRECISION_REAL
711 3219456 : h1 = _mm_set1_pd(hh[nb-1]);
712 : #endif
713 : #ifdef SINGLE_PRECISION_REAL
714 603648 : h1 = _mm_set1_ps(hh[nb-1]);
715 : #endif
716 : #endif
717 :
718 : #ifdef HAVE_SPARC64_SSE
719 : #ifdef DOUBLE_PRECISION_REAL
720 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
721 : #endif
722 : #ifdef SINGLE_PRECISION_REAL
723 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
724 : #endif
725 : #endif
726 :
727 :
728 3823104 : q1 = _SSE_LOAD(&q[(nb+2)*ldq]);
729 3823104 : q2 = _SSE_LOAD(&q[((nb+2)*ldq)+offset]);
730 3823104 : q3 = _SSE_LOAD(&q[((nb+2)*ldq)+2*offset]);
731 :
732 3823104 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
733 3823104 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
734 3823104 : x3 = _SSE_ADD(x3, _SSE_MUL(q3,h1));
735 :
736 : /////////////////////////////////////////////////////
737 : // Rank-1 update of Q [6 x nb+3]
738 : /////////////////////////////////////////////////////
739 :
740 : #ifdef HAVE_SSE_INTRINSICS
741 : #ifdef DOUBLE_PRECISION_REAL
742 3219456 : __SSE_DATATYPE tau1 = _mm_set1_pd(hh[0]);
743 : #endif
744 : #ifdef SINGLE_PRECISION_REAL
745 603648 : __m128 tau1 = _mm_set1_ps(hh[0]);
746 : #endif
747 : #endif
748 :
749 : #ifdef HAVE_SPARC64_SSE
750 : #ifdef DOUBLE_PRECISION_REAL
751 : __SSE_DATATYPE tau1 = _mm_set_pd(hh[0], hh[0]);
752 : #endif
753 : #ifdef SINGLE_PRECISION_REAL
754 : __m128 tau1 = _mm_set_ps(hh[0], hh[0]);
755 : #endif
756 : #endif
757 :
758 :
759 1911552 : h1 = tau1;
760 1911552 : x1 = _SSE_MUL(x1, h1);
761 1911552 : x2 = _SSE_MUL(x2, h1);
762 1911552 : x3 = _SSE_MUL(x3, h1);
763 :
764 : #ifdef HAVE_SSE_INTRINSICS
765 : #ifdef DOUBLE_PRECISION_REAL
766 3219456 : __SSE_DATATYPE tau2 = _mm_set1_pd(hh[ldh]);
767 1609728 : __SSE_DATATYPE vs_1_2 = _mm_set1_pd(s_1_2);
768 : #endif
769 : #ifdef SINGLE_PRECISION_REAL
770 603648 : __m128 tau2 = _mm_set1_ps(hh[ldh]);
771 301824 : __m128 vs_1_2 = _mm_set1_ps(s_1_2);
772 : #endif
773 : #endif
774 :
775 : #ifdef HAVE_SPARC64_SSE
776 : #ifdef DOUBLE_PRECISION_REAL
777 : __SSE_DATATYPE tau2 = _mm_set_pd(hh[ldh], hh[ldh]);
778 : __SSE_DATATYPE vs_1_2 = _mm_set_pd(s_1_2, s_1_2);
779 : #endif
780 : #ifdef SINGLE_PRECISION_REAL
781 : __m128 tau2 = _mm_set_ps(hh[ldh], hh[ldh]);
782 : __m128 vs_1_2 = _mm_set_ps(s_1_2, s_1_2);
783 : #endif
784 : #endif
785 :
786 :
787 1911552 : h1 = tau2;
788 1911552 : h2 = _SSE_MUL(h1, vs_1_2);
789 :
790 5734656 : y1 = _SSE_SUB(_SSE_MUL(y1,h1), _SSE_MUL(x1,h2));
791 5734656 : y2 = _SSE_SUB(_SSE_MUL(y2,h1), _SSE_MUL(x2,h2));
792 5734656 : y3 = _SSE_SUB(_SSE_MUL(y3,h1), _SSE_MUL(x3,h2));
793 :
794 : #ifdef HAVE_SSE_INTRINSICS
795 : #ifdef DOUBLE_PRECISION_REAL
796 3219456 : __SSE_DATATYPE tau3 = _mm_set1_pd(hh[ldh*2]);
797 1609728 : __SSE_DATATYPE vs_1_3 = _mm_set1_pd(s_1_3);
798 1609728 : __SSE_DATATYPE vs_2_3 = _mm_set1_pd(s_2_3);
799 : #endif
800 : #ifdef SINGLE_PRECISION_REAL
801 603648 : __m128 tau3 = _mm_set1_ps(hh[ldh*2]);
802 301824 : __m128 vs_1_3 = _mm_set1_ps(s_1_3);
803 301824 : __m128 vs_2_3 = _mm_set1_ps(s_2_3);
804 : #endif
805 : #endif
806 :
807 : #ifdef HAVE_SPARC64_SSE
808 : #ifdef DOUBLE_PRECISION_REAL
809 : __SSE_DATATYPE tau3 = _mm_set_pd(hh[ldh*2], hh[ldh*2]);
810 : __SSE_DATATYPE vs_1_3 = _mm_set_pd(s_1_3, s_1_3);
811 : __SSE_DATATYPE vs_2_3 = _mm_set_pd(s_2_3, s_2_3);
812 : #endif
813 : #ifdef SINGLE_PRECISION_REAL
814 : __m128 tau3 = _mm_set_ps(hh[ldh*2], hh[ldh*2]);
815 : __m128 vs_1_3 = _mm_set_ps(s_1_3, s_1_3);
816 : __m128 vs_2_3 = _mm_set_ps(s_2_3, s_2_3);
817 : #endif
818 : #endif
819 :
820 :
821 1911552 : h1 = tau3;
822 1911552 : h2 = _SSE_MUL(h1, vs_1_3);
823 1911552 : h3 = _SSE_MUL(h1, vs_2_3);
824 :
825 9557760 : z1 = _SSE_SUB(_SSE_MUL(z1,h1), _SSE_ADD(_SSE_MUL(y1,h3), _SSE_MUL(x1,h2)));
826 9557760 : z2 = _SSE_SUB(_SSE_MUL(z2,h1), _SSE_ADD(_SSE_MUL(y2,h3), _SSE_MUL(x2,h2)));
827 9557760 : z3 = _SSE_SUB(_SSE_MUL(z3,h1), _SSE_ADD(_SSE_MUL(y3,h3), _SSE_MUL(x3,h2)));
828 :
829 : #ifdef HAVE_SSE_INTRINSICS
830 : #ifdef DOUBLE_PRECISION_REAL
831 3219456 : __SSE_DATATYPE tau4 = _mm_set1_pd(hh[ldh*3]);
832 1609728 : __SSE_DATATYPE vs_1_4 = _mm_set1_pd(s_1_4);
833 1609728 : __SSE_DATATYPE vs_2_4 = _mm_set1_pd(s_2_4);
834 1609728 : __SSE_DATATYPE vs_3_4 = _mm_set1_pd(s_3_4);
835 : #endif
836 : #ifdef SINGLE_PRECISION_REAL
837 603648 : __m128 tau4 = _mm_set1_ps(hh[ldh*3]);
838 301824 : __m128 vs_1_4 = _mm_set1_ps(s_1_4);
839 301824 : __m128 vs_2_4 = _mm_set1_ps(s_2_4);
840 301824 : __m128 vs_3_4 = _mm_set1_ps(s_3_4);
841 : #endif
842 : #endif
843 :
844 : #ifdef HAVE_SPARC64_SSE
845 : #ifdef DOUBLE_PRECISION_REAL
846 : __SSE_DATATYPE tau4 = _mm_set_pd(hh[ldh*3], hh[ldh*3]);
847 : __SSE_DATATYPE vs_1_4 = _mm_set_pd(s_1_4, s_1_4);
848 : __SSE_DATATYPE vs_2_4 = _mm_set_pd(s_2_4, s_2_4);
849 : __SSE_DATATYPE vs_3_4 = _mm_set_pd(s_3_4, s_3_4);
850 : #endif
851 : #ifdef SINGLE_PRECISION_REAL
852 : __m128 tau4 = _mm_set_ps(hh[ldh*3], hh[ldh*3]);
853 : __m128 vs_1_4 = _mm_set_ps(s_1_4, s_1_4);
854 : __m128 vs_2_4 = _mm_set_ps(s_2_4, s_2_4);
855 : __m128 vs_3_4 = _mm_set_ps(s_3_4, s_3_4);
856 : #endif
857 : #endif
858 :
859 1911552 : h1 = tau4;
860 1911552 : h2 = _SSE_MUL(h1, vs_1_4);
861 1911552 : h3 = _SSE_MUL(h1, vs_2_4);
862 1911552 : h4 = _SSE_MUL(h1, vs_3_4);
863 :
864 13380864 : w1 = _SSE_SUB(_SSE_MUL(w1,h1), _SSE_ADD(_SSE_MUL(z1,h4), _SSE_ADD(_SSE_MUL(y1,h3), _SSE_MUL(x1,h2))));
865 13380864 : w2 = _SSE_SUB(_SSE_MUL(w2,h1), _SSE_ADD(_SSE_MUL(z2,h4), _SSE_ADD(_SSE_MUL(y2,h3), _SSE_MUL(x2,h2))));
866 13380864 : w3 = _SSE_SUB(_SSE_MUL(w3,h1), _SSE_ADD(_SSE_MUL(z3,h4), _SSE_ADD(_SSE_MUL(y3,h3), _SSE_MUL(x3,h2))));
867 :
868 1911552 : q1 = _SSE_LOAD(&q[0]);
869 3823104 : q2 = _SSE_LOAD(&q[offset]);
870 3823104 : q3 = _SSE_LOAD(&q[2*offset]);
871 1911552 : q1 = _SSE_SUB(q1, w1);
872 1911552 : q2 = _SSE_SUB(q2, w2);
873 1911552 : q3 = _SSE_SUB(q3, w3);
874 : _SSE_STORE(&q[0],q1);
875 1911552 : _SSE_STORE(&q[offset],q2);
876 1911552 : _SSE_STORE(&q[2*offset],q3);
877 :
878 : #ifdef HAVE_SSE_INTRINSICS
879 : #ifdef DOUBLE_PRECISION_REAL
880 3219456 : h4 = _mm_set1_pd(hh[(ldh*3)+1]);
881 : #endif
882 : #ifdef SINGLE_PRECISION_REAL
883 603648 : h4 = _mm_set1_ps(hh[(ldh*3)+1]);
884 : #endif
885 : #endif
886 :
887 : #ifdef HAVE_SPARC64_SSE
888 : #ifdef DOUBLE_PRECISION_REAL
889 : h4 = _mm_set_pd(hh[(ldh*3)+1], hh[(ldh*3)+1]);
890 : #endif
891 : #ifdef SINGLE_PRECISION_REAL
892 : h4 = _mm_set_ps(hh[(ldh*3)+1], hh[(ldh*3)+1]);
893 : #endif
894 : #endif
895 :
896 :
897 3823104 : q1 = _SSE_LOAD(&q[ldq]);
898 3823104 : q2 = _SSE_LOAD(&q[ldq+offset]);
899 3823104 : q3 = _SSE_LOAD(&q[ldq+2*offset]);
900 :
901 5734656 : q1 = _SSE_SUB(q1, _SSE_ADD(z1, _SSE_MUL(w1, h4)));
902 5734656 : q2 = _SSE_SUB(q2, _SSE_ADD(z2, _SSE_MUL(w2, h4)));
903 5734656 : q3 = _SSE_SUB(q3, _SSE_ADD(z3, _SSE_MUL(w3, h4)));
904 :
905 1911552 : _SSE_STORE(&q[ldq],q1);
906 1911552 : _SSE_STORE(&q[ldq+offset],q2);
907 1911552 : _SSE_STORE(&q[ldq+2*offset],q3);
908 :
909 : #ifdef HAVE_SSE_INTRINSICS
910 : #ifdef DOUBLE_PRECISION_REAL
911 3219456 : h4 = _mm_set1_pd(hh[(ldh*3)+2]);
912 : #endif
913 : #ifdef SINGLE_PRECISION_REAL
914 603648 : h4 = _mm_set1_ps(hh[(ldh*3)+2]);
915 : #endif
916 : #endif
917 :
918 : #ifdef HAVE_SPARC64_SSE
919 : #ifdef DOUBLE_PRECISION_REAL
920 : h4 = _mm_set_pd(hh[(ldh*3)+2], hh[(ldh*3)+2]);
921 : #endif
922 : #ifdef SINGLE_PRECISION_REAL
923 : h4 = _mm_set_ps(hh[(ldh*3)+2], hh[(ldh*3)+2]);
924 : #endif
925 : #endif
926 :
927 3823104 : q1 = _SSE_LOAD(&q[ldq*2]);
928 3823104 : q2 = _SSE_LOAD(&q[(ldq*2)+offset]);
929 3823104 : q3 = _SSE_LOAD(&q[(ldq*2)+2*offset]);
930 1911552 : q1 = _SSE_SUB(q1, y1);
931 1911552 : q2 = _SSE_SUB(q2, y2);
932 1911552 : q3 = _SSE_SUB(q3, y3);
933 :
934 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(w1, h4));
935 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(w2, h4));
936 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(w3, h4));
937 :
938 : #ifdef HAVE_SSE_INTRINSICS
939 : #ifdef DOUBLE_PRECISION_REAL
940 3219456 : h3 = _mm_set1_pd(hh[(ldh*2)+1]);
941 : #endif
942 : #ifdef SINGLE_PRECISION_REAL
943 603648 : h3 = _mm_set1_ps(hh[(ldh*2)+1]);
944 : #endif
945 : #endif
946 :
947 : #ifdef HAVE_SPARC64_SSE
948 : #ifdef DOUBLE_PRECISION_REAL
949 : h3 = _mm_set_pd(hh[(ldh*2)+1], hh[(ldh*2)+1]);
950 : #endif
951 : #ifdef SINGLE_PRECISION_REAL
952 : h3 = _mm_set_ps(hh[(ldh*2)+1], hh[(ldh*2)+1]);
953 : #endif
954 : #endif
955 :
956 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(z1, h3));
957 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(z2, h3));
958 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(z3, h3));
959 :
960 1911552 : _SSE_STORE(&q[ldq*2],q1);
961 1911552 : _SSE_STORE(&q[(ldq*2)+offset],q2);
962 1911552 : _SSE_STORE(&q[(ldq*2)+2*offset],q3);
963 :
964 : #ifdef HAVE_SSE_INTRINSICS
965 : #ifdef DOUBLE_PRECISION_REAL
966 3219456 : h4 = _mm_set1_pd(hh[(ldh*3)+3]);
967 : #endif
968 : #ifdef SINGLE_PRECISION_REAL
969 603648 : h4 = _mm_set1_ps(hh[(ldh*3)+3]);
970 : #endif
971 : #endif
972 :
973 : #ifdef HAVE_SPARC64_SSE
974 : #ifdef DOUBLE_PRECISION_REAL
975 : h4 = _mm_set_pd(hh[(ldh*3)+3], hh[(ldh*3)+3]);
976 : #endif
977 : #ifdef SINGLE_PRECISION_REAL
978 : h4 = _mm_set_ps(hh[(ldh*3)+3], hh[(ldh*3)+3]);
979 : #endif
980 : #endif
981 :
982 3823104 : q1 = _SSE_LOAD(&q[ldq*3]);
983 3823104 : q2 = _SSE_LOAD(&q[(ldq*3)+offset]);
984 3823104 : q3 = _SSE_LOAD(&q[(ldq*3)+2*offset]);
985 1911552 : q1 = _SSE_SUB(q1, x1);
986 1911552 : q2 = _SSE_SUB(q2, x2);
987 1911552 : q3 = _SSE_SUB(q3, x3);
988 :
989 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(w1, h4));
990 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(w2, h4));
991 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(w3, h4));
992 :
993 : #ifdef HAVE_SSE_INTRINSICS
994 : #ifdef DOUBLE_PRECISION_REAL
995 3219456 : h2 = _mm_set1_pd(hh[ldh+1]);
996 : #endif
997 : #ifdef SINGLE_PRECISION_REAL
998 603648 : h2 = _mm_set1_ps(hh[ldh+1]);
999 : #endif
1000 : #endif
1001 :
1002 : #ifdef HAVE_SPARC64_SSE
1003 : #ifdef DOUBLE_PRECISION_REAL
1004 : h2 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1005 : #endif
1006 : #ifdef SINGLE_PRECISION_REAL
1007 : h2 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
1008 : #endif
1009 : #endif
1010 :
1011 :
1012 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(y1, h2));
1013 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(y2, h2));
1014 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(y3, h2));
1015 :
1016 : #ifdef HAVE_SSE_INTRINSICS
1017 : #ifdef DOUBLE_PRECISION_REAL
1018 3219456 : h3 = _mm_set1_pd(hh[(ldh*2)+2]);
1019 : #endif
1020 : #ifdef SINGLE_PRECISION_REAL
1021 603648 : h3 = _mm_set1_ps(hh[(ldh*2)+2]);
1022 : #endif
1023 : #endif
1024 :
1025 : #ifdef HAVE_SPARC64_SSE
1026 : #ifdef DOUBLE_PRECISION_REAL
1027 : h3 = _mm_set_pd(hh[(ldh*2)+2], hh[(ldh*2)+2]);
1028 : #endif
1029 : #ifdef SINGLE_PRECISION_REAL
1030 : h3 = _mm_set_ps(hh[(ldh*2)+2], hh[(ldh*2)+2]);
1031 : #endif
1032 : #endif
1033 :
1034 :
1035 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(z1, h3));
1036 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(z2, h3));
1037 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(z3, h3));
1038 1911552 : _SSE_STORE(&q[ldq*3], q1);
1039 1911552 : _SSE_STORE(&q[(ldq*3)+offset], q2);
1040 1911552 : _SSE_STORE(&q[(ldq*3)+2*offset], q3);
1041 :
1042 116604672 : for (i = 4; i < nb; i++)
1043 : {
1044 : #ifdef HAVE_SSE_INTRINSICS
1045 : #ifdef DOUBLE_PRECISION_REAL
1046 193167360 : h1 = _mm_set1_pd(hh[i-3]);
1047 : #endif
1048 : #ifdef SINGLE_PRECISION_REAL
1049 36218880 : h1 = _mm_set1_ps(hh[i-3]);
1050 : #endif
1051 : #endif
1052 :
1053 : #ifdef HAVE_SPARC64_SSE
1054 : #ifdef DOUBLE_PRECISION_REAL
1055 : h1 = _mm_set_pd(hh[i-3], hh[i-3]);
1056 : #endif
1057 : #ifdef SINGLE_PRECISION_REAL
1058 : h1 = _mm_set_ps(hh[i-3], hh[i-3]);
1059 : #endif
1060 : #endif
1061 :
1062 229386240 : q1 = _SSE_LOAD(&q[i*ldq]);
1063 229386240 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
1064 229386240 : q3 = _SSE_LOAD(&q[(i*ldq)+2*offset]);
1065 :
1066 229386240 : q1 = _SSE_SUB(q1, _SSE_MUL(x1,h1));
1067 229386240 : q2 = _SSE_SUB(q2, _SSE_MUL(x2,h1));
1068 229386240 : q3 = _SSE_SUB(q3, _SSE_MUL(x3,h1));
1069 :
1070 : #ifdef HAVE_SSE_INTRINSICS
1071 : #ifdef DOUBLE_PRECISION_REAL
1072 193167360 : h2 = _mm_set1_pd(hh[ldh+i-2]);
1073 : #endif
1074 : #ifdef SINGLE_PRECISION_REAL
1075 36218880 : h2 = _mm_set1_ps(hh[ldh+i-2]);
1076 : #endif
1077 : #endif
1078 :
1079 : #ifdef HAVE_SPARC64_SSE
1080 : #ifdef DOUBLE_PRECISION_REAL
1081 : h2 = _mm_set_pd(hh[ldh+i-2], hh[ldh+i-2]);
1082 : #endif
1083 : #ifdef SINGLE_PRECISION_REAL
1084 : h2 = _mm_set_ps(hh[ldh+i-2], hh[ldh+i-2]);
1085 : #endif
1086 : #endif
1087 :
1088 229386240 : q1 = _SSE_SUB(q1, _SSE_MUL(y1,h2));
1089 229386240 : q2 = _SSE_SUB(q2, _SSE_MUL(y2,h2));
1090 229386240 : q3 = _SSE_SUB(q3, _SSE_MUL(y3,h2));
1091 :
1092 : #ifdef HAVE_SSE_INTRINSICS
1093 : #ifdef DOUBLE_PRECISION_REAL
1094 193167360 : h3 = _mm_set1_pd(hh[(ldh*2)+i-1]);
1095 : #endif
1096 : #ifdef SINGLE_PRECISION_REAL
1097 36218880 : h3 = _mm_set1_ps(hh[(ldh*2)+i-1]);
1098 : #endif
1099 : #endif
1100 :
1101 : #ifdef HAVE_SPARC64_SSE
1102 : #ifdef DOUBLE_PRECISION_REAL
1103 : h3 = _mm_set_pd(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
1104 : #endif
1105 : #ifdef SINGLE_PRECISION_REAL
1106 : h3 = _mm_set_ps(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
1107 : #endif
1108 : #endif
1109 :
1110 229386240 : q1 = _SSE_SUB(q1, _SSE_MUL(z1,h3));
1111 229386240 : q2 = _SSE_SUB(q2, _SSE_MUL(z2,h3));
1112 229386240 : q3 = _SSE_SUB(q3, _SSE_MUL(z3,h3));
1113 :
1114 : #ifdef HAVE_SSE_INTRINSICS
1115 : #ifdef DOUBLE_PRECISION_REAL
1116 193167360 : h4 = _mm_set1_pd(hh[(ldh*3)+i]);
1117 : #endif
1118 : #ifdef SINGLE_PRECISION_REAL
1119 36218880 : h4 = _mm_set1_ps(hh[(ldh*3)+i]);
1120 : #endif
1121 : #endif
1122 :
1123 : #ifdef HAVE_SPRC64_SSE
1124 : #ifdef DOUBLE_PRECISION_REAL
1125 : h4 = _mm_set_pd(hh[(ldh*3)+i], hh[(ldh*3)+i]);
1126 : #endif
1127 : #ifdef SINGLE_PRECISION_REAL
1128 : h4 = _mm_set_ps(hh[(ldh*3)+i], hh[(ldh*3)+i]);
1129 : #endif
1130 : #endif
1131 :
1132 229386240 : q1 = _SSE_SUB(q1, _SSE_MUL(w1,h4));
1133 229386240 : q2 = _SSE_SUB(q2, _SSE_MUL(w2,h4));
1134 229386240 : q3 = _SSE_SUB(q3, _SSE_MUL(w3,h4));
1135 :
1136 114693120 : _SSE_STORE(&q[i*ldq],q1);
1137 114693120 : _SSE_STORE(&q[(i*ldq)+offset],q2);
1138 114693120 : _SSE_STORE(&q[(i*ldq)+2*offset],q3);
1139 : }
1140 :
1141 : #ifdef HAVE_SSE_INTRINSICS
1142 : #ifdef DOUBLE_PRECISION_REAL
1143 3219456 : h1 = _mm_set1_pd(hh[nb-3]);
1144 : #endif
1145 : #ifdef SINGLE_PRECISION_REAL
1146 603648 : h1 = _mm_set1_ps(hh[nb-3]);
1147 : #endif
1148 : #endif
1149 :
1150 : #ifdef HAVE_SPARC64_SSE
1151 : #ifdef DOUBLE_PRECISION_REAL
1152 : h1 = _mm_set_pd(hh[nb-3], hh[nb-3]);
1153 : #endif
1154 : #ifdef SINGLE_PRECISION_REAL
1155 : h1 = _mm_set_ps(hh[nb-3], hh[nb-3]);
1156 : #endif
1157 : #endif
1158 :
1159 3823104 : q1 = _SSE_LOAD(&q[nb*ldq]);
1160 3823104 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
1161 3823104 : q3 = _SSE_LOAD(&q[(nb*ldq)+2*offset]);
1162 :
1163 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(x1, h1));
1164 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(x2, h1));
1165 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(x3, h1));
1166 :
1167 : #ifdef HAVE_SSE_INTRINSICS
1168 : #ifdef DOUBLE_PRECISION_REAL
1169 3219456 : h2 = _mm_set1_pd(hh[ldh+nb-2]);
1170 : #endif
1171 : #ifdef SINGLE_PRECISION_REAL
1172 603648 : h2 = _mm_set1_ps(hh[ldh+nb-2]);
1173 : #endif
1174 : #endif
1175 :
1176 : #ifdef HAVE_SPARC64_SSE
1177 : #ifdef DOUBLE_PRECISION_REAL
1178 : h2 = _mm_set_pd(hh[ldh+nb-2], hh[ldh+nb-2]);
1179 : #endif
1180 : #ifdef SINGLE_PRECISION_REAL
1181 : h2 = _mm_set_ps(hh[ldh+nb-2], hh[ldh+nb-2]);
1182 : #endif
1183 : #endif
1184 :
1185 :
1186 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(y1, h2));
1187 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(y2, h2));
1188 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(y3, h2));
1189 :
1190 : #ifdef HAVE_SSE_INTRINSICS
1191 : #ifdef DOUBLE_PRECISION_REAL
1192 3219456 : h3 = _mm_set1_pd(hh[(ldh*2)+nb-1]);
1193 : #endif
1194 : #ifdef SINGLE_PRECISION_REAL
1195 603648 : h3 = _mm_set1_ps(hh[(ldh*2)+nb-1]);
1196 : #endif
1197 : #endif
1198 :
1199 : #ifdef HAVE_SPARC64_SSE
1200 : #ifdef DOUBLE_PRECISION_REAL
1201 : h3 = _mm_set_pd(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
1202 : #endif
1203 : #ifdef SINGLE_PRECISION_REAL
1204 : h3 = _mm_set_ps(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
1205 : #endif
1206 : #endif
1207 :
1208 :
1209 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(z1, h3));
1210 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(z2, h3));
1211 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(z3, h3));
1212 :
1213 1911552 : _SSE_STORE(&q[nb*ldq],q1);
1214 1911552 : _SSE_STORE(&q[(nb*ldq)+offset],q2);
1215 1911552 : _SSE_STORE(&q[(nb*ldq)+2*offset],q3);
1216 :
1217 : #ifdef HAVE_SSE_INTRINSICS
1218 : #ifdef DOUBLE_PRECISION_REAL
1219 3219456 : h1 = _mm_set1_pd(hh[nb-2]);
1220 : #endif
1221 : #ifdef SINGLE_PRECISION_REAL
1222 603648 : h1 = _mm_set1_ps(hh[nb-2]);
1223 : #endif
1224 : #endif
1225 :
1226 : #ifdef HAVE_SPARC64_SSE
1227 : #ifdef DOUBLE_PRECISION_REAL
1228 : h1 = _mm_set_pd(hh[nb-2], hh[nb-2]);
1229 : #endif
1230 : #ifdef SINGLE_PRECISION_REAL
1231 : h1 = _mm_set_ps(hh[nb-2], hh[nb-2]);
1232 : #endif
1233 : #endif
1234 :
1235 3823104 : q1 = _SSE_LOAD(&q[(nb+1)*ldq]);
1236 3823104 : q2 = _SSE_LOAD(&q[((nb+1)*ldq)+offset]);
1237 3823104 : q3 = _SSE_LOAD(&q[((nb+1)*ldq)+2*offset]);
1238 :
1239 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(x1, h1));
1240 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(x2, h1));
1241 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(x3, h1));
1242 :
1243 : #ifdef HAVE_SSE_INTRINSICS
1244 : #ifdef DOUBLE_PRECISION_REAL
1245 3219456 : h2 = _mm_set1_pd(hh[ldh+nb-1]);
1246 : #endif
1247 : #ifdef SINGLE_PRECISION_REAL
1248 603648 : h2 = _mm_set1_ps(hh[ldh+nb-1]);
1249 : #endif
1250 : #endif
1251 :
1252 : #ifdef HAVE_SPARC64_SSE
1253 : #ifdef DOUBLE_PRECISION_REAL
1254 : h2 = _mm_set_pd(hh[ldh+nb-1], hh[ldh+nb-1]);
1255 : #endif
1256 : #ifdef SINGLE_PRECISION_REAL
1257 : h2 = _mm_set_ps(hh[ldh+nb-1], hh[ldh+nb-1]);
1258 : #endif
1259 : #endif
1260 :
1261 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(y1, h2));
1262 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(y2, h2));
1263 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(y3, h2));
1264 :
1265 1911552 : _SSE_STORE(&q[(nb+1)*ldq],q1);
1266 1911552 : _SSE_STORE(&q[((nb+1)*ldq)+offset],q2);
1267 1911552 : _SSE_STORE(&q[((nb+1)*ldq)+2*offset],q3);
1268 :
1269 : #ifdef HAVE_SSE_INTRINSICS
1270 : #ifdef DOUBLE_PRECISION_REAL
1271 3219456 : h1 = _mm_set1_pd(hh[nb-1]);
1272 : #endif
1273 : #ifdef SINGLE_PRECISION_REAL
1274 603648 : h1 = _mm_set1_ps(hh[nb-1]);
1275 : #endif
1276 : #endif
1277 :
1278 : #ifdef HAVE_SPARC64_SSE
1279 : #ifdef DOUBLE_PRECISION_REAL
1280 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
1281 : #endif
1282 : #ifdef SINGLE_PRECISION_REAL
1283 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
1284 : #endif
1285 : #endif
1286 :
1287 3823104 : q1 = _SSE_LOAD(&q[(nb+2)*ldq]);
1288 3823104 : q2 = _SSE_LOAD(&q[((nb+2)*ldq)+offset]);
1289 3823104 : q3 = _SSE_LOAD(&q[((nb+2)*ldq)+2*offset]);
1290 :
1291 3823104 : q1 = _SSE_SUB(q1, _SSE_MUL(x1, h1));
1292 3823104 : q2 = _SSE_SUB(q2, _SSE_MUL(x2, h1));
1293 3823104 : q3 = _SSE_SUB(q3, _SSE_MUL(x3, h1));
1294 :
1295 1911552 : _SSE_STORE(&q[(nb+2)*ldq],q1);
1296 1911552 : _SSE_STORE(&q[((nb+2)*ldq)+offset],q2);
1297 1911552 : _SSE_STORE(&q[((nb+2)*ldq)+2*offset],q3);
1298 : }
1299 :
1300 : /**
1301 : * Unrolled kernel that computes
1302 : #ifdef DOUBLE_PRECISION_REAL
1303 : * 4 rows of Q simultaneously, a
1304 : #endif
1305 : #ifdef SINGLE_PRECISION_REAL
1306 : * 8 rows of Q simultaneously, a
1307 : #endif
1308 : * matrix Vector product with two householder
1309 : * vectors + a rank 1 update is performed
1310 : */
1311 : #ifdef HAVE_SSE_INTRINSICS
1312 : #ifdef DOUBLE_PRECISION_REAL
1313 : __forceinline void hh_trafo_kernel_4_SSE_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4)
1314 : #endif
1315 : #ifdef SINGLE_PRECISION_REAL
1316 : __forceinline void hh_trafo_kernel_8_SSE_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4)
1317 : #endif
1318 : #endif
1319 : #ifdef HAVE_SPARC64_SSE
1320 : #ifdef DOUBLE_PRECISION_REAL
1321 : __forceinline void hh_trafo_kernel_4_SPARC64_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4)
1322 : #endif
1323 : #ifdef SINGLE_PRECISION_REAL
1324 : __forceinline void hh_trafo_kernel_8_SPARC64_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4)
1325 : #endif
1326 : #endif
1327 :
1328 : {
1329 : /////////////////////////////////////////////////////
1330 : // Matrix Vector Multiplication, Q [4 x nb+3] * hh
1331 : // hh contains four householder vectors
1332 : /////////////////////////////////////////////////////
1333 : int i;
1334 :
1335 637184 : __SSE_DATATYPE a1_1 = _SSE_LOAD(&q[ldq*3]);
1336 637184 : __SSE_DATATYPE a2_1 = _SSE_LOAD(&q[ldq*2]);
1337 637184 : __SSE_DATATYPE a3_1 = _SSE_LOAD(&q[ldq]);
1338 318592 : __SSE_DATATYPE a4_1 = _SSE_LOAD(&q[0]);
1339 :
1340 : #ifdef HAVE_SSE_INTRINSICS
1341 : #ifdef DOUBLE_PRECISION_REAL
1342 536576 : __SSE_DATATYPE h_2_1 = _mm_set1_pd(hh[ldh+1]);
1343 536576 : __SSE_DATATYPE h_3_2 = _mm_set1_pd(hh[(ldh*2)+1]);
1344 536576 : __SSE_DATATYPE h_3_1 = _mm_set1_pd(hh[(ldh*2)+2]);
1345 536576 : __SSE_DATATYPE h_4_3 = _mm_set1_pd(hh[(ldh*3)+1]);
1346 536576 : __SSE_DATATYPE h_4_2 = _mm_set1_pd(hh[(ldh*3)+2]);
1347 536576 : __SSE_DATATYPE h_4_1 = _mm_set1_pd(hh[(ldh*3)+3]);
1348 : #endif
1349 :
1350 : #ifdef SINGLE_PRECISION_REAL
1351 100608 : __m128 h_2_1 = _mm_set1_ps(hh[ldh+1] ); // h_2_1 contains four times hh[ldh+1]
1352 100608 : __m128 h_3_2 = _mm_set1_ps(hh[(ldh*2)+1]);
1353 100608 : __m128 h_3_1 = _mm_set1_ps(hh[(ldh*2)+2]);
1354 100608 : __m128 h_4_3 = _mm_set1_ps(hh[(ldh*3)+1]);
1355 100608 : __m128 h_4_2 = _mm_set1_ps(hh[(ldh*3)+2]);
1356 100608 : __m128 h_4_1 = _mm_set1_ps(hh[(ldh*3)+3]);
1357 : #endif
1358 : #endif
1359 :
1360 : #ifdef HAVE_SPARC64_SSE
1361 : #ifdef DOUBLE_PRECISION_REAL
1362 : __SSE_DATATYPE h_2_1 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1363 : __SSE_DATATYPE h_3_2 = _mm_set_pd(hh[(ldh*2)+1], hh[(ldh*2)+1]);
1364 : __SSE_DATATYPE h_3_1 = _mm_set_pd(hh[(ldh*2)+2], hh[(ldh*2)+2]);
1365 : __SSE_DATATYPE h_4_3 = _mm_set_pd(hh[(ldh*3)+1], hh[(ldh*3)+1]);
1366 : __SSE_DATATYPE h_4_2 = _mm_set_pd(hh[(ldh*3)+2], hh[(ldh*3)+2]);
1367 : __SSE_DATATYPE h_4_1 = _mm_set_pd(hh[(ldh*3)+3], hh[(ldh*3)+3]);
1368 : #endif
1369 :
1370 : #ifdef SINGLE_PRECISION_REAL
1371 : __m128 h_2_1 = _mm_set_ps(hh[ldh+1], hh[ldh+1]); // h_2_1 contains four times hh[ldh+1]
1372 : __m128 h_3_2 = _mm_set_ps(hh[(ldh*2)+1], hh[(ldh*2)+1]);
1373 : __m128 h_3_1 = _mm_set_ps(hh[(ldh*2)+2], hh[(ldh*2)+2]);
1374 : __m128 h_4_3 = _mm_set_ps(hh[(ldh*3)+1], hh[(ldh*3)+1]);
1375 : __m128 h_4_2 = _mm_set_ps(hh[(ldh*3)+2], hh[(ldh*3)+2]);
1376 : __m128 h_4_1 = _mm_set_ps(hh[(ldh*3)+3], hh[(ldh*3)+3]);
1377 : #endif
1378 : #endif
1379 :
1380 :
1381 637184 : __SSE_DATATYPE w1 = _SSE_ADD(a4_1, _SSE_MUL(a3_1, h_4_3));
1382 637184 : w1 = _SSE_ADD(w1, _SSE_MUL(a2_1, h_4_2));
1383 637184 : w1 = _SSE_ADD(w1, _SSE_MUL(a1_1, h_4_1));
1384 637184 : __SSE_DATATYPE z1 = _SSE_ADD(a3_1, _SSE_MUL(a2_1, h_3_2));
1385 637184 : z1 = _SSE_ADD(z1, _SSE_MUL(a1_1, h_3_1));
1386 637184 : __SSE_DATATYPE y1 = _SSE_ADD(a2_1, _SSE_MUL(a1_1, h_2_1));
1387 318592 : __SSE_DATATYPE x1 = a1_1;
1388 :
1389 637184 : __SSE_DATATYPE a1_2 = _SSE_LOAD(&q[(ldq*3)+offset]);
1390 637184 : __SSE_DATATYPE a2_2 = _SSE_LOAD(&q[(ldq*2)+offset]);
1391 637184 : __SSE_DATATYPE a3_2 = _SSE_LOAD(&q[ldq+offset]);
1392 637184 : __SSE_DATATYPE a4_2 = _SSE_LOAD(&q[0+offset]);
1393 :
1394 637184 : __SSE_DATATYPE w2 = _SSE_ADD(a4_2, _SSE_MUL(a3_2, h_4_3));
1395 637184 : w2 = _SSE_ADD(w2, _SSE_MUL(a2_2, h_4_2));
1396 637184 : w2 = _SSE_ADD(w2, _SSE_MUL(a1_2, h_4_1));
1397 637184 : __SSE_DATATYPE z2 = _SSE_ADD(a3_2, _SSE_MUL(a2_2, h_3_2));
1398 637184 : z2 = _SSE_ADD(z2, _SSE_MUL(a1_2, h_3_1));
1399 637184 : __SSE_DATATYPE y2 = _SSE_ADD(a2_2, _SSE_MUL(a1_2, h_2_1));
1400 318592 : __SSE_DATATYPE x2 = a1_2;
1401 :
1402 : __SSE_DATATYPE q1;
1403 : __SSE_DATATYPE q2;
1404 :
1405 : __SSE_DATATYPE h1;
1406 : __SSE_DATATYPE h2;
1407 : __SSE_DATATYPE h3;
1408 : __SSE_DATATYPE h4;
1409 :
1410 19434112 : for(i = 4; i < nb; i++)
1411 : {
1412 :
1413 : #ifdef HAVE_SSE_INTRINSICS
1414 : #ifdef DOUBLE_PRECISION_REAL
1415 32194560 : h1 = _mm_set1_pd(hh[i-3]);
1416 32194560 : h2 = _mm_set1_pd(hh[ldh+i-2]);
1417 32194560 : h3 = _mm_set1_pd(hh[(ldh*2)+i-1]);
1418 32194560 : h4 = _mm_set1_pd(hh[(ldh*3)+i]);
1419 : #endif
1420 :
1421 : #ifdef SINGLE_PRECISION_REAL
1422 6036480 : h1 = _mm_set1_ps(hh[i-3]);
1423 6036480 : h2 = _mm_set1_ps(hh[ldh+i-2]);
1424 6036480 : h3 = _mm_set1_ps(hh[(ldh*2)+i-1]);
1425 6036480 : h4 = _mm_set1_ps(hh[(ldh*3)+i]);
1426 : #endif
1427 : #endif
1428 :
1429 : #ifdef HAVE_SPARC64_SSE
1430 : #ifdef DOUBLE_PRECISION_REAL
1431 : h1 = _mm_set_pd(hh[i-3], hh[i-3]);
1432 : h2 = _mm_set_pd(hh[ldh+i-2], hh[ldh+i-2]);
1433 : h3 = _mm_set_pd(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
1434 : h4 = _mm_set_pd(hh[(ldh*3)+i], hh[(ldh*3)+i]);
1435 : #endif
1436 :
1437 : #ifdef SINGLE_PRECISION_REAL
1438 : h1 = _mm_set_ps(hh[i-3], hh[i-3]);
1439 : h2 = _mm_set_ps(hh[ldh+i-2], hh[ldh+i-2]);
1440 : h3 = _mm_set_ps(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
1441 : h4 = _mm_set_ps(hh[(ldh*3)+i], hh[(ldh*3)+i]);
1442 : #endif
1443 : #endif
1444 :
1445 38231040 : q1 = _SSE_LOAD(&q[i*ldq]);
1446 :
1447 38231040 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
1448 38231040 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
1449 38231040 : z1 = _SSE_ADD(z1, _SSE_MUL(q1,h3));
1450 38231040 : w1 = _SSE_ADD(w1, _SSE_MUL(q1,h4));
1451 :
1452 38231040 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
1453 :
1454 38231040 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
1455 38231040 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
1456 38231040 : z2 = _SSE_ADD(z2, _SSE_MUL(q2,h3));
1457 38231040 : w2 = _SSE_ADD(w2, _SSE_MUL(q2,h4));
1458 : }
1459 :
1460 : #ifdef HAVE_SSE_INTRINSICS
1461 : #ifdef DOUBLE_PRECISION_REAL
1462 536576 : h1 = _mm_set1_pd(hh[nb-3]);
1463 536576 : h2 = _mm_set1_pd(hh[ldh+nb-2]);
1464 536576 : h3 = _mm_set1_pd(hh[(ldh*2)+nb-1]);
1465 : #endif
1466 : #ifdef SINGLE_PRECISION_REAL
1467 100608 : h1 = _mm_set1_ps(hh[nb-3]);
1468 100608 : h2 = _mm_set1_ps(hh[ldh+nb-2]);
1469 100608 : h3 = _mm_set1_ps(hh[(ldh*2)+nb-1]);
1470 : #endif
1471 : #endif
1472 :
1473 : #ifdef HAVE_SPARC64_SSE
1474 : #ifdef DOUBLE_PRECISION_REAL
1475 : h1 = _mm_set_pd(hh[nb-3], hh[nb-3]);
1476 : h2 = _mm_set_pd(hh[ldh+nb-2], hh[ldh+nb-2]);
1477 : h3 = _mm_set_pd(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
1478 : #endif
1479 : #ifdef SINGLE_PRECISION_REAL
1480 : h1 = _mm_set_ps(hh[nb-3], hh[nb-3]);
1481 : h2 = _mm_set_ps(hh[ldh+nb-2], hh[ldh+nb-2]);
1482 : h3 = _mm_set_ps(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
1483 : #endif
1484 : #endif
1485 :
1486 637184 : q1 = _SSE_LOAD(&q[nb*ldq]);
1487 637184 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
1488 :
1489 637184 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
1490 637184 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
1491 637184 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
1492 637184 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
1493 637184 : z1 = _SSE_ADD(z1, _SSE_MUL(q1,h3));
1494 637184 : z2 = _SSE_ADD(z2, _SSE_MUL(q2,h3));
1495 :
1496 : #ifdef HAVE_SSE_INTRINSICS
1497 : #ifdef DOUBLE_PRECISION_REAL
1498 536576 : h1 = _mm_set1_pd(hh[nb-2]);
1499 536576 : h2 = _mm_set1_pd(hh[(ldh*1)+nb-1]);
1500 : #endif
1501 : #ifdef SINGLE_PRECISION_REAL
1502 100608 : h1 = _mm_set1_ps(hh[nb-2]);
1503 100608 : h2 = _mm_set1_ps(hh[(ldh*1)+nb-1]);
1504 : #endif
1505 : #endif
1506 :
1507 : #ifdef HAVE_SPARC64_SSE
1508 : #ifdef DOUBLE_PRECISION_REAL
1509 : h1 = _mm_set_pd(hh[nb-2], hh[nb-2]);
1510 : h2 = _mm_set_pd(hh[(ldh*1)+nb-1], hh[(ldh*1)+nb-1]);
1511 : #endif
1512 : #ifdef SINGLE_PRECISION_REAL
1513 : h1 = _mm_set_ps(hh[nb-2], hh[nb-2]);
1514 : h2 = _mm_set_ps(hh[(ldh*1)+nb-1], hh[(ldh*1)+nb-1]);
1515 : #endif
1516 : #endif
1517 :
1518 :
1519 637184 : q1 = _SSE_LOAD(&q[(nb+1)*ldq]);
1520 637184 : q2 = _SSE_LOAD(&q[((nb+1)*ldq)+offset]);
1521 :
1522 637184 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
1523 637184 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
1524 637184 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
1525 637184 : y2 = _SSE_ADD(y2, _SSE_MUL(q2,h2));
1526 :
1527 : #ifdef HAVE_SSE_INTRINSICS
1528 : #ifdef DOUBLE_PRECISION_REAL
1529 536576 : h1 = _mm_set1_pd(hh[nb-1]);
1530 : #endif
1531 : #ifdef SINGLE_PRECISION_REAL
1532 100608 : h1 = _mm_set1_ps(hh[nb-1]);
1533 : #endif
1534 : #endif
1535 :
1536 : #ifdef HAVE_SPARC64_SSE
1537 : #ifdef DOUBLE_PRECISION_REAL
1538 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
1539 : #endif
1540 : #ifdef SINGLE_PRECISION_REAL
1541 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
1542 : #endif
1543 : #endif
1544 :
1545 :
1546 637184 : q1 = _SSE_LOAD(&q[(nb+2)*ldq]);
1547 637184 : q2 = _SSE_LOAD(&q[((nb+2)*ldq)+offset]);
1548 :
1549 637184 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
1550 637184 : x2 = _SSE_ADD(x2, _SSE_MUL(q2,h1));
1551 :
1552 : /////////////////////////////////////////////////////
1553 : // Rank-1 update of Q [4 x nb+3]
1554 : /////////////////////////////////////////////////////
1555 :
1556 : #ifdef HAVE_SSE_INTRINSICS
1557 : #ifdef DOUBLE_PRECISION_REAL
1558 536576 : __SSE_DATATYPE tau1 = _mm_set1_pd(hh[0]);
1559 536576 : __SSE_DATATYPE tau2 = _mm_set1_pd(hh[ldh]);
1560 536576 : __SSE_DATATYPE tau3 = _mm_set1_pd(hh[ldh*2]);
1561 536576 : __SSE_DATATYPE tau4 = _mm_set1_pd(hh[ldh*3]);
1562 :
1563 268288 : __SSE_DATATYPE vs_1_2 = _mm_set1_pd(s_1_2);
1564 268288 : __SSE_DATATYPE vs_1_3 = _mm_set1_pd(s_1_3);
1565 268288 : __SSE_DATATYPE vs_2_3 = _mm_set1_pd(s_2_3);
1566 268288 : __SSE_DATATYPE vs_1_4 = _mm_set1_pd(s_1_4);
1567 268288 : __SSE_DATATYPE vs_2_4 = _mm_set1_pd(s_2_4);
1568 268288 : __SSE_DATATYPE vs_3_4 = _mm_set1_pd(s_3_4);
1569 : #endif
1570 :
1571 : #ifdef SINGLE_PRECISION_REAL
1572 100608 : __m128 tau1 = _mm_set1_ps(hh[0]);
1573 100608 : __m128 tau2 = _mm_set1_ps(hh[ldh]);
1574 100608 : __m128 tau3 = _mm_set1_ps(hh[ldh*2]);
1575 100608 : __m128 tau4 = _mm_set1_ps(hh[ldh*3]);
1576 :
1577 50304 : __m128 vs_1_2 = _mm_set1_ps(s_1_2);
1578 50304 : __m128 vs_1_3 = _mm_set1_ps(s_1_3);
1579 50304 : __m128 vs_2_3 = _mm_set1_ps(s_2_3);
1580 50304 : __m128 vs_1_4 = _mm_set1_ps(s_1_4);
1581 50304 : __m128 vs_2_4 = _mm_set1_ps(s_2_4);
1582 50304 : __m128 vs_3_4 = _mm_set1_ps(s_3_4);
1583 : #endif
1584 : #endif
1585 :
1586 : #ifdef HAVE_SPARC64_SSE
1587 : #ifdef DOUBLE_PRECISION_REAL
1588 : __SSE_DATATYPE tau1 = _mm_set_pd(hh[0], hh[0]);
1589 : __SSE_DATATYPE tau2 = _mm_set_pd(hh[ldh], hh[ldh]);
1590 : __SSE_DATATYPE tau3 = _mm_set_pd(hh[ldh*2], hh[ldh*2]);
1591 : __SSE_DATATYPE tau4 = _mm_set_pd(hh[ldh*3], hh[ldh*3]);
1592 :
1593 : __SSE_DATATYPE vs_1_2 = _mm_set_pd(s_1_2, s_1_2);
1594 : __SSE_DATATYPE vs_1_3 = _mm_set_pd(s_1_3, s_1_3);
1595 : __SSE_DATATYPE vs_2_3 = _mm_set_pd(s_2_3, s_2_3);
1596 : __SSE_DATATYPE vs_1_4 = _mm_set_pd(s_1_4, s_1_4);
1597 : __SSE_DATATYPE vs_2_4 = _mm_set_pd(s_2_4, s_2_4);
1598 : __SSE_DATATYPE vs_3_4 = _mm_set_pd(s_3_4, s_3_4);
1599 : #endif
1600 :
1601 : #ifdef SINGLE_PRECISION_REAL
1602 : __m128 tau1 = _mm_set_ps(hh[0], hh[0]);
1603 : __m128 tau2 = _mm_set_ps(hh[ldh], hh[ldh]);
1604 : __m128 tau3 = _mm_set_ps(hh[ldh*2], hh[ldh*2]);
1605 : __m128 tau4 = _mm_set_ps(hh[ldh*3], hh[ldh*3]);
1606 :
1607 : __m128 vs_1_2 = _mm_set_ps(s_1_2, s_1_2);
1608 : __m128 vs_1_3 = _mm_set_ps(s_1_3, s_1_3);
1609 : __m128 vs_2_3 = _mm_set_ps(s_2_3, s_2_3);
1610 : __m128 vs_1_4 = _mm_set_ps(s_1_4, s_1_4);
1611 : __m128 vs_2_4 = _mm_set_ps(s_2_4, s_2_4);
1612 : __m128 vs_3_4 = _mm_set_ps(s_3_4, s_3_4);
1613 : #endif
1614 : #endif
1615 :
1616 :
1617 318592 : h1 = tau1;
1618 318592 : x1 = _SSE_MUL(x1, h1);
1619 318592 : x2 = _SSE_MUL(x2, h1);
1620 :
1621 318592 : h1 = tau2;
1622 318592 : h2 = _SSE_MUL(h1, vs_1_2);
1623 :
1624 955776 : y1 = _SSE_SUB(_SSE_MUL(y1,h1), _SSE_MUL(x1,h2));
1625 955776 : y2 = _SSE_SUB(_SSE_MUL(y2,h1), _SSE_MUL(x2,h2));
1626 :
1627 318592 : h1 = tau3;
1628 318592 : h2 = _SSE_MUL(h1, vs_1_3);
1629 318592 : h3 = _SSE_MUL(h1, vs_2_3);
1630 :
1631 1592960 : z1 = _SSE_SUB(_SSE_MUL(z1,h1), _SSE_ADD(_SSE_MUL(y1,h3), _SSE_MUL(x1,h2)));
1632 1592960 : z2 = _SSE_SUB(_SSE_MUL(z2,h1), _SSE_ADD(_SSE_MUL(y2,h3), _SSE_MUL(x2,h2)));
1633 :
1634 318592 : h1 = tau4;
1635 318592 : h2 = _SSE_MUL(h1, vs_1_4);
1636 318592 : h3 = _SSE_MUL(h1, vs_2_4);
1637 318592 : h4 = _SSE_MUL(h1, vs_3_4);
1638 :
1639 2230144 : w1 = _SSE_SUB(_SSE_MUL(w1,h1), _SSE_ADD(_SSE_MUL(z1,h4), _SSE_ADD(_SSE_MUL(y1,h3), _SSE_MUL(x1,h2))));
1640 2230144 : w2 = _SSE_SUB(_SSE_MUL(w2,h1), _SSE_ADD(_SSE_MUL(z2,h4), _SSE_ADD(_SSE_MUL(y2,h3), _SSE_MUL(x2,h2))));
1641 :
1642 318592 : q1 = _SSE_LOAD(&q[0]);
1643 637184 : q2 = _SSE_LOAD(&q[offset]);
1644 318592 : q1 = _SSE_SUB(q1, w1);
1645 318592 : q2 = _SSE_SUB(q2, w2);
1646 : _SSE_STORE(&q[0],q1);
1647 318592 : _SSE_STORE(&q[offset],q2);
1648 :
1649 : #ifdef HAVE_SSE_INTRINSICS
1650 : #ifdef DOUBLE_PRECISION_REAL
1651 536576 : h4 = _mm_set1_pd(hh[(ldh*3)+1]);
1652 : #endif
1653 : #ifdef SINGLE_PRECISION_REAL
1654 100608 : h4 = _mm_set1_ps(hh[(ldh*3)+1]);
1655 : #endif
1656 : #endif
1657 :
1658 : #ifdef HAVE_SPARC64_SSE
1659 : #ifdef DOUBLE_PRECISION_REAL
1660 : h4 = _mm_set_pd(hh[(ldh*3)+1], hh[(ldh*3)+1]);
1661 : #endif
1662 : #ifdef SINGLE_PRECISION_REAL
1663 : h4 = _mm_set_ps(hh[(ldh*3)+1], hh[(ldh*3)+1]);
1664 : #endif
1665 : #endif
1666 :
1667 637184 : q1 = _SSE_LOAD(&q[ldq]);
1668 637184 : q2 = _SSE_LOAD(&q[ldq+offset]);
1669 :
1670 955776 : q1 = _SSE_SUB(q1, _SSE_ADD(z1, _SSE_MUL(w1, h4)));
1671 955776 : q2 = _SSE_SUB(q2, _SSE_ADD(z2, _SSE_MUL(w2, h4)));
1672 :
1673 318592 : _SSE_STORE(&q[ldq],q1);
1674 318592 : _SSE_STORE(&q[ldq+offset],q2);
1675 :
1676 : #ifdef HAVE_SSE_INTRINSICS
1677 : #ifdef DOUBLE_PRECISION_REAL
1678 536576 : h3 = _mm_set1_pd(hh[(ldh*2)+1]);
1679 536576 : h4 = _mm_set1_pd(hh[(ldh*3)+2]);
1680 : #endif
1681 : #ifdef SINGLE_PRECISION_REAL
1682 100608 : h3 = _mm_set1_ps(hh[(ldh*2)+1]);
1683 100608 : h4 = _mm_set1_ps(hh[(ldh*3)+2]);
1684 : #endif
1685 : #endif
1686 :
1687 : #ifdef HAVE_SPARC64_SSE
1688 : #ifdef DOUBLE_PRECISION_REAL
1689 : h3 = _mm_set_pd(hh[(ldh*2)+1], hh[(ldh*2)+1]);
1690 : h4 = _mm_set_pd(hh[(ldh*3)+2], hh[(ldh*3)+2]);
1691 : #endif
1692 : #ifdef SINGLE_PRECISION_REAL
1693 : h3 = _mm_set_ps(hh[(ldh*2)+1], hh[(ldh*2)+1]);
1694 : h4 = _mm_set_ps(hh[(ldh*3)+2], hh[(ldh*3)+2]);
1695 : #endif
1696 : #endif
1697 :
1698 :
1699 637184 : q1 = _SSE_LOAD(&q[ldq*2]);
1700 637184 : q2 = _SSE_LOAD(&q[(ldq*2)+offset]);
1701 :
1702 1592960 : q1 = _SSE_SUB(q1, _SSE_ADD(y1, _SSE_ADD(_SSE_MUL(z1, h3), _SSE_MUL(w1, h4))));
1703 1592960 : q2 = _SSE_SUB(q2, _SSE_ADD(y2, _SSE_ADD(_SSE_MUL(z2, h3), _SSE_MUL(w2, h4))));
1704 318592 : _SSE_STORE(&q[ldq*2],q1);
1705 318592 : _SSE_STORE(&q[(ldq*2)+offset],q2);
1706 :
1707 : #ifdef HAVE_SSE_INTRINSICS
1708 : #ifdef DOUBLE_PRECISION_REAL
1709 536576 : h2 = _mm_set1_pd(hh[ldh+1]);
1710 536576 : h3 = _mm_set1_pd(hh[(ldh*2)+2]);
1711 536576 : h4 = _mm_set1_pd(hh[(ldh*3)+3]);
1712 : #endif
1713 : #ifdef SINGLE_PRECISION_REAL
1714 100608 : h2 = _mm_set1_ps(hh[ldh+1]);
1715 100608 : h3 = _mm_set1_ps(hh[(ldh*2)+2]);
1716 100608 : h4 = _mm_set1_ps(hh[(ldh*3)+3]);
1717 :
1718 : #endif
1719 : #endif
1720 :
1721 : #ifdef HAVE_SPARC64_SSE
1722 : #ifdef DOUBLE_PRECISION_REAL
1723 : h2 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1724 : h3 = _mm_set_pd(hh[(ldh*2)+2], hh[(ldh*2)+2]);
1725 : h4 = _mm_set_pd(hh[(ldh*3)+3], hh[(ldh*3)+3]);
1726 : #endif
1727 : #ifdef SINGLE_PRECISION_REAL
1728 : h2 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
1729 : h3 = _mm_set_ps(hh[(ldh*2)+2], hh[(ldh*2)+2]);
1730 : h4 = _mm_set_ps(hh[(ldh*3)+3], hh[(ldh*3)+3]);
1731 :
1732 : #endif
1733 : #endif
1734 :
1735 637184 : q1 = _SSE_LOAD(&q[ldq*3]);
1736 637184 : q2 = _SSE_LOAD(&q[(ldq*3)+offset]);
1737 :
1738 2230144 : q1 = _SSE_SUB(q1, _SSE_ADD(x1, _SSE_ADD(_SSE_MUL(y1, h2), _SSE_ADD(_SSE_MUL(z1, h3), _SSE_MUL(w1, h4)))));
1739 2230144 : q2 = _SSE_SUB(q2, _SSE_ADD(x2, _SSE_ADD(_SSE_MUL(y2, h2), _SSE_ADD(_SSE_MUL(z2, h3), _SSE_MUL(w2, h4)))));
1740 :
1741 318592 : _SSE_STORE(&q[ldq*3], q1);
1742 318592 : _SSE_STORE(&q[(ldq*3)+offset], q2);
1743 :
1744 19434112 : for (i = 4; i < nb; i++)
1745 : {
1746 : #ifdef HAVE_SSE_INTRINSICS
1747 : #ifdef DOUBLE_PRECISION_REAL
1748 32194560 : h1 = _mm_set1_pd(hh[i-3]);
1749 32194560 : h2 = _mm_set1_pd(hh[ldh+i-2]);
1750 32194560 : h3 = _mm_set1_pd(hh[(ldh*2)+i-1]);
1751 32194560 : h4 = _mm_set1_pd(hh[(ldh*3)+i]);
1752 : #endif
1753 : #ifdef SINGLE_PRECISION_REAL
1754 6036480 : h1 = _mm_set1_ps(hh[i-3]);
1755 6036480 : h2 = _mm_set1_ps(hh[ldh+i-2]);
1756 6036480 : h3 = _mm_set1_ps(hh[(ldh*2)+i-1]);
1757 6036480 : h4 = _mm_set1_ps(hh[(ldh*3)+i]);
1758 : #endif
1759 : #endif
1760 :
1761 : #ifdef HAVE_SPARC64_SSE
1762 : #ifdef DOUBLE_PRECISION_REAL
1763 : h1 = _mm_set_pd(hh[i-3], hh[i-3]);
1764 : h2 = _mm_set_pd(hh[ldh+i-2], hh[ldh+i-2]);
1765 : h3 = _mm_set_pd(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
1766 : h4 = _mm_set_pd(hh[(ldh*3)+i], hh[(ldh*3)+i]);
1767 : #endif
1768 : #ifdef SINGLE_PRECISION_REAL
1769 : h1 = _mm_set_ps(hh[i-3], hh[i-3]);
1770 : h2 = _mm_set_ps(hh[ldh+i-2], hh[ldh+i-2]);
1771 : h3 = _mm_set_ps(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
1772 : h4 = _mm_set_ps(hh[(ldh*3)+i], hh[(ldh*3)+i]);
1773 : #endif
1774 : #endif
1775 :
1776 :
1777 38231040 : q1 = _SSE_LOAD(&q[i*ldq]);
1778 :
1779 152924160 : q1 = _SSE_SUB(q1, _SSE_ADD(_SSE_ADD(_SSE_MUL(w1, h4), _SSE_MUL(z1, h3)), _SSE_ADD(_SSE_MUL(x1,h1), _SSE_MUL(y1, h2))));
1780 :
1781 19115520 : _SSE_STORE(&q[i*ldq],q1);
1782 :
1783 38231040 : q2 = _SSE_LOAD(&q[(i*ldq)+offset]);
1784 :
1785 152924160 : q2 = _SSE_SUB(q2, _SSE_ADD(_SSE_ADD(_SSE_MUL(w2, h4), _SSE_MUL(z2, h3)), _SSE_ADD(_SSE_MUL(x2,h1), _SSE_MUL(y2, h2))));
1786 :
1787 19115520 : _SSE_STORE(&q[(i*ldq)+offset],q2);
1788 : }
1789 :
1790 : #ifdef HAVE_SSE_INTRINSICS
1791 : #ifdef DOUBLE_PRECISION_REAL
1792 536576 : h1 = _mm_set1_pd(hh[nb-3]);
1793 536576 : h2 = _mm_set1_pd(hh[ldh+nb-2]);
1794 536576 : h3 = _mm_set1_pd(hh[(ldh*2)+nb-1]);
1795 : #endif
1796 : #ifdef SINGLE_PRECISION_REAL
1797 100608 : h1 = _mm_set1_ps(hh[nb-3]);
1798 100608 : h2 = _mm_set1_ps(hh[ldh+nb-2]);
1799 100608 : h3 = _mm_set1_ps(hh[(ldh*2)+nb-1]);
1800 : #endif
1801 : #endif
1802 :
1803 : #ifdef HAVE_SPARC64_SSE
1804 : #ifdef DOUBLE_PRECISION_REAL
1805 : h1 = _mm_set_pd(hh[nb-3], hh[nb-3]);
1806 : h2 = _mm_set_pd(hh[ldh+nb-2], hh[ldh+nb-2]);
1807 : h3 = _mm_set_pd(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
1808 : #endif
1809 : #ifdef SINGLE_PRECISION_REAL
1810 : h1 = _mm_set_ps(hh[nb-3], hh[nb-3]);
1811 : h2 = _mm_set_ps(hh[ldh+nb-2], hh[ldh+nb-2]);
1812 : h3 = _mm_set_ps(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
1813 : #endif
1814 : #endif
1815 :
1816 :
1817 637184 : q1 = _SSE_LOAD(&q[nb*ldq]);
1818 637184 : q2 = _SSE_LOAD(&q[(nb*ldq)+offset]);
1819 :
1820 1911552 : q1 = _SSE_SUB(q1, _SSE_ADD(_SSE_ADD(_SSE_MUL(z1, h3), _SSE_MUL(y1, h2)) , _SSE_MUL(x1, h1)));
1821 1911552 : q2 = _SSE_SUB(q2, _SSE_ADD(_SSE_ADD(_SSE_MUL(z2, h3), _SSE_MUL(y2, h2)) , _SSE_MUL(x2, h1)));
1822 :
1823 318592 : _SSE_STORE(&q[nb*ldq],q1);
1824 318592 : _SSE_STORE(&q[(nb*ldq)+offset],q2);
1825 :
1826 : #ifdef HAVE_SSE_INTRINSICS
1827 : #ifdef DOUBLE_PRECISION_REAL
1828 536576 : h1 = _mm_set1_pd(hh[nb-2]);
1829 536576 : h2 = _mm_set1_pd(hh[ldh+nb-1]);
1830 : #endif
1831 : #ifdef SINGLE_PRECISION_REAL
1832 100608 : h1 = _mm_set1_ps(hh[nb-2]);
1833 100608 : h2 = _mm_set1_ps(hh[ldh+nb-1]);
1834 : #endif
1835 : #endif
1836 :
1837 : #ifdef HAVE_SPARC64_SSE
1838 : #ifdef DOUBLE_PRECISION_REAL
1839 : h1 = _mm_set_pd(hh[nb-2], hh[nb-2]);
1840 : h2 = _mm_set_pd(hh[ldh+nb-1], hh[ldh+nb-1]);
1841 : #endif
1842 : #ifdef SINGLE_PRECISION_REAL
1843 : h1 = _mm_set_ps(hh[nb-2], hh[nb-2]);
1844 : h2 = _mm_set_ps(hh[ldh+nb-1], hh[ldh+nb-1]);
1845 : #endif
1846 : #endif
1847 :
1848 :
1849 637184 : q1 = _SSE_LOAD(&q[(nb+1)*ldq]);
1850 637184 : q2 = _SSE_LOAD(&q[((nb+1)*ldq)+offset]);
1851 :
1852 1274368 : q1 = _SSE_SUB(q1, _SSE_ADD( _SSE_MUL(y1, h2) , _SSE_MUL(x1, h1)));
1853 1274368 : q2 = _SSE_SUB(q2, _SSE_ADD( _SSE_MUL(y2, h2) , _SSE_MUL(x2, h1)));
1854 :
1855 318592 : _SSE_STORE(&q[(nb+1)*ldq],q1);
1856 318592 : _SSE_STORE(&q[((nb+1)*ldq)+offset],q2);
1857 :
1858 : #ifdef HAVE_SSE_INTRINSICS
1859 : #ifdef DOUBLE_PRECISION_REAL
1860 536576 : h1 = _mm_set1_pd(hh[nb-1]);
1861 : #endif
1862 : #ifdef SINGLE_PRECISION_REAL
1863 100608 : h1 = _mm_set1_ps(hh[nb-1]);
1864 : #endif
1865 : #endif
1866 :
1867 : #ifdef HAVE_SPARC64_SSE
1868 : #ifdef DOUBLE_PRECISION_REAL
1869 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
1870 : #endif
1871 : #ifdef SINGLE_PRECISION_REAL
1872 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
1873 : #endif
1874 : #endif
1875 :
1876 637184 : q1 = _SSE_LOAD(&q[(nb+2)*ldq]);
1877 637184 : q2 = _SSE_LOAD(&q[((nb+2)*ldq)+offset]);
1878 :
1879 637184 : q1 = _SSE_SUB(q1, _SSE_MUL(x1, h1));
1880 637184 : q2 = _SSE_SUB(q2, _SSE_MUL(x2, h1));
1881 :
1882 318592 : _SSE_STORE(&q[(nb+2)*ldq],q1);
1883 318592 : _SSE_STORE(&q[((nb+2)*ldq)+offset],q2);
1884 : }
1885 : /**
1886 : * Unrolled kernel that computes
1887 : #ifdef DOUBLE_PRECISION_REAL
1888 : * 2 rows of Q simultaneously, a
1889 : #endif
1890 : #ifdef SINGLE_PRECISION_REAL
1891 : * 4 rows of Q simultaneously, a
1892 : #endif
1893 : * matrix Vector product with two householder
1894 : * vectors + a rank 1 update is performed
1895 : */
1896 : #ifdef HAVE_SSE_INTRINSICS
1897 : #ifdef DOUBLE_PRECISION_REAL
1898 : __forceinline void hh_trafo_kernel_2_SSE_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4)
1899 : #endif
1900 : #ifdef SINGLE_PRECISION_REAL
1901 : __forceinline void hh_trafo_kernel_4_SSE_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4)
1902 : #endif
1903 : #endif
1904 : #ifdef HAVE_SPARC64_SSE
1905 : #ifdef DOUBLE_PRECISION_REAL
1906 : __forceinline void hh_trafo_kernel_2_SPARC64_4hv_double(double* q, double* hh, int nb, int ldq, int ldh, double s_1_2, double s_1_3, double s_2_3, double s_1_4, double s_2_4, double s_3_4)
1907 : #endif
1908 : #ifdef SINGLE_PRECISION_REAL
1909 : __forceinline void hh_trafo_kernel_4_SPARC64_4hv_single(float* q, float* hh, int nb, int ldq, int ldh, float s_1_2, float s_1_3, float s_2_3, float s_1_4, float s_2_4, float s_3_4)
1910 : #endif
1911 : #endif
1912 :
1913 : {
1914 : /////////////////////////////////////////////////////
1915 : // Matrix Vector Multiplication, Q [2 x nb+3] * hh
1916 : // hh contains four householder vectors
1917 : /////////////////////////////////////////////////////
1918 : int i;
1919 :
1920 0 : __SSE_DATATYPE a1_1 = _SSE_LOAD(&q[ldq*3]);
1921 0 : __SSE_DATATYPE a2_1 = _SSE_LOAD(&q[ldq*2]);
1922 0 : __SSE_DATATYPE a3_1 = _SSE_LOAD(&q[ldq]);
1923 0 : __SSE_DATATYPE a4_1 = _SSE_LOAD(&q[0]);
1924 :
1925 : #ifdef HAVE_SSE_INTRINSICS
1926 : #ifdef DOUBLE_PRECISION_REAL
1927 0 : __SSE_DATATYPE h_2_1 = _mm_set1_pd(hh[ldh+1]);
1928 0 : __SSE_DATATYPE h_3_2 = _mm_set1_pd(hh[(ldh*2)+1]);
1929 0 : __SSE_DATATYPE h_3_1 = _mm_set1_pd(hh[(ldh*2)+2]);
1930 0 : __SSE_DATATYPE h_4_3 = _mm_set1_pd(hh[(ldh*3)+1]);
1931 0 : __SSE_DATATYPE h_4_2 = _mm_set1_pd(hh[(ldh*3)+2]);
1932 0 : __SSE_DATATYPE h_4_1 = _mm_set1_pd(hh[(ldh*3)+3]);
1933 : #endif
1934 : #ifdef SINGLE_PRECISION_REAL
1935 0 : __m128 h_2_1 = _mm_set1_ps(hh[ldh+1] ); // h_2_1 contains four times hh[ldh+1]
1936 0 : __m128 h_3_2 = _mm_set1_ps(hh[(ldh*2)+1]);
1937 0 : __m128 h_3_1 = _mm_set1_ps(hh[(ldh*2)+2]);
1938 0 : __m128 h_4_3 = _mm_set1_ps(hh[(ldh*3)+1]);
1939 0 : __m128 h_4_2 = _mm_set1_ps(hh[(ldh*3)+2]);
1940 0 : __m128 h_4_1 = _mm_set1_ps(hh[(ldh*3)+3]);
1941 : #endif
1942 : #endif
1943 :
1944 : #ifdef HAVE_SPARC64_SSE
1945 : #ifdef DOUBLE_PRECISION_REAL
1946 : __SSE_DATATYPE h_2_1 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
1947 : __SSE_DATATYPE h_3_2 = _mm_set_pd(hh[(ldh*2)+1], hh[(ldh*2)+1]);
1948 : __SSE_DATATYPE h_3_1 = _mm_set_pd(hh[(ldh*2)+2], hh[(ldh*2)+2]);
1949 : __SSE_DATATYPE h_4_3 = _mm_set_pd(hh[(ldh*3)+1], hh[(ldh*3)+1]);
1950 : __SSE_DATATYPE h_4_2 = _mm_set_pd(hh[(ldh*3)+2], hh[(ldh*3)+2]);
1951 : __SSE_DATATYPE h_4_1 = _mm_set_pd(hh[(ldh*3)+3], hh[(ldh*3)+3]);
1952 : #endif
1953 : #ifdef SINGLE_PRECISION_REAL
1954 : __m128 h_2_1 = _mm_set_ps(hh[ldh+1], hh[ldh+1]); // h_2_1 contains four times hh[ldh+1]
1955 : __m128 h_3_2 = _mm_set_ps(hh[(ldh*2)+1], hh[(ldh*2)+1]);
1956 : __m128 h_3_1 = _mm_set_ps(hh[(ldh*2)+2], hh[(ldh*2)+2]);
1957 : __m128 h_4_3 = _mm_set_ps(hh[(ldh*3)+1], hh[(ldh*3)+1]);
1958 : __m128 h_4_2 = _mm_set_ps(hh[(ldh*3)+2], hh[(ldh*3)+2]);
1959 : __m128 h_4_1 = _mm_set_ps(hh[(ldh*3)+3], hh[(ldh*3)+3]);
1960 : #endif
1961 : #endif
1962 :
1963 0 : __SSE_DATATYPE w1 = _SSE_ADD(a4_1, _SSE_MUL(a3_1, h_4_3));
1964 0 : w1 = _SSE_ADD(w1, _SSE_MUL(a2_1, h_4_2));
1965 0 : w1 = _SSE_ADD(w1, _SSE_MUL(a1_1, h_4_1));
1966 0 : __SSE_DATATYPE z1 = _SSE_ADD(a3_1, _SSE_MUL(a2_1, h_3_2));
1967 0 : z1 = _SSE_ADD(z1, _SSE_MUL(a1_1, h_3_1));
1968 0 : __SSE_DATATYPE y1 = _SSE_ADD(a2_1, _SSE_MUL(a1_1, h_2_1));
1969 0 : __SSE_DATATYPE x1 = a1_1;
1970 :
1971 : __SSE_DATATYPE q1;
1972 :
1973 : __SSE_DATATYPE h1;
1974 : __SSE_DATATYPE h2;
1975 : __SSE_DATATYPE h3;
1976 : __SSE_DATATYPE h4;
1977 :
1978 0 : for(i = 4; i < nb; i++)
1979 : {
1980 : #ifdef HAVE_SSE_INTRINSICS
1981 : #ifdef DOUBLE_PRECISION_REAL
1982 0 : h1 = _mm_set1_pd(hh[i-3]);
1983 0 : h2 = _mm_set1_pd(hh[ldh+i-2]);
1984 0 : h3 = _mm_set1_pd(hh[(ldh*2)+i-1]);
1985 0 : h4 = _mm_set1_pd(hh[(ldh*3)+i]);
1986 : #endif
1987 : #ifdef SINGLE_PRECISION_REAL
1988 0 : h1 = _mm_set1_ps(hh[i-3]);
1989 0 : h2 = _mm_set1_ps(hh[ldh+i-2]);
1990 0 : h3 = _mm_set1_ps(hh[(ldh*2)+i-1]);
1991 0 : h4 = _mm_set1_ps(hh[(ldh*3)+i]);
1992 : #endif
1993 : #endif
1994 :
1995 : #ifdef HAVE_SPARC64_SSE
1996 : #ifdef DOUBLE_PRECISION_REAL
1997 : h1 = _mm_set_pd(hh[i-3], hh[i-3]);
1998 : h2 = _mm_set_pd(hh[ldh+i-2], hh[ldh+i-2]);
1999 : h3 = _mm_set_pd(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
2000 : h4 = _mm_set_pd(hh[(ldh*3)+i], hh[(ldh*3)+i]);
2001 : #endif
2002 : #ifdef SINGLE_PRECISION_REAL
2003 : h1 = _mm_set_ps(hh[i-3], hh[i-3]);
2004 : h2 = _mm_set_ps(hh[ldh+i-2], hh[ldh+i-2]);
2005 : h3 = _mm_set_ps(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
2006 : h4 = _mm_set_ps(hh[(ldh*3)+i], hh[(ldh*3)+i]);
2007 : #endif
2008 : #endif
2009 :
2010 0 : q1 = _SSE_LOAD(&q[i*ldq]);
2011 :
2012 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
2013 0 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
2014 0 : z1 = _SSE_ADD(z1, _SSE_MUL(q1,h3));
2015 0 : w1 = _SSE_ADD(w1, _SSE_MUL(q1,h4));
2016 : }
2017 :
2018 : #ifdef HAVE_SSE_INTRINSICS
2019 : #ifdef DOUBLE_PRECISION_REAL
2020 0 : h1 = _mm_set1_pd(hh[nb-3]);
2021 0 : h2 = _mm_set1_pd(hh[ldh+nb-2]);
2022 0 : h3 = _mm_set1_pd(hh[(ldh*2)+nb-1]);
2023 : #endif
2024 : #ifdef SINGLE_PRECISION_REAL
2025 0 : h1 = _mm_set1_ps(hh[nb-3]);
2026 0 : h2 = _mm_set1_ps(hh[ldh+nb-2]);
2027 0 : h3 = _mm_set1_ps(hh[(ldh*2)+nb-1]);
2028 : #endif
2029 : #endif
2030 :
2031 : #ifdef HAVE_SPARC64_SSE
2032 : #ifdef DOUBLE_PRECISION_REAL
2033 : h1 = _mm_set_pd(hh[nb-3], hh[nb-3]);
2034 : h2 = _mm_set_pd(hh[ldh+nb-2], hh[ldh+nb-2]);
2035 : h3 = _mm_set_pd(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
2036 : #endif
2037 : #ifdef SINGLE_PRECISION_REAL
2038 : h1 = _mm_set_ps(hh[nb-3], hh[nb-3]);
2039 : h2 = _mm_set_ps(hh[ldh+nb-2], hh[ldh+nb-2]);
2040 : h3 = _mm_set_ps(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
2041 : #endif
2042 : #endif
2043 :
2044 0 : q1 = _SSE_LOAD(&q[nb*ldq]);
2045 :
2046 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
2047 0 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
2048 0 : z1 = _SSE_ADD(z1, _SSE_MUL(q1,h3));
2049 :
2050 : #ifdef HAVE_SSE_INTRINSICS
2051 : #ifdef DOUBLE_PRECISION_REAL
2052 0 : h1 = _mm_set1_pd(hh[nb-2]);
2053 0 : h2 = _mm_set1_pd(hh[(ldh*1)+nb-1]);
2054 : #endif
2055 : #ifdef SINGLE_PRECISION_REAL
2056 0 : h1 = _mm_set1_ps(hh[nb-2]);
2057 0 : h2 = _mm_set1_ps(hh[(ldh*1)+nb-1]);
2058 : #endif
2059 : #endif
2060 :
2061 : #ifdef HAVE_SPARC64_SSE
2062 : #ifdef DOUBLE_PRECISION_REAL
2063 : h1 = _mm_set_pd(hh[nb-2], hh[nb-2]);
2064 : h2 = _mm_set_pd(hh[(ldh*1)+nb-1], hh[(ldh*1)+nb-1]);
2065 : #endif
2066 : #ifdef SINGLE_PRECISION_REAL
2067 : h1 = _mm_set_ps(hh[nb-2], hh[nb-2]);
2068 : h2 = _mm_set_ps(hh[(ldh*1)+nb-1], hh[(ldh*1)+nb-1]);
2069 : #endif
2070 : #endif
2071 :
2072 0 : q1 = _SSE_LOAD(&q[(nb+1)*ldq]);
2073 :
2074 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
2075 0 : y1 = _SSE_ADD(y1, _SSE_MUL(q1,h2));
2076 :
2077 : #ifdef HAVE_SSE_INTRINSICS
2078 : #ifdef DOUBLE_PRECISION_REAL
2079 0 : h1 = _mm_set1_pd(hh[nb-1]);
2080 : #endif
2081 : #ifdef SINGLE_PRECISION_REAL
2082 0 : h1 = _mm_set1_ps(hh[nb-1]);
2083 : #endif
2084 : #endif
2085 :
2086 : #ifdef HAVE_SPARC64_SSE
2087 : #ifdef DOUBLE_PRECISION_REAL
2088 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
2089 : #endif
2090 : #ifdef SINGLE_PRECISION_REAL
2091 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
2092 : #endif
2093 : #endif
2094 :
2095 0 : q1 = _SSE_LOAD(&q[(nb+2)*ldq]);
2096 :
2097 0 : x1 = _SSE_ADD(x1, _SSE_MUL(q1,h1));
2098 : /////////////////////////////////////////////////////
2099 : // Rank-1 update of Q [2 x nb+3]
2100 : /////////////////////////////////////////////////////
2101 :
2102 : #ifdef HAVE_SSE_INTRINSICS
2103 : #ifdef DOUBLE_PRECISION_REAL
2104 0 : __SSE_DATATYPE tau1 = _mm_set1_pd(hh[0]);
2105 0 : __SSE_DATATYPE tau2 = _mm_set1_pd(hh[ldh]);
2106 0 : __SSE_DATATYPE tau3 = _mm_set1_pd(hh[ldh*2]);
2107 0 : __SSE_DATATYPE tau4 = _mm_set1_pd(hh[ldh*3]);
2108 :
2109 0 : __SSE_DATATYPE vs_1_2 = _mm_set1_pd(s_1_2);
2110 0 : __SSE_DATATYPE vs_1_3 = _mm_set1_pd(s_1_3);
2111 0 : __SSE_DATATYPE vs_2_3 = _mm_set1_pd(s_2_3);
2112 0 : __SSE_DATATYPE vs_1_4 = _mm_set1_pd(s_1_4);
2113 0 : __SSE_DATATYPE vs_2_4 = _mm_set1_pd(s_2_4);
2114 0 : __SSE_DATATYPE vs_3_4 = _mm_set1_pd(s_3_4);
2115 : #endif
2116 : #ifdef SINGLE_PRECISION_REAL
2117 0 : __m128 tau1 = _mm_set1_ps(hh[0]);
2118 0 : __m128 tau2 = _mm_set1_ps(hh[ldh]);
2119 0 : __m128 tau3 = _mm_set1_ps(hh[ldh*2]);
2120 0 : __m128 tau4 = _mm_set1_ps(hh[ldh*3]);
2121 :
2122 0 : __m128 vs_1_2 = _mm_set1_ps(s_1_2);
2123 0 : __m128 vs_1_3 = _mm_set1_ps(s_1_3);
2124 0 : __m128 vs_2_3 = _mm_set1_ps(s_2_3);
2125 0 : __m128 vs_1_4 = _mm_set1_ps(s_1_4);
2126 0 : __m128 vs_2_4 = _mm_set1_ps(s_2_4);
2127 0 : __m128 vs_3_4 = _mm_set1_ps(s_3_4);
2128 : #endif
2129 : #endif
2130 :
2131 : #ifdef HAVE_SPARC64_SSE
2132 : #ifdef DOUBLE_PRECISION_REAL
2133 : __SSE_DATATYPE tau1 = _mm_set_pd(hh[0], hh[0]);
2134 : __SSE_DATATYPE tau2 = _mm_set_pd(hh[ldh], hh[ldh]);
2135 : __SSE_DATATYPE tau3 = _mm_set_pd(hh[ldh*2], hh[ldh*2]);
2136 : __SSE_DATATYPE tau4 = _mm_set_pd(hh[ldh*3], hh[ldh*3]);
2137 :
2138 : __SSE_DATATYPE vs_1_2 = _mm_set_pd(s_1_2, s_1_2);
2139 : __SSE_DATATYPE vs_1_3 = _mm_set_pd(s_1_3, s_1_3);
2140 : __SSE_DATATYPE vs_2_3 = _mm_set_pd(s_2_3, s_2_3);
2141 : __SSE_DATATYPE vs_1_4 = _mm_set_pd(s_1_4, s_1_4);
2142 : __SSE_DATATYPE vs_2_4 = _mm_set_pd(s_2_4, s_2_4);
2143 : __SSE_DATATYPE vs_3_4 = _mm_set_pd(s_3_4, s_3_4);
2144 : #endif
2145 : #ifdef SINGLE_PRECISION_REAL
2146 : __m128 tau1 = _mm_set1_ps(hh[0], hh[0]);
2147 : __m128 tau2 = _mm_set1_ps(hh[ldh], hh[ldh]);
2148 : __m128 tau3 = _mm_set1_ps(hh[ldh*2], hh[ldh*2]);
2149 : __m128 tau4 = _mm_set1_ps(hh[ldh*3], hh[ldh*3]);
2150 :
2151 : __m128 vs_1_2 = _mm_set1_ps(s_1_2, s_1_2);
2152 : __m128 vs_1_3 = _mm_set1_ps(s_1_3, s_1_3);
2153 : __m128 vs_2_3 = _mm_set1_ps(s_2_3, s_2_3);
2154 : __m128 vs_1_4 = _mm_set1_ps(s_1_4, s_1_4);
2155 : __m128 vs_2_4 = _mm_set1_ps(s_2_4, s_2_4);
2156 : __m128 vs_3_4 = _mm_set1_ps(s_3_4, s_3_4);
2157 : #endif
2158 : #endif
2159 :
2160 :
2161 0 : h1 = tau1;
2162 0 : x1 = _SSE_MUL(x1, h1);
2163 :
2164 0 : h1 = tau2;
2165 0 : h2 = _SSE_MUL(h1, vs_1_2);
2166 :
2167 0 : y1 = _SSE_SUB(_SSE_MUL(y1,h1), _SSE_MUL(x1,h2));
2168 :
2169 0 : h1 = tau3;
2170 0 : h2 = _SSE_MUL(h1, vs_1_3);
2171 0 : h3 = _SSE_MUL(h1, vs_2_3);
2172 :
2173 0 : z1 = _SSE_SUB(_SSE_MUL(z1,h1), _SSE_ADD(_SSE_MUL(y1,h3), _SSE_MUL(x1,h2)));
2174 :
2175 0 : h1 = tau4;
2176 0 : h2 = _SSE_MUL(h1, vs_1_4);
2177 0 : h3 = _SSE_MUL(h1, vs_2_4);
2178 0 : h4 = _SSE_MUL(h1, vs_3_4);
2179 :
2180 0 : w1 = _SSE_SUB(_SSE_MUL(w1,h1), _SSE_ADD(_SSE_MUL(z1,h4), _SSE_ADD(_SSE_MUL(y1,h3), _SSE_MUL(x1,h2))));
2181 :
2182 0 : q1 = _SSE_LOAD(&q[0]);
2183 0 : q1 = _SSE_SUB(q1, w1);
2184 : _SSE_STORE(&q[0],q1);
2185 :
2186 : #ifdef HAVE_SSE_INTRINSICS
2187 : #ifdef DOUBLE_PRECISION_REAL
2188 0 : h4 = _mm_set1_pd(hh[(ldh*3)+1]);
2189 : #endif
2190 : #ifdef SINGLE_PRECISION_REAL
2191 0 : h4 = _mm_set1_ps(hh[(ldh*3)+1]);
2192 : #endif
2193 : #endif
2194 :
2195 : #ifdef HAVE_SPARC64_SSE
2196 : #ifdef DOUBLE_PRECISION_REAL
2197 : h4 = _mm_set_pd(hh[(ldh*3)+1], hh[(ldh*3)+1]);
2198 : #endif
2199 : #ifdef SINGLE_PRECISION_REAL
2200 : h4 = _mm_set_ps(hh[(ldh*3)+1], hh[(ldh*3)+1]);
2201 : #endif
2202 : #endif
2203 :
2204 0 : q1 = _SSE_LOAD(&q[ldq]);
2205 :
2206 0 : q1 = _SSE_SUB(q1, _SSE_ADD(z1, _SSE_MUL(w1, h4)));
2207 :
2208 0 : _SSE_STORE(&q[ldq],q1);
2209 :
2210 : #ifdef HAVE_SSE_INTRINSICS
2211 : #ifdef DOUBLE_PRECISION_REAL
2212 0 : h3 = _mm_set1_pd(hh[(ldh*2)+1]);
2213 0 : h4 = _mm_set1_pd(hh[(ldh*3)+2]);
2214 : #endif
2215 : #ifdef SINGLE_PRECISION_REAL
2216 0 : h3 = _mm_set1_ps(hh[(ldh*2)+1]);
2217 0 : h4 = _mm_set1_ps(hh[(ldh*3)+2]);
2218 : #endif
2219 : #endif
2220 :
2221 : #ifdef HAVE_SPARC64_SSE
2222 : #ifdef DOUBLE_PRECISION_REAL
2223 : h3 = _mm_set_pd(hh[(ldh*2)+1], hh[(ldh*2)+1]);
2224 : h4 = _mm_set_pd(hh[(ldh*3)+2], hh[(ldh*3)+2]);
2225 : #endif
2226 : #ifdef SINGLE_PRECISION_REAL
2227 : h3 = _mm_set_ps(hh[(ldh*2)+1], hh[(ldh*2)+1]);
2228 : h4 = _mm_set_ps(hh[(ldh*3)+2], hh[(ldh*3)+2]);
2229 : #endif
2230 : #endif
2231 :
2232 0 : q1 = _SSE_LOAD(&q[ldq*2]);
2233 :
2234 0 : q1 = _SSE_SUB(q1, _SSE_ADD(y1, _SSE_ADD(_SSE_MUL(z1, h3), _SSE_MUL(w1, h4))));
2235 :
2236 0 : _SSE_STORE(&q[ldq*2],q1);
2237 :
2238 : #ifdef HAVE_SSE_INTRINSICS
2239 : #ifdef DOUBLE_PRECISION_REAL
2240 0 : h2 = _mm_set1_pd(hh[ldh+1]);
2241 0 : h3 = _mm_set1_pd(hh[(ldh*2)+2]);
2242 0 : h4 = _mm_set1_pd(hh[(ldh*3)+3]);
2243 : #endif
2244 : #ifdef SINGLE_PRECISION_REAL
2245 0 : h2 = _mm_set1_ps(hh[ldh+1]);
2246 0 : h3 = _mm_set1_ps(hh[(ldh*2)+2]);
2247 0 : h4 = _mm_set1_ps(hh[(ldh*3)+3]);
2248 : #endif
2249 : #endif
2250 :
2251 : #ifdef HAVE_SPARC64_SSE
2252 : #ifdef DOUBLE_PRECISION_REAL
2253 : h2 = _mm_set_pd(hh[ldh+1], hh[ldh+1]);
2254 : h3 = _mm_set_pd(hh[(ldh*2)+2], hh[(ldh*2)+2]);
2255 : h4 = _mm_set_pd(hh[(ldh*3)+3], hh[(ldh*3)+3]);
2256 : #endif
2257 : #ifdef SINGLE_PRECISION_REAL
2258 : h2 = _mm_set_ps(hh[ldh+1], hh[ldh+1]);
2259 : h3 = _mm_set_ps(hh[(ldh*2)+2], hh[(ldh*2)+2]);
2260 : h4 = _mm_set_ps(hh[(ldh*3)+3], hh[(ldh*3)+3]);
2261 : #endif
2262 : #endif
2263 :
2264 0 : q1 = _SSE_LOAD(&q[ldq*3]);
2265 :
2266 0 : q1 = _SSE_SUB(q1, _SSE_ADD(x1, _SSE_ADD(_SSE_MUL(y1, h2), _SSE_ADD(_SSE_MUL(z1, h3), _SSE_MUL(w1, h4)))));
2267 :
2268 0 : _SSE_STORE(&q[ldq*3], q1);
2269 :
2270 0 : for (i = 4; i < nb; i++)
2271 : {
2272 : #ifdef HAVE_SSE_INTRINSICS
2273 : #ifdef DOUBLE_PRECISION_REAL
2274 0 : h1 = _mm_set1_pd(hh[i-3]);
2275 0 : h2 = _mm_set1_pd(hh[ldh+i-2]);
2276 0 : h3 = _mm_set1_pd(hh[(ldh*2)+i-1]);
2277 0 : h4 = _mm_set1_pd(hh[(ldh*3)+i]);
2278 : #endif
2279 : #ifdef SINGLE_PRECISION_REAL
2280 0 : h1 = _mm_set1_ps(hh[i-3]);
2281 0 : h2 = _mm_set1_ps(hh[ldh+i-2]);
2282 0 : h3 = _mm_set1_ps(hh[(ldh*2)+i-1]);
2283 0 : h4 = _mm_set1_ps(hh[(ldh*3)+i]);
2284 : #endif
2285 : #endif
2286 :
2287 : #ifdef HAVE_SPARC64_SSE
2288 : #ifdef DOUBLE_PRECISION_REAL
2289 : h1 = _mm_set_pd(hh[i-3], hh[i-3]);
2290 : h2 = _mm_set_pd(hh[ldh+i-2], hh[ldh+i-2]);
2291 : h3 = _mm_set_pd(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
2292 : h4 = _mm_set_pd(hh[(ldh*3)+i], hh[(ldh*3)+i]);
2293 : #endif
2294 : #ifdef SINGLE_PRECISION_REAL
2295 : h1 = _mm_set_ps(hh[i-3], hh[i-3]);
2296 : h2 = _mm_set_ps(hh[ldh+i-2], hh[ldh+i-2]);
2297 : h3 = _mm_set_ps(hh[(ldh*2)+i-1], hh[(ldh*2)+i-1]);
2298 : h4 = _mm_set_ps(hh[(ldh*3)+i], hh[(ldh*3)+i]);
2299 : #endif
2300 : #endif
2301 :
2302 0 : q1 = _SSE_LOAD(&q[i*ldq]);
2303 :
2304 0 : q1 = _SSE_SUB(q1, _SSE_ADD(_SSE_ADD(_SSE_MUL(w1, h4), _SSE_MUL(z1, h3)), _SSE_ADD(_SSE_MUL(x1,h1), _SSE_MUL(y1, h2))));
2305 :
2306 0 : _SSE_STORE(&q[i*ldq],q1);
2307 : }
2308 :
2309 : #ifdef HAVE_SSE_INTRINSICS
2310 : #ifdef DOUBLE_PRECISION_REAL
2311 0 : h1 = _mm_set1_pd(hh[nb-3]);
2312 0 : h2 = _mm_set1_pd(hh[ldh+nb-2]);
2313 0 : h3 = _mm_set1_pd(hh[(ldh*2)+nb-1]);
2314 : #endif
2315 : #ifdef SINGLE_PRECISION_REAL
2316 0 : h1 = _mm_set1_ps(hh[nb-3]);
2317 0 : h2 = _mm_set1_ps(hh[ldh+nb-2]);
2318 0 : h3 = _mm_set1_ps(hh[(ldh*2)+nb-1]);
2319 : #endif
2320 : #endif
2321 :
2322 : #ifdef HAVE_SPARC64_SSE
2323 : #ifdef DOUBLE_PRECISION_REAL
2324 : h1 = _mm_set_pd(hh[nb-3], hh[nb-3]);
2325 : h2 = _mm_set_pd(hh[ldh+nb-2], hh[ldh+nb-2]);
2326 : h3 = _mm_set_pd(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
2327 : #endif
2328 : #ifdef SINGLE_PRECISION_REAL
2329 : h1 = _mm_set_ps(hh[nb-3], hh[nb-3]);
2330 : h2 = _mm_set_ps(hh[ldh+nb-2], hh[ldh+nb-2]);
2331 : h3 = _mm_set_ps(hh[(ldh*2)+nb-1], hh[(ldh*2)+nb-1]);
2332 : #endif
2333 : #endif
2334 :
2335 :
2336 0 : q1 = _SSE_LOAD(&q[nb*ldq]);
2337 0 : q1 = _SSE_SUB(q1, _SSE_ADD(_SSE_ADD(_SSE_MUL(z1, h3), _SSE_MUL(y1, h2)) , _SSE_MUL(x1, h1)));
2338 :
2339 0 : _SSE_STORE(&q[nb*ldq],q1);
2340 :
2341 : #ifdef HAVE_SSE_INTRINSICS
2342 : #ifdef DOUBLE_PRECISION_REAL
2343 0 : h1 = _mm_set1_pd(hh[nb-2]);
2344 0 : h2 = _mm_set1_pd(hh[ldh+nb-1]);
2345 : #endif
2346 : #ifdef SINGLE_PRECISION_REAL
2347 0 : h1 = _mm_set1_ps(hh[nb-2]);
2348 0 : h2 = _mm_set1_ps(hh[ldh+nb-1]);
2349 : #endif
2350 : #endif
2351 :
2352 : #ifdef HAVE_SPARC64_SSE
2353 : #ifdef DOUBLE_PRECISION_REAL
2354 : h1 = _mm_set_pd(hh[nb-2], hh[nb-2]);
2355 : h2 = _mm_set_pd(hh[ldh+nb-1], hh[ldh+nb-1]);
2356 : #endif
2357 : #ifdef SINGLE_PRECISION_REAL
2358 : h1 = _mm_set_ps(hh[nb-2], hh[nb-2]);
2359 : h2 = _mm_set_ps(hh[ldh+nb-1], hh[ldh+nb-1]);
2360 : #endif
2361 : #endif
2362 :
2363 0 : q1 = _SSE_LOAD(&q[(nb+1)*ldq]);
2364 :
2365 0 : q1 = _SSE_SUB(q1, _SSE_ADD( _SSE_MUL(y1, h2) , _SSE_MUL(x1, h1)));
2366 :
2367 0 : _SSE_STORE(&q[(nb+1)*ldq],q1);
2368 :
2369 : #ifdef HAVE_SSE_INTRINSICS
2370 : #ifdef DOUBLE_PRECISION_REAL
2371 0 : h1 = _mm_set1_pd(hh[nb-1]);
2372 : #endif
2373 : #ifdef SINGLE_PRECISION_REAL
2374 0 : h1 = _mm_set1_ps(hh[nb-1]);
2375 : #endif
2376 : #endif
2377 :
2378 : #ifdef HAVE_SPARC64_SSE
2379 : #ifdef DOUBLE_PRECISION_REAL
2380 : h1 = _mm_set_pd(hh[nb-1], hh[nb-1]);
2381 : #endif
2382 : #ifdef SINGLE_PRECISION_REAL
2383 : h1 = _mm_set_ps(hh[nb-1], hh[nb-1]);
2384 : #endif
2385 : #endif
2386 :
2387 0 : q1 = _SSE_LOAD(&q[(nb+2)*ldq]);
2388 :
2389 0 : q1 = _SSE_SUB(q1, _SSE_MUL(x1, h1));
2390 :
2391 0 : _SSE_STORE(&q[(nb+2)*ldq],q1);
2392 : }
2393 :
|
