1 #include <blas_extended.h>
2 #include <blas_extended_private.h>
3 #include <blas_fpu.h>
BLAS_zhemv2_c_c_x(enum blas_order_type order,enum blas_uplo_type uplo,int n,const void * alpha,const void * a,int lda,const void * x_head,const void * x_tail,int incx,const void * beta,const void * y,int incy,enum blas_prec_type prec)4 void BLAS_zhemv2_c_c_x(enum blas_order_type order, enum blas_uplo_type uplo,
5 int n, const void *alpha, const void *a, int lda,
6 const void *x_head, const void *x_tail, int incx,
7 const void *beta, const void *y, int incy,
8 enum blas_prec_type prec)
9
10 /*
11 * Purpose
12 * =======
13 *
14 * This routines computes the matrix product:
15 *
16 * y <- alpha * A * (x_head + x_tail) + beta * y
17 *
18 * where A is a complex Hermitian matrix.
19 *
20 * Arguments
21 * =========
22 *
23 * order (input) enum blas_order_type
24 * Storage format of input symmetric matrix A.
25 *
26 * uplo (input) enum blas_uplo_type
27 * Determines which half of matrix A (upper or lower triangle)
28 * is accessed.
29 *
30 * n (input) int
31 * Dimension of A and size of vectors x, y.
32 *
33 * alpha (input) const void*
34 *
35 * a (input) void*
36 * Matrix A.
37 *
38 * lda (input) int
39 * Leading dimension of matrix A.
40 *
41 * x_head (input) void*
42 * Vector x_head
43 *
44 * x_tail (input) void*
45 * Vector x_tail
46 *
47 * incx (input) int
48 * Stride for vector x.
49 *
50 * beta (input) const void*
51 *
52 * y (input) void*
53 * Vector y.
54 *
55 * incy (input) int
56 * Stride for vector y.
57 *
58 * prec (input) enum blas_prec_type
59 * Specifies the internal precision to be used.
60 * = blas_prec_single: single precision.
61 * = blas_prec_double: double precision.
62 * = blas_prec_extra : anything at least 1.5 times as accurate
63 * than double, and wider than 80-bits.
64 * We use double-double in our implementation.
65 *
66 */
67 {
68 /* Routine name */
69 const char routine_name[] = "BLAS_zhemv2_c_c_x";
70 switch (prec) {
71
72 case blas_prec_single:{
73
74 int i, j;
75 int xi, yi, xi0, yi0;
76 int aij, ai;
77 int incai;
78 int incaij, incaij2;
79
80 const float *a_i = (float *) a;
81 const float *x_head_i = (float *) x_head;
82 const float *x_tail_i = (float *) x_tail;
83 double *y_i = (double *) y;
84 double *alpha_i = (double *) alpha;
85 double *beta_i = (double *) beta;
86 float a_elem[2];
87 float x_elem[2];
88 double y_elem[2];
89 float diag_elem;
90 double prod1[2];
91 double prod2[2];
92 double sum1[2];
93 double sum2[2];
94 double tmp1[2];
95 double tmp2[2];
96 double tmp3[2];
97
98
99
100 /* Test for no-op */
101 if (n <= 0) {
102 return;
103 }
104 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0
105 && (beta_i[0] == 1.0 && beta_i[1] == 0.0)) {
106 return;
107 }
108
109 /* Check for error conditions. */
110 if (n < 0) {
111 BLAS_error(routine_name, -3, n, NULL);
112 }
113 if (lda < n) {
114 BLAS_error(routine_name, -6, n, NULL);
115 }
116 if (incx == 0) {
117 BLAS_error(routine_name, -9, incx, NULL);
118 }
119 if (incy == 0) {
120 BLAS_error(routine_name, -12, incy, NULL);
121 }
122
123 if ((order == blas_colmajor && uplo == blas_upper) ||
124 (order == blas_rowmajor && uplo == blas_lower)) {
125 incai = lda;
126 incaij = 1;
127 incaij2 = lda;
128 } else {
129 incai = 1;
130 incaij = lda;
131 incaij2 = 1;
132 }
133
134 incx *= 2;
135 incy *= 2;
136 incai *= 2;
137 incaij *= 2;
138 incaij2 *= 2;
139 xi0 = (incx > 0) ? 0 : ((-n + 1) * incx);
140 yi0 = (incy > 0) ? 0 : ((-n + 1) * incy);
141
142
143
144 /* The most general form, y <--- alpha * A * (x_head + x_tail) + beta * y */
145 if (uplo == blas_lower) {
146 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
147 sum1[0] = sum1[1] = 0.0;
148 sum2[0] = sum2[1] = 0.0;
149 for (j = 0, aij = ai, xi = xi0; j < i;
150 j++, aij += incaij, xi += incx) {
151 a_elem[0] = a_i[aij];
152 a_elem[1] = a_i[aij + 1];
153 x_elem[0] = x_head_i[xi];
154 x_elem[1] = x_head_i[xi + 1];
155 {
156 prod1[0] =
157 (double) a_elem[0] * x_elem[0] -
158 (double) a_elem[1] * x_elem[1];
159 prod1[1] =
160 (double) a_elem[0] * x_elem[1] +
161 (double) a_elem[1] * x_elem[0];
162 }
163 sum1[0] = sum1[0] + prod1[0];
164 sum1[1] = sum1[1] + prod1[1];
165 x_elem[0] = x_tail_i[xi];
166 x_elem[1] = x_tail_i[xi + 1];
167 {
168 prod2[0] =
169 (double) a_elem[0] * x_elem[0] -
170 (double) a_elem[1] * x_elem[1];
171 prod2[1] =
172 (double) a_elem[0] * x_elem[1] +
173 (double) a_elem[1] * x_elem[0];
174 }
175 sum2[0] = sum2[0] + prod2[0];
176 sum2[1] = sum2[1] + prod2[1];
177 }
178
179 diag_elem = a_i[aij];
180 x_elem[0] = x_head_i[xi];
181 x_elem[1] = x_head_i[xi + 1];
182 {
183 prod1[0] = (double) x_elem[0] * diag_elem;
184 prod1[1] = (double) x_elem[1] * diag_elem;
185 }
186 sum1[0] = sum1[0] + prod1[0];
187 sum1[1] = sum1[1] + prod1[1];
188 x_elem[0] = x_tail_i[xi];
189 x_elem[1] = x_tail_i[xi + 1];
190 {
191 prod2[0] = (double) x_elem[0] * diag_elem;
192 prod2[1] = (double) x_elem[1] * diag_elem;
193 }
194 sum2[0] = sum2[0] + prod2[0];
195 sum2[1] = sum2[1] + prod2[1];
196 j++;
197 aij += incaij2;
198 xi += incx;
199
200 for (; j < n; j++, aij += incaij2, xi += incx) {
201 a_elem[0] = a_i[aij];
202 a_elem[1] = a_i[aij + 1];
203 a_elem[1] = -a_elem[1];
204 x_elem[0] = x_head_i[xi];
205 x_elem[1] = x_head_i[xi + 1];
206 {
207 prod1[0] =
208 (double) a_elem[0] * x_elem[0] -
209 (double) a_elem[1] * x_elem[1];
210 prod1[1] =
211 (double) a_elem[0] * x_elem[1] +
212 (double) a_elem[1] * x_elem[0];
213 }
214 sum1[0] = sum1[0] + prod1[0];
215 sum1[1] = sum1[1] + prod1[1];
216 x_elem[0] = x_tail_i[xi];
217 x_elem[1] = x_tail_i[xi + 1];
218 {
219 prod2[0] =
220 (double) a_elem[0] * x_elem[0] -
221 (double) a_elem[1] * x_elem[1];
222 prod2[1] =
223 (double) a_elem[0] * x_elem[1] +
224 (double) a_elem[1] * x_elem[0];
225 }
226 sum2[0] = sum2[0] + prod2[0];
227 sum2[1] = sum2[1] + prod2[1];
228 }
229 sum1[0] = sum1[0] + sum2[0];
230 sum1[1] = sum1[1] + sum2[1];
231 {
232 tmp1[0] =
233 (double) sum1[0] * alpha_i[0] - (double) sum1[1] * alpha_i[1];
234 tmp1[1] =
235 (double) sum1[0] * alpha_i[1] + (double) sum1[1] * alpha_i[0];
236 }
237 y_elem[0] = y_i[yi];
238 y_elem[1] = y_i[yi + 1];
239 {
240 tmp2[0] =
241 (double) y_elem[0] * beta_i[0] - (double) y_elem[1] * beta_i[1];
242 tmp2[1] =
243 (double) y_elem[0] * beta_i[1] + (double) y_elem[1] * beta_i[0];
244 }
245 tmp3[0] = tmp1[0] + tmp2[0];
246 tmp3[1] = tmp1[1] + tmp2[1];
247 y_i[yi] = tmp3[0];
248 y_i[yi + 1] = tmp3[1];
249 }
250 } else {
251 /* uplo == blas_upper */
252 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
253 sum1[0] = sum1[1] = 0.0;
254 sum2[0] = sum2[1] = 0.0;
255
256 for (j = 0, aij = ai, xi = xi0; j < i;
257 j++, aij += incaij, xi += incx) {
258 a_elem[0] = a_i[aij];
259 a_elem[1] = a_i[aij + 1];
260 a_elem[1] = -a_elem[1];
261 x_elem[0] = x_head_i[xi];
262 x_elem[1] = x_head_i[xi + 1];
263 {
264 prod1[0] =
265 (double) a_elem[0] * x_elem[0] -
266 (double) a_elem[1] * x_elem[1];
267 prod1[1] =
268 (double) a_elem[0] * x_elem[1] +
269 (double) a_elem[1] * x_elem[0];
270 }
271 sum1[0] = sum1[0] + prod1[0];
272 sum1[1] = sum1[1] + prod1[1];
273 x_elem[0] = x_tail_i[xi];
274 x_elem[1] = x_tail_i[xi + 1];
275 {
276 prod2[0] =
277 (double) a_elem[0] * x_elem[0] -
278 (double) a_elem[1] * x_elem[1];
279 prod2[1] =
280 (double) a_elem[0] * x_elem[1] +
281 (double) a_elem[1] * x_elem[0];
282 }
283 sum2[0] = sum2[0] + prod2[0];
284 sum2[1] = sum2[1] + prod2[1];
285 }
286
287 diag_elem = a_i[aij];
288 x_elem[0] = x_head_i[xi];
289 x_elem[1] = x_head_i[xi + 1];
290 {
291 prod1[0] = (double) x_elem[0] * diag_elem;
292 prod1[1] = (double) x_elem[1] * diag_elem;
293 }
294 sum1[0] = sum1[0] + prod1[0];
295 sum1[1] = sum1[1] + prod1[1];
296 x_elem[0] = x_tail_i[xi];
297 x_elem[1] = x_tail_i[xi + 1];
298 {
299 prod2[0] = (double) x_elem[0] * diag_elem;
300 prod2[1] = (double) x_elem[1] * diag_elem;
301 }
302 sum2[0] = sum2[0] + prod2[0];
303 sum2[1] = sum2[1] + prod2[1];
304 j++;
305 aij += incaij2;
306 xi += incx;
307
308 for (; j < n; j++, aij += incaij2, xi += incx) {
309 a_elem[0] = a_i[aij];
310 a_elem[1] = a_i[aij + 1];
311 x_elem[0] = x_head_i[xi];
312 x_elem[1] = x_head_i[xi + 1];
313 {
314 prod1[0] =
315 (double) a_elem[0] * x_elem[0] -
316 (double) a_elem[1] * x_elem[1];
317 prod1[1] =
318 (double) a_elem[0] * x_elem[1] +
319 (double) a_elem[1] * x_elem[0];
320 }
321 sum1[0] = sum1[0] + prod1[0];
322 sum1[1] = sum1[1] + prod1[1];
323 x_elem[0] = x_tail_i[xi];
324 x_elem[1] = x_tail_i[xi + 1];
325 {
326 prod2[0] =
327 (double) a_elem[0] * x_elem[0] -
328 (double) a_elem[1] * x_elem[1];
329 prod2[1] =
330 (double) a_elem[0] * x_elem[1] +
331 (double) a_elem[1] * x_elem[0];
332 }
333 sum2[0] = sum2[0] + prod2[0];
334 sum2[1] = sum2[1] + prod2[1];
335 }
336 sum1[0] = sum1[0] + sum2[0];
337 sum1[1] = sum1[1] + sum2[1];
338 {
339 tmp1[0] =
340 (double) sum1[0] * alpha_i[0] - (double) sum1[1] * alpha_i[1];
341 tmp1[1] =
342 (double) sum1[0] * alpha_i[1] + (double) sum1[1] * alpha_i[0];
343 }
344 y_elem[0] = y_i[yi];
345 y_elem[1] = y_i[yi + 1];
346 {
347 tmp2[0] =
348 (double) y_elem[0] * beta_i[0] - (double) y_elem[1] * beta_i[1];
349 tmp2[1] =
350 (double) y_elem[0] * beta_i[1] + (double) y_elem[1] * beta_i[0];
351 }
352 tmp3[0] = tmp1[0] + tmp2[0];
353 tmp3[1] = tmp1[1] + tmp2[1];
354 y_i[yi] = tmp3[0];
355 y_i[yi + 1] = tmp3[1];
356 }
357 }
358
359
360
361 break;
362 }
363
364 case blas_prec_double:
365 case blas_prec_indigenous:{
366
367 int i, j;
368 int xi, yi, xi0, yi0;
369 int aij, ai;
370 int incai;
371 int incaij, incaij2;
372
373 const float *a_i = (float *) a;
374 const float *x_head_i = (float *) x_head;
375 const float *x_tail_i = (float *) x_tail;
376 double *y_i = (double *) y;
377 double *alpha_i = (double *) alpha;
378 double *beta_i = (double *) beta;
379 float a_elem[2];
380 float x_elem[2];
381 double y_elem[2];
382 float diag_elem;
383 double prod1[2];
384 double prod2[2];
385 double sum1[2];
386 double sum2[2];
387 double tmp1[2];
388 double tmp2[2];
389 double tmp3[2];
390
391
392
393 /* Test for no-op */
394 if (n <= 0) {
395 return;
396 }
397 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0
398 && (beta_i[0] == 1.0 && beta_i[1] == 0.0)) {
399 return;
400 }
401
402 /* Check for error conditions. */
403 if (n < 0) {
404 BLAS_error(routine_name, -3, n, NULL);
405 }
406 if (lda < n) {
407 BLAS_error(routine_name, -6, n, NULL);
408 }
409 if (incx == 0) {
410 BLAS_error(routine_name, -9, incx, NULL);
411 }
412 if (incy == 0) {
413 BLAS_error(routine_name, -12, incy, NULL);
414 }
415
416 if ((order == blas_colmajor && uplo == blas_upper) ||
417 (order == blas_rowmajor && uplo == blas_lower)) {
418 incai = lda;
419 incaij = 1;
420 incaij2 = lda;
421 } else {
422 incai = 1;
423 incaij = lda;
424 incaij2 = 1;
425 }
426
427 incx *= 2;
428 incy *= 2;
429 incai *= 2;
430 incaij *= 2;
431 incaij2 *= 2;
432 xi0 = (incx > 0) ? 0 : ((-n + 1) * incx);
433 yi0 = (incy > 0) ? 0 : ((-n + 1) * incy);
434
435
436
437 /* The most general form, y <--- alpha * A * (x_head + x_tail) + beta * y */
438 if (uplo == blas_lower) {
439 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
440 sum1[0] = sum1[1] = 0.0;
441 sum2[0] = sum2[1] = 0.0;
442 for (j = 0, aij = ai, xi = xi0; j < i;
443 j++, aij += incaij, xi += incx) {
444 a_elem[0] = a_i[aij];
445 a_elem[1] = a_i[aij + 1];
446 x_elem[0] = x_head_i[xi];
447 x_elem[1] = x_head_i[xi + 1];
448 {
449 prod1[0] =
450 (double) a_elem[0] * x_elem[0] -
451 (double) a_elem[1] * x_elem[1];
452 prod1[1] =
453 (double) a_elem[0] * x_elem[1] +
454 (double) a_elem[1] * x_elem[0];
455 }
456 sum1[0] = sum1[0] + prod1[0];
457 sum1[1] = sum1[1] + prod1[1];
458 x_elem[0] = x_tail_i[xi];
459 x_elem[1] = x_tail_i[xi + 1];
460 {
461 prod2[0] =
462 (double) a_elem[0] * x_elem[0] -
463 (double) a_elem[1] * x_elem[1];
464 prod2[1] =
465 (double) a_elem[0] * x_elem[1] +
466 (double) a_elem[1] * x_elem[0];
467 }
468 sum2[0] = sum2[0] + prod2[0];
469 sum2[1] = sum2[1] + prod2[1];
470 }
471
472 diag_elem = a_i[aij];
473 x_elem[0] = x_head_i[xi];
474 x_elem[1] = x_head_i[xi + 1];
475 {
476 prod1[0] = (double) x_elem[0] * diag_elem;
477 prod1[1] = (double) x_elem[1] * diag_elem;
478 }
479 sum1[0] = sum1[0] + prod1[0];
480 sum1[1] = sum1[1] + prod1[1];
481 x_elem[0] = x_tail_i[xi];
482 x_elem[1] = x_tail_i[xi + 1];
483 {
484 prod2[0] = (double) x_elem[0] * diag_elem;
485 prod2[1] = (double) x_elem[1] * diag_elem;
486 }
487 sum2[0] = sum2[0] + prod2[0];
488 sum2[1] = sum2[1] + prod2[1];
489 j++;
490 aij += incaij2;
491 xi += incx;
492
493 for (; j < n; j++, aij += incaij2, xi += incx) {
494 a_elem[0] = a_i[aij];
495 a_elem[1] = a_i[aij + 1];
496 a_elem[1] = -a_elem[1];
497 x_elem[0] = x_head_i[xi];
498 x_elem[1] = x_head_i[xi + 1];
499 {
500 prod1[0] =
501 (double) a_elem[0] * x_elem[0] -
502 (double) a_elem[1] * x_elem[1];
503 prod1[1] =
504 (double) a_elem[0] * x_elem[1] +
505 (double) a_elem[1] * x_elem[0];
506 }
507 sum1[0] = sum1[0] + prod1[0];
508 sum1[1] = sum1[1] + prod1[1];
509 x_elem[0] = x_tail_i[xi];
510 x_elem[1] = x_tail_i[xi + 1];
511 {
512 prod2[0] =
513 (double) a_elem[0] * x_elem[0] -
514 (double) a_elem[1] * x_elem[1];
515 prod2[1] =
516 (double) a_elem[0] * x_elem[1] +
517 (double) a_elem[1] * x_elem[0];
518 }
519 sum2[0] = sum2[0] + prod2[0];
520 sum2[1] = sum2[1] + prod2[1];
521 }
522 sum1[0] = sum1[0] + sum2[0];
523 sum1[1] = sum1[1] + sum2[1];
524 {
525 tmp1[0] =
526 (double) sum1[0] * alpha_i[0] - (double) sum1[1] * alpha_i[1];
527 tmp1[1] =
528 (double) sum1[0] * alpha_i[1] + (double) sum1[1] * alpha_i[0];
529 }
530 y_elem[0] = y_i[yi];
531 y_elem[1] = y_i[yi + 1];
532 {
533 tmp2[0] =
534 (double) y_elem[0] * beta_i[0] - (double) y_elem[1] * beta_i[1];
535 tmp2[1] =
536 (double) y_elem[0] * beta_i[1] + (double) y_elem[1] * beta_i[0];
537 }
538 tmp3[0] = tmp1[0] + tmp2[0];
539 tmp3[1] = tmp1[1] + tmp2[1];
540 y_i[yi] = tmp3[0];
541 y_i[yi + 1] = tmp3[1];
542 }
543 } else {
544 /* uplo == blas_upper */
545 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
546 sum1[0] = sum1[1] = 0.0;
547 sum2[0] = sum2[1] = 0.0;
548
549 for (j = 0, aij = ai, xi = xi0; j < i;
550 j++, aij += incaij, xi += incx) {
551 a_elem[0] = a_i[aij];
552 a_elem[1] = a_i[aij + 1];
553 a_elem[1] = -a_elem[1];
554 x_elem[0] = x_head_i[xi];
555 x_elem[1] = x_head_i[xi + 1];
556 {
557 prod1[0] =
558 (double) a_elem[0] * x_elem[0] -
559 (double) a_elem[1] * x_elem[1];
560 prod1[1] =
561 (double) a_elem[0] * x_elem[1] +
562 (double) a_elem[1] * x_elem[0];
563 }
564 sum1[0] = sum1[0] + prod1[0];
565 sum1[1] = sum1[1] + prod1[1];
566 x_elem[0] = x_tail_i[xi];
567 x_elem[1] = x_tail_i[xi + 1];
568 {
569 prod2[0] =
570 (double) a_elem[0] * x_elem[0] -
571 (double) a_elem[1] * x_elem[1];
572 prod2[1] =
573 (double) a_elem[0] * x_elem[1] +
574 (double) a_elem[1] * x_elem[0];
575 }
576 sum2[0] = sum2[0] + prod2[0];
577 sum2[1] = sum2[1] + prod2[1];
578 }
579
580 diag_elem = a_i[aij];
581 x_elem[0] = x_head_i[xi];
582 x_elem[1] = x_head_i[xi + 1];
583 {
584 prod1[0] = (double) x_elem[0] * diag_elem;
585 prod1[1] = (double) x_elem[1] * diag_elem;
586 }
587 sum1[0] = sum1[0] + prod1[0];
588 sum1[1] = sum1[1] + prod1[1];
589 x_elem[0] = x_tail_i[xi];
590 x_elem[1] = x_tail_i[xi + 1];
591 {
592 prod2[0] = (double) x_elem[0] * diag_elem;
593 prod2[1] = (double) x_elem[1] * diag_elem;
594 }
595 sum2[0] = sum2[0] + prod2[0];
596 sum2[1] = sum2[1] + prod2[1];
597 j++;
598 aij += incaij2;
599 xi += incx;
600
601 for (; j < n; j++, aij += incaij2, xi += incx) {
602 a_elem[0] = a_i[aij];
603 a_elem[1] = a_i[aij + 1];
604 x_elem[0] = x_head_i[xi];
605 x_elem[1] = x_head_i[xi + 1];
606 {
607 prod1[0] =
608 (double) a_elem[0] * x_elem[0] -
609 (double) a_elem[1] * x_elem[1];
610 prod1[1] =
611 (double) a_elem[0] * x_elem[1] +
612 (double) a_elem[1] * x_elem[0];
613 }
614 sum1[0] = sum1[0] + prod1[0];
615 sum1[1] = sum1[1] + prod1[1];
616 x_elem[0] = x_tail_i[xi];
617 x_elem[1] = x_tail_i[xi + 1];
618 {
619 prod2[0] =
620 (double) a_elem[0] * x_elem[0] -
621 (double) a_elem[1] * x_elem[1];
622 prod2[1] =
623 (double) a_elem[0] * x_elem[1] +
624 (double) a_elem[1] * x_elem[0];
625 }
626 sum2[0] = sum2[0] + prod2[0];
627 sum2[1] = sum2[1] + prod2[1];
628 }
629 sum1[0] = sum1[0] + sum2[0];
630 sum1[1] = sum1[1] + sum2[1];
631 {
632 tmp1[0] =
633 (double) sum1[0] * alpha_i[0] - (double) sum1[1] * alpha_i[1];
634 tmp1[1] =
635 (double) sum1[0] * alpha_i[1] + (double) sum1[1] * alpha_i[0];
636 }
637 y_elem[0] = y_i[yi];
638 y_elem[1] = y_i[yi + 1];
639 {
640 tmp2[0] =
641 (double) y_elem[0] * beta_i[0] - (double) y_elem[1] * beta_i[1];
642 tmp2[1] =
643 (double) y_elem[0] * beta_i[1] + (double) y_elem[1] * beta_i[0];
644 }
645 tmp3[0] = tmp1[0] + tmp2[0];
646 tmp3[1] = tmp1[1] + tmp2[1];
647 y_i[yi] = tmp3[0];
648 y_i[yi + 1] = tmp3[1];
649 }
650 }
651
652
653
654 break;
655 }
656
657 case blas_prec_extra:{
658
659 int i, j;
660 int xi, yi, xi0, yi0;
661 int aij, ai;
662 int incai;
663 int incaij, incaij2;
664
665 const float *a_i = (float *) a;
666 const float *x_head_i = (float *) x_head;
667 const float *x_tail_i = (float *) x_tail;
668 double *y_i = (double *) y;
669 double *alpha_i = (double *) alpha;
670 double *beta_i = (double *) beta;
671 float a_elem[2];
672 float x_elem[2];
673 double y_elem[2];
674 float diag_elem;
675 double head_prod1[2], tail_prod1[2];
676 double head_prod2[2], tail_prod2[2];
677 double head_sum1[2], tail_sum1[2];
678 double head_sum2[2], tail_sum2[2];
679 double head_tmp1[2], tail_tmp1[2];
680 double head_tmp2[2], tail_tmp2[2];
681 double head_tmp3[2], tail_tmp3[2];
682
683 FPU_FIX_DECL;
684
685 /* Test for no-op */
686 if (n <= 0) {
687 return;
688 }
689 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0
690 && (beta_i[0] == 1.0 && beta_i[1] == 0.0)) {
691 return;
692 }
693
694 /* Check for error conditions. */
695 if (n < 0) {
696 BLAS_error(routine_name, -3, n, NULL);
697 }
698 if (lda < n) {
699 BLAS_error(routine_name, -6, n, NULL);
700 }
701 if (incx == 0) {
702 BLAS_error(routine_name, -9, incx, NULL);
703 }
704 if (incy == 0) {
705 BLAS_error(routine_name, -12, incy, NULL);
706 }
707
708 if ((order == blas_colmajor && uplo == blas_upper) ||
709 (order == blas_rowmajor && uplo == blas_lower)) {
710 incai = lda;
711 incaij = 1;
712 incaij2 = lda;
713 } else {
714 incai = 1;
715 incaij = lda;
716 incaij2 = 1;
717 }
718
719 incx *= 2;
720 incy *= 2;
721 incai *= 2;
722 incaij *= 2;
723 incaij2 *= 2;
724 xi0 = (incx > 0) ? 0 : ((-n + 1) * incx);
725 yi0 = (incy > 0) ? 0 : ((-n + 1) * incy);
726
727 FPU_FIX_START;
728
729 /* The most general form, y <--- alpha * A * (x_head + x_tail) + beta * y */
730 if (uplo == blas_lower) {
731 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
732 head_sum1[0] = head_sum1[1] = tail_sum1[0] = tail_sum1[1] = 0.0;
733 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] = 0.0;
734 for (j = 0, aij = ai, xi = xi0; j < i;
735 j++, aij += incaij, xi += incx) {
736 a_elem[0] = a_i[aij];
737 a_elem[1] = a_i[aij + 1];
738 x_elem[0] = x_head_i[xi];
739 x_elem[1] = x_head_i[xi + 1];
740 {
741 double head_e1, tail_e1;
742 double d1;
743 double d2;
744 /* Real part */
745 d1 = (double) a_elem[0] * x_elem[0];
746 d2 = (double) -a_elem[1] * x_elem[1];
747 {
748 /* Compute double-double = double + double. */
749 double e, t1, t2;
750
751 /* Knuth trick. */
752 t1 = d1 + d2;
753 e = t1 - d1;
754 t2 = ((d2 - e) + (d1 - (t1 - e)));
755
756 /* The result is t1 + t2, after normalization. */
757 head_e1 = t1 + t2;
758 tail_e1 = t2 - (head_e1 - t1);
759 }
760 head_prod1[0] = head_e1;
761 tail_prod1[0] = tail_e1;
762 /* imaginary part */
763 d1 = (double) a_elem[0] * x_elem[1];
764 d2 = (double) a_elem[1] * x_elem[0];
765 {
766 /* Compute double-double = double + double. */
767 double e, t1, t2;
768
769 /* Knuth trick. */
770 t1 = d1 + d2;
771 e = t1 - d1;
772 t2 = ((d2 - e) + (d1 - (t1 - e)));
773
774 /* The result is t1 + t2, after normalization. */
775 head_e1 = t1 + t2;
776 tail_e1 = t2 - (head_e1 - t1);
777 }
778 head_prod1[1] = head_e1;
779 tail_prod1[1] = tail_e1;
780 }
781 {
782 double head_t, tail_t;
783 double head_a, tail_a;
784 double head_b, tail_b;
785 /* Real part */
786 head_a = head_sum1[0];
787 tail_a = tail_sum1[0];
788 head_b = head_prod1[0];
789 tail_b = tail_prod1[0];
790 {
791 /* Compute double-double = double-double + double-double. */
792 double bv;
793 double s1, s2, t1, t2;
794
795 /* Add two hi words. */
796 s1 = head_a + head_b;
797 bv = s1 - head_a;
798 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
799
800 /* Add two lo words. */
801 t1 = tail_a + tail_b;
802 bv = t1 - tail_a;
803 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
804
805 s2 += t1;
806
807 /* Renormalize (s1, s2) to (t1, s2) */
808 t1 = s1 + s2;
809 s2 = s2 - (t1 - s1);
810
811 t2 += s2;
812
813 /* Renormalize (t1, t2) */
814 head_t = t1 + t2;
815 tail_t = t2 - (head_t - t1);
816 }
817 head_sum1[0] = head_t;
818 tail_sum1[0] = tail_t;
819 /* Imaginary part */
820 head_a = head_sum1[1];
821 tail_a = tail_sum1[1];
822 head_b = head_prod1[1];
823 tail_b = tail_prod1[1];
824 {
825 /* Compute double-double = double-double + double-double. */
826 double bv;
827 double s1, s2, t1, t2;
828
829 /* Add two hi words. */
830 s1 = head_a + head_b;
831 bv = s1 - head_a;
832 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
833
834 /* Add two lo words. */
835 t1 = tail_a + tail_b;
836 bv = t1 - tail_a;
837 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
838
839 s2 += t1;
840
841 /* Renormalize (s1, s2) to (t1, s2) */
842 t1 = s1 + s2;
843 s2 = s2 - (t1 - s1);
844
845 t2 += s2;
846
847 /* Renormalize (t1, t2) */
848 head_t = t1 + t2;
849 tail_t = t2 - (head_t - t1);
850 }
851 head_sum1[1] = head_t;
852 tail_sum1[1] = tail_t;
853 }
854 x_elem[0] = x_tail_i[xi];
855 x_elem[1] = x_tail_i[xi + 1];
856 {
857 double head_e1, tail_e1;
858 double d1;
859 double d2;
860 /* Real part */
861 d1 = (double) a_elem[0] * x_elem[0];
862 d2 = (double) -a_elem[1] * x_elem[1];
863 {
864 /* Compute double-double = double + double. */
865 double e, t1, t2;
866
867 /* Knuth trick. */
868 t1 = d1 + d2;
869 e = t1 - d1;
870 t2 = ((d2 - e) + (d1 - (t1 - e)));
871
872 /* The result is t1 + t2, after normalization. */
873 head_e1 = t1 + t2;
874 tail_e1 = t2 - (head_e1 - t1);
875 }
876 head_prod2[0] = head_e1;
877 tail_prod2[0] = tail_e1;
878 /* imaginary part */
879 d1 = (double) a_elem[0] * x_elem[1];
880 d2 = (double) a_elem[1] * x_elem[0];
881 {
882 /* Compute double-double = double + double. */
883 double e, t1, t2;
884
885 /* Knuth trick. */
886 t1 = d1 + d2;
887 e = t1 - d1;
888 t2 = ((d2 - e) + (d1 - (t1 - e)));
889
890 /* The result is t1 + t2, after normalization. */
891 head_e1 = t1 + t2;
892 tail_e1 = t2 - (head_e1 - t1);
893 }
894 head_prod2[1] = head_e1;
895 tail_prod2[1] = tail_e1;
896 }
897 {
898 double head_t, tail_t;
899 double head_a, tail_a;
900 double head_b, tail_b;
901 /* Real part */
902 head_a = head_sum2[0];
903 tail_a = tail_sum2[0];
904 head_b = head_prod2[0];
905 tail_b = tail_prod2[0];
906 {
907 /* Compute double-double = double-double + double-double. */
908 double bv;
909 double s1, s2, t1, t2;
910
911 /* Add two hi words. */
912 s1 = head_a + head_b;
913 bv = s1 - head_a;
914 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
915
916 /* Add two lo words. */
917 t1 = tail_a + tail_b;
918 bv = t1 - tail_a;
919 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
920
921 s2 += t1;
922
923 /* Renormalize (s1, s2) to (t1, s2) */
924 t1 = s1 + s2;
925 s2 = s2 - (t1 - s1);
926
927 t2 += s2;
928
929 /* Renormalize (t1, t2) */
930 head_t = t1 + t2;
931 tail_t = t2 - (head_t - t1);
932 }
933 head_sum2[0] = head_t;
934 tail_sum2[0] = tail_t;
935 /* Imaginary part */
936 head_a = head_sum2[1];
937 tail_a = tail_sum2[1];
938 head_b = head_prod2[1];
939 tail_b = tail_prod2[1];
940 {
941 /* Compute double-double = double-double + double-double. */
942 double bv;
943 double s1, s2, t1, t2;
944
945 /* Add two hi words. */
946 s1 = head_a + head_b;
947 bv = s1 - head_a;
948 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
949
950 /* Add two lo words. */
951 t1 = tail_a + tail_b;
952 bv = t1 - tail_a;
953 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
954
955 s2 += t1;
956
957 /* Renormalize (s1, s2) to (t1, s2) */
958 t1 = s1 + s2;
959 s2 = s2 - (t1 - s1);
960
961 t2 += s2;
962
963 /* Renormalize (t1, t2) */
964 head_t = t1 + t2;
965 tail_t = t2 - (head_t - t1);
966 }
967 head_sum2[1] = head_t;
968 tail_sum2[1] = tail_t;
969 }
970 }
971
972 diag_elem = a_i[aij];
973 x_elem[0] = x_head_i[xi];
974 x_elem[1] = x_head_i[xi + 1];
975 {
976 head_prod1[0] = (double) x_elem[0] * diag_elem;
977 tail_prod1[0] = 0.0;
978 head_prod1[1] = (double) x_elem[1] * diag_elem;
979 tail_prod1[1] = 0.0;
980 }
981 {
982 double head_t, tail_t;
983 double head_a, tail_a;
984 double head_b, tail_b;
985 /* Real part */
986 head_a = head_sum1[0];
987 tail_a = tail_sum1[0];
988 head_b = head_prod1[0];
989 tail_b = tail_prod1[0];
990 {
991 /* Compute double-double = double-double + double-double. */
992 double bv;
993 double s1, s2, t1, t2;
994
995 /* Add two hi words. */
996 s1 = head_a + head_b;
997 bv = s1 - head_a;
998 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
999
1000 /* Add two lo words. */
1001 t1 = tail_a + tail_b;
1002 bv = t1 - tail_a;
1003 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1004
1005 s2 += t1;
1006
1007 /* Renormalize (s1, s2) to (t1, s2) */
1008 t1 = s1 + s2;
1009 s2 = s2 - (t1 - s1);
1010
1011 t2 += s2;
1012
1013 /* Renormalize (t1, t2) */
1014 head_t = t1 + t2;
1015 tail_t = t2 - (head_t - t1);
1016 }
1017 head_sum1[0] = head_t;
1018 tail_sum1[0] = tail_t;
1019 /* Imaginary part */
1020 head_a = head_sum1[1];
1021 tail_a = tail_sum1[1];
1022 head_b = head_prod1[1];
1023 tail_b = tail_prod1[1];
1024 {
1025 /* Compute double-double = double-double + double-double. */
1026 double bv;
1027 double s1, s2, t1, t2;
1028
1029 /* Add two hi words. */
1030 s1 = head_a + head_b;
1031 bv = s1 - head_a;
1032 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1033
1034 /* Add two lo words. */
1035 t1 = tail_a + tail_b;
1036 bv = t1 - tail_a;
1037 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1038
1039 s2 += t1;
1040
1041 /* Renormalize (s1, s2) to (t1, s2) */
1042 t1 = s1 + s2;
1043 s2 = s2 - (t1 - s1);
1044
1045 t2 += s2;
1046
1047 /* Renormalize (t1, t2) */
1048 head_t = t1 + t2;
1049 tail_t = t2 - (head_t - t1);
1050 }
1051 head_sum1[1] = head_t;
1052 tail_sum1[1] = tail_t;
1053 }
1054 x_elem[0] = x_tail_i[xi];
1055 x_elem[1] = x_tail_i[xi + 1];
1056 {
1057 head_prod2[0] = (double) x_elem[0] * diag_elem;
1058 tail_prod2[0] = 0.0;
1059 head_prod2[1] = (double) x_elem[1] * diag_elem;
1060 tail_prod2[1] = 0.0;
1061 }
1062 {
1063 double head_t, tail_t;
1064 double head_a, tail_a;
1065 double head_b, tail_b;
1066 /* Real part */
1067 head_a = head_sum2[0];
1068 tail_a = tail_sum2[0];
1069 head_b = head_prod2[0];
1070 tail_b = tail_prod2[0];
1071 {
1072 /* Compute double-double = double-double + double-double. */
1073 double bv;
1074 double s1, s2, t1, t2;
1075
1076 /* Add two hi words. */
1077 s1 = head_a + head_b;
1078 bv = s1 - head_a;
1079 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1080
1081 /* Add two lo words. */
1082 t1 = tail_a + tail_b;
1083 bv = t1 - tail_a;
1084 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1085
1086 s2 += t1;
1087
1088 /* Renormalize (s1, s2) to (t1, s2) */
1089 t1 = s1 + s2;
1090 s2 = s2 - (t1 - s1);
1091
1092 t2 += s2;
1093
1094 /* Renormalize (t1, t2) */
1095 head_t = t1 + t2;
1096 tail_t = t2 - (head_t - t1);
1097 }
1098 head_sum2[0] = head_t;
1099 tail_sum2[0] = tail_t;
1100 /* Imaginary part */
1101 head_a = head_sum2[1];
1102 tail_a = tail_sum2[1];
1103 head_b = head_prod2[1];
1104 tail_b = tail_prod2[1];
1105 {
1106 /* Compute double-double = double-double + double-double. */
1107 double bv;
1108 double s1, s2, t1, t2;
1109
1110 /* Add two hi words. */
1111 s1 = head_a + head_b;
1112 bv = s1 - head_a;
1113 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1114
1115 /* Add two lo words. */
1116 t1 = tail_a + tail_b;
1117 bv = t1 - tail_a;
1118 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1119
1120 s2 += t1;
1121
1122 /* Renormalize (s1, s2) to (t1, s2) */
1123 t1 = s1 + s2;
1124 s2 = s2 - (t1 - s1);
1125
1126 t2 += s2;
1127
1128 /* Renormalize (t1, t2) */
1129 head_t = t1 + t2;
1130 tail_t = t2 - (head_t - t1);
1131 }
1132 head_sum2[1] = head_t;
1133 tail_sum2[1] = tail_t;
1134 }
1135 j++;
1136 aij += incaij2;
1137 xi += incx;
1138
1139 for (; j < n; j++, aij += incaij2, xi += incx) {
1140 a_elem[0] = a_i[aij];
1141 a_elem[1] = a_i[aij + 1];
1142 a_elem[1] = -a_elem[1];
1143 x_elem[0] = x_head_i[xi];
1144 x_elem[1] = x_head_i[xi + 1];
1145 {
1146 double head_e1, tail_e1;
1147 double d1;
1148 double d2;
1149 /* Real part */
1150 d1 = (double) a_elem[0] * x_elem[0];
1151 d2 = (double) -a_elem[1] * x_elem[1];
1152 {
1153 /* Compute double-double = double + double. */
1154 double e, t1, t2;
1155
1156 /* Knuth trick. */
1157 t1 = d1 + d2;
1158 e = t1 - d1;
1159 t2 = ((d2 - e) + (d1 - (t1 - e)));
1160
1161 /* The result is t1 + t2, after normalization. */
1162 head_e1 = t1 + t2;
1163 tail_e1 = t2 - (head_e1 - t1);
1164 }
1165 head_prod1[0] = head_e1;
1166 tail_prod1[0] = tail_e1;
1167 /* imaginary part */
1168 d1 = (double) a_elem[0] * x_elem[1];
1169 d2 = (double) a_elem[1] * x_elem[0];
1170 {
1171 /* Compute double-double = double + double. */
1172 double e, t1, t2;
1173
1174 /* Knuth trick. */
1175 t1 = d1 + d2;
1176 e = t1 - d1;
1177 t2 = ((d2 - e) + (d1 - (t1 - e)));
1178
1179 /* The result is t1 + t2, after normalization. */
1180 head_e1 = t1 + t2;
1181 tail_e1 = t2 - (head_e1 - t1);
1182 }
1183 head_prod1[1] = head_e1;
1184 tail_prod1[1] = tail_e1;
1185 }
1186 {
1187 double head_t, tail_t;
1188 double head_a, tail_a;
1189 double head_b, tail_b;
1190 /* Real part */
1191 head_a = head_sum1[0];
1192 tail_a = tail_sum1[0];
1193 head_b = head_prod1[0];
1194 tail_b = tail_prod1[0];
1195 {
1196 /* Compute double-double = double-double + double-double. */
1197 double bv;
1198 double s1, s2, t1, t2;
1199
1200 /* Add two hi words. */
1201 s1 = head_a + head_b;
1202 bv = s1 - head_a;
1203 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1204
1205 /* Add two lo words. */
1206 t1 = tail_a + tail_b;
1207 bv = t1 - tail_a;
1208 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1209
1210 s2 += t1;
1211
1212 /* Renormalize (s1, s2) to (t1, s2) */
1213 t1 = s1 + s2;
1214 s2 = s2 - (t1 - s1);
1215
1216 t2 += s2;
1217
1218 /* Renormalize (t1, t2) */
1219 head_t = t1 + t2;
1220 tail_t = t2 - (head_t - t1);
1221 }
1222 head_sum1[0] = head_t;
1223 tail_sum1[0] = tail_t;
1224 /* Imaginary part */
1225 head_a = head_sum1[1];
1226 tail_a = tail_sum1[1];
1227 head_b = head_prod1[1];
1228 tail_b = tail_prod1[1];
1229 {
1230 /* Compute double-double = double-double + double-double. */
1231 double bv;
1232 double s1, s2, t1, t2;
1233
1234 /* Add two hi words. */
1235 s1 = head_a + head_b;
1236 bv = s1 - head_a;
1237 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1238
1239 /* Add two lo words. */
1240 t1 = tail_a + tail_b;
1241 bv = t1 - tail_a;
1242 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1243
1244 s2 += t1;
1245
1246 /* Renormalize (s1, s2) to (t1, s2) */
1247 t1 = s1 + s2;
1248 s2 = s2 - (t1 - s1);
1249
1250 t2 += s2;
1251
1252 /* Renormalize (t1, t2) */
1253 head_t = t1 + t2;
1254 tail_t = t2 - (head_t - t1);
1255 }
1256 head_sum1[1] = head_t;
1257 tail_sum1[1] = tail_t;
1258 }
1259 x_elem[0] = x_tail_i[xi];
1260 x_elem[1] = x_tail_i[xi + 1];
1261 {
1262 double head_e1, tail_e1;
1263 double d1;
1264 double d2;
1265 /* Real part */
1266 d1 = (double) a_elem[0] * x_elem[0];
1267 d2 = (double) -a_elem[1] * x_elem[1];
1268 {
1269 /* Compute double-double = double + double. */
1270 double e, t1, t2;
1271
1272 /* Knuth trick. */
1273 t1 = d1 + d2;
1274 e = t1 - d1;
1275 t2 = ((d2 - e) + (d1 - (t1 - e)));
1276
1277 /* The result is t1 + t2, after normalization. */
1278 head_e1 = t1 + t2;
1279 tail_e1 = t2 - (head_e1 - t1);
1280 }
1281 head_prod2[0] = head_e1;
1282 tail_prod2[0] = tail_e1;
1283 /* imaginary part */
1284 d1 = (double) a_elem[0] * x_elem[1];
1285 d2 = (double) a_elem[1] * x_elem[0];
1286 {
1287 /* Compute double-double = double + double. */
1288 double e, t1, t2;
1289
1290 /* Knuth trick. */
1291 t1 = d1 + d2;
1292 e = t1 - d1;
1293 t2 = ((d2 - e) + (d1 - (t1 - e)));
1294
1295 /* The result is t1 + t2, after normalization. */
1296 head_e1 = t1 + t2;
1297 tail_e1 = t2 - (head_e1 - t1);
1298 }
1299 head_prod2[1] = head_e1;
1300 tail_prod2[1] = tail_e1;
1301 }
1302 {
1303 double head_t, tail_t;
1304 double head_a, tail_a;
1305 double head_b, tail_b;
1306 /* Real part */
1307 head_a = head_sum2[0];
1308 tail_a = tail_sum2[0];
1309 head_b = head_prod2[0];
1310 tail_b = tail_prod2[0];
1311 {
1312 /* Compute double-double = double-double + double-double. */
1313 double bv;
1314 double s1, s2, t1, t2;
1315
1316 /* Add two hi words. */
1317 s1 = head_a + head_b;
1318 bv = s1 - head_a;
1319 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1320
1321 /* Add two lo words. */
1322 t1 = tail_a + tail_b;
1323 bv = t1 - tail_a;
1324 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1325
1326 s2 += t1;
1327
1328 /* Renormalize (s1, s2) to (t1, s2) */
1329 t1 = s1 + s2;
1330 s2 = s2 - (t1 - s1);
1331
1332 t2 += s2;
1333
1334 /* Renormalize (t1, t2) */
1335 head_t = t1 + t2;
1336 tail_t = t2 - (head_t - t1);
1337 }
1338 head_sum2[0] = head_t;
1339 tail_sum2[0] = tail_t;
1340 /* Imaginary part */
1341 head_a = head_sum2[1];
1342 tail_a = tail_sum2[1];
1343 head_b = head_prod2[1];
1344 tail_b = tail_prod2[1];
1345 {
1346 /* Compute double-double = double-double + double-double. */
1347 double bv;
1348 double s1, s2, t1, t2;
1349
1350 /* Add two hi words. */
1351 s1 = head_a + head_b;
1352 bv = s1 - head_a;
1353 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1354
1355 /* Add two lo words. */
1356 t1 = tail_a + tail_b;
1357 bv = t1 - tail_a;
1358 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1359
1360 s2 += t1;
1361
1362 /* Renormalize (s1, s2) to (t1, s2) */
1363 t1 = s1 + s2;
1364 s2 = s2 - (t1 - s1);
1365
1366 t2 += s2;
1367
1368 /* Renormalize (t1, t2) */
1369 head_t = t1 + t2;
1370 tail_t = t2 - (head_t - t1);
1371 }
1372 head_sum2[1] = head_t;
1373 tail_sum2[1] = tail_t;
1374 }
1375 }
1376 {
1377 double head_t, tail_t;
1378 double head_a, tail_a;
1379 double head_b, tail_b;
1380 /* Real part */
1381 head_a = head_sum1[0];
1382 tail_a = tail_sum1[0];
1383 head_b = head_sum2[0];
1384 tail_b = tail_sum2[0];
1385 {
1386 /* Compute double-double = double-double + double-double. */
1387 double bv;
1388 double s1, s2, t1, t2;
1389
1390 /* Add two hi words. */
1391 s1 = head_a + head_b;
1392 bv = s1 - head_a;
1393 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1394
1395 /* Add two lo words. */
1396 t1 = tail_a + tail_b;
1397 bv = t1 - tail_a;
1398 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1399
1400 s2 += t1;
1401
1402 /* Renormalize (s1, s2) to (t1, s2) */
1403 t1 = s1 + s2;
1404 s2 = s2 - (t1 - s1);
1405
1406 t2 += s2;
1407
1408 /* Renormalize (t1, t2) */
1409 head_t = t1 + t2;
1410 tail_t = t2 - (head_t - t1);
1411 }
1412 head_sum1[0] = head_t;
1413 tail_sum1[0] = tail_t;
1414 /* Imaginary part */
1415 head_a = head_sum1[1];
1416 tail_a = tail_sum1[1];
1417 head_b = head_sum2[1];
1418 tail_b = tail_sum2[1];
1419 {
1420 /* Compute double-double = double-double + double-double. */
1421 double bv;
1422 double s1, s2, t1, t2;
1423
1424 /* Add two hi words. */
1425 s1 = head_a + head_b;
1426 bv = s1 - head_a;
1427 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1428
1429 /* Add two lo words. */
1430 t1 = tail_a + tail_b;
1431 bv = t1 - tail_a;
1432 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1433
1434 s2 += t1;
1435
1436 /* Renormalize (s1, s2) to (t1, s2) */
1437 t1 = s1 + s2;
1438 s2 = s2 - (t1 - s1);
1439
1440 t2 += s2;
1441
1442 /* Renormalize (t1, t2) */
1443 head_t = t1 + t2;
1444 tail_t = t2 - (head_t - t1);
1445 }
1446 head_sum1[1] = head_t;
1447 tail_sum1[1] = tail_t;
1448 }
1449 {
1450 /* Compute complex-extra = complex-extra * complex-double. */
1451 double head_a0, tail_a0;
1452 double head_a1, tail_a1;
1453 double head_t1, tail_t1;
1454 double head_t2, tail_t2;
1455 head_a0 = head_sum1[0];
1456 tail_a0 = tail_sum1[0];
1457 head_a1 = head_sum1[1];
1458 tail_a1 = tail_sum1[1];
1459 /* real part */
1460 {
1461 /* Compute double-double = double-double * double. */
1462 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1463
1464 con = head_a0 * split;
1465 a11 = con - head_a0;
1466 a11 = con - a11;
1467 a21 = head_a0 - a11;
1468 con = alpha_i[0] * split;
1469 b1 = con - alpha_i[0];
1470 b1 = con - b1;
1471 b2 = alpha_i[0] - b1;
1472
1473 c11 = head_a0 * alpha_i[0];
1474 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1475
1476 c2 = tail_a0 * alpha_i[0];
1477 t1 = c11 + c2;
1478 t2 = (c2 - (t1 - c11)) + c21;
1479
1480 head_t1 = t1 + t2;
1481 tail_t1 = t2 - (head_t1 - t1);
1482 }
1483 {
1484 /* Compute double-double = double-double * double. */
1485 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1486
1487 con = head_a1 * split;
1488 a11 = con - head_a1;
1489 a11 = con - a11;
1490 a21 = head_a1 - a11;
1491 con = alpha_i[1] * split;
1492 b1 = con - alpha_i[1];
1493 b1 = con - b1;
1494 b2 = alpha_i[1] - b1;
1495
1496 c11 = head_a1 * alpha_i[1];
1497 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1498
1499 c2 = tail_a1 * alpha_i[1];
1500 t1 = c11 + c2;
1501 t2 = (c2 - (t1 - c11)) + c21;
1502
1503 head_t2 = t1 + t2;
1504 tail_t2 = t2 - (head_t2 - t1);
1505 }
1506 head_t2 = -head_t2;
1507 tail_t2 = -tail_t2;
1508 {
1509 /* Compute double-double = double-double + double-double. */
1510 double bv;
1511 double s1, s2, t1, t2;
1512
1513 /* Add two hi words. */
1514 s1 = head_t1 + head_t2;
1515 bv = s1 - head_t1;
1516 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1517
1518 /* Add two lo words. */
1519 t1 = tail_t1 + tail_t2;
1520 bv = t1 - tail_t1;
1521 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1522
1523 s2 += t1;
1524
1525 /* Renormalize (s1, s2) to (t1, s2) */
1526 t1 = s1 + s2;
1527 s2 = s2 - (t1 - s1);
1528
1529 t2 += s2;
1530
1531 /* Renormalize (t1, t2) */
1532 head_t1 = t1 + t2;
1533 tail_t1 = t2 - (head_t1 - t1);
1534 }
1535 head_tmp1[0] = head_t1;
1536 tail_tmp1[0] = tail_t1;
1537 /* imaginary part */
1538 {
1539 /* Compute double-double = double-double * double. */
1540 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1541
1542 con = head_a1 * split;
1543 a11 = con - head_a1;
1544 a11 = con - a11;
1545 a21 = head_a1 - a11;
1546 con = alpha_i[0] * split;
1547 b1 = con - alpha_i[0];
1548 b1 = con - b1;
1549 b2 = alpha_i[0] - b1;
1550
1551 c11 = head_a1 * alpha_i[0];
1552 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1553
1554 c2 = tail_a1 * alpha_i[0];
1555 t1 = c11 + c2;
1556 t2 = (c2 - (t1 - c11)) + c21;
1557
1558 head_t1 = t1 + t2;
1559 tail_t1 = t2 - (head_t1 - t1);
1560 }
1561 {
1562 /* Compute double-double = double-double * double. */
1563 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1564
1565 con = head_a0 * split;
1566 a11 = con - head_a0;
1567 a11 = con - a11;
1568 a21 = head_a0 - a11;
1569 con = alpha_i[1] * split;
1570 b1 = con - alpha_i[1];
1571 b1 = con - b1;
1572 b2 = alpha_i[1] - b1;
1573
1574 c11 = head_a0 * alpha_i[1];
1575 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1576
1577 c2 = tail_a0 * alpha_i[1];
1578 t1 = c11 + c2;
1579 t2 = (c2 - (t1 - c11)) + c21;
1580
1581 head_t2 = t1 + t2;
1582 tail_t2 = t2 - (head_t2 - t1);
1583 }
1584 {
1585 /* Compute double-double = double-double + double-double. */
1586 double bv;
1587 double s1, s2, t1, t2;
1588
1589 /* Add two hi words. */
1590 s1 = head_t1 + head_t2;
1591 bv = s1 - head_t1;
1592 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1593
1594 /* Add two lo words. */
1595 t1 = tail_t1 + tail_t2;
1596 bv = t1 - tail_t1;
1597 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1598
1599 s2 += t1;
1600
1601 /* Renormalize (s1, s2) to (t1, s2) */
1602 t1 = s1 + s2;
1603 s2 = s2 - (t1 - s1);
1604
1605 t2 += s2;
1606
1607 /* Renormalize (t1, t2) */
1608 head_t1 = t1 + t2;
1609 tail_t1 = t2 - (head_t1 - t1);
1610 }
1611 head_tmp1[1] = head_t1;
1612 tail_tmp1[1] = tail_t1;
1613 }
1614
1615 y_elem[0] = y_i[yi];
1616 y_elem[1] = y_i[yi + 1];
1617 {
1618 /* Compute complex-extra = complex-double * complex-double. */
1619 double head_t1, tail_t1;
1620 double head_t2, tail_t2;
1621 /* Real part */
1622 {
1623 /* Compute double_double = double * double. */
1624 double a1, a2, b1, b2, con;
1625
1626 con = y_elem[0] * split;
1627 a1 = con - y_elem[0];
1628 a1 = con - a1;
1629 a2 = y_elem[0] - a1;
1630 con = beta_i[0] * split;
1631 b1 = con - beta_i[0];
1632 b1 = con - b1;
1633 b2 = beta_i[0] - b1;
1634
1635 head_t1 = y_elem[0] * beta_i[0];
1636 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1637 }
1638 {
1639 /* Compute double_double = double * double. */
1640 double a1, a2, b1, b2, con;
1641
1642 con = y_elem[1] * split;
1643 a1 = con - y_elem[1];
1644 a1 = con - a1;
1645 a2 = y_elem[1] - a1;
1646 con = beta_i[1] * split;
1647 b1 = con - beta_i[1];
1648 b1 = con - b1;
1649 b2 = beta_i[1] - b1;
1650
1651 head_t2 = y_elem[1] * beta_i[1];
1652 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1653 }
1654 head_t2 = -head_t2;
1655 tail_t2 = -tail_t2;
1656 {
1657 /* Compute double-double = double-double + double-double. */
1658 double bv;
1659 double s1, s2, t1, t2;
1660
1661 /* Add two hi words. */
1662 s1 = head_t1 + head_t2;
1663 bv = s1 - head_t1;
1664 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1665
1666 /* Add two lo words. */
1667 t1 = tail_t1 + tail_t2;
1668 bv = t1 - tail_t1;
1669 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1670
1671 s2 += t1;
1672
1673 /* Renormalize (s1, s2) to (t1, s2) */
1674 t1 = s1 + s2;
1675 s2 = s2 - (t1 - s1);
1676
1677 t2 += s2;
1678
1679 /* Renormalize (t1, t2) */
1680 head_t1 = t1 + t2;
1681 tail_t1 = t2 - (head_t1 - t1);
1682 }
1683 head_tmp2[0] = head_t1;
1684 tail_tmp2[0] = tail_t1;
1685 /* Imaginary part */
1686 {
1687 /* Compute double_double = double * double. */
1688 double a1, a2, b1, b2, con;
1689
1690 con = y_elem[1] * split;
1691 a1 = con - y_elem[1];
1692 a1 = con - a1;
1693 a2 = y_elem[1] - a1;
1694 con = beta_i[0] * split;
1695 b1 = con - beta_i[0];
1696 b1 = con - b1;
1697 b2 = beta_i[0] - b1;
1698
1699 head_t1 = y_elem[1] * beta_i[0];
1700 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1701 }
1702 {
1703 /* Compute double_double = double * double. */
1704 double a1, a2, b1, b2, con;
1705
1706 con = y_elem[0] * split;
1707 a1 = con - y_elem[0];
1708 a1 = con - a1;
1709 a2 = y_elem[0] - a1;
1710 con = beta_i[1] * split;
1711 b1 = con - beta_i[1];
1712 b1 = con - b1;
1713 b2 = beta_i[1] - b1;
1714
1715 head_t2 = y_elem[0] * beta_i[1];
1716 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1717 }
1718 {
1719 /* Compute double-double = double-double + double-double. */
1720 double bv;
1721 double s1, s2, t1, t2;
1722
1723 /* Add two hi words. */
1724 s1 = head_t1 + head_t2;
1725 bv = s1 - head_t1;
1726 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1727
1728 /* Add two lo words. */
1729 t1 = tail_t1 + tail_t2;
1730 bv = t1 - tail_t1;
1731 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1732
1733 s2 += t1;
1734
1735 /* Renormalize (s1, s2) to (t1, s2) */
1736 t1 = s1 + s2;
1737 s2 = s2 - (t1 - s1);
1738
1739 t2 += s2;
1740
1741 /* Renormalize (t1, t2) */
1742 head_t1 = t1 + t2;
1743 tail_t1 = t2 - (head_t1 - t1);
1744 }
1745 head_tmp2[1] = head_t1;
1746 tail_tmp2[1] = tail_t1;
1747 }
1748 {
1749 double head_t, tail_t;
1750 double head_a, tail_a;
1751 double head_b, tail_b;
1752 /* Real part */
1753 head_a = head_tmp1[0];
1754 tail_a = tail_tmp1[0];
1755 head_b = head_tmp2[0];
1756 tail_b = tail_tmp2[0];
1757 {
1758 /* Compute double-double = double-double + double-double. */
1759 double bv;
1760 double s1, s2, t1, t2;
1761
1762 /* Add two hi words. */
1763 s1 = head_a + head_b;
1764 bv = s1 - head_a;
1765 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1766
1767 /* Add two lo words. */
1768 t1 = tail_a + tail_b;
1769 bv = t1 - tail_a;
1770 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1771
1772 s2 += t1;
1773
1774 /* Renormalize (s1, s2) to (t1, s2) */
1775 t1 = s1 + s2;
1776 s2 = s2 - (t1 - s1);
1777
1778 t2 += s2;
1779
1780 /* Renormalize (t1, t2) */
1781 head_t = t1 + t2;
1782 tail_t = t2 - (head_t - t1);
1783 }
1784 head_tmp3[0] = head_t;
1785 tail_tmp3[0] = tail_t;
1786 /* Imaginary part */
1787 head_a = head_tmp1[1];
1788 tail_a = tail_tmp1[1];
1789 head_b = head_tmp2[1];
1790 tail_b = tail_tmp2[1];
1791 {
1792 /* Compute double-double = double-double + double-double. */
1793 double bv;
1794 double s1, s2, t1, t2;
1795
1796 /* Add two hi words. */
1797 s1 = head_a + head_b;
1798 bv = s1 - head_a;
1799 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1800
1801 /* Add two lo words. */
1802 t1 = tail_a + tail_b;
1803 bv = t1 - tail_a;
1804 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1805
1806 s2 += t1;
1807
1808 /* Renormalize (s1, s2) to (t1, s2) */
1809 t1 = s1 + s2;
1810 s2 = s2 - (t1 - s1);
1811
1812 t2 += s2;
1813
1814 /* Renormalize (t1, t2) */
1815 head_t = t1 + t2;
1816 tail_t = t2 - (head_t - t1);
1817 }
1818 head_tmp3[1] = head_t;
1819 tail_tmp3[1] = tail_t;
1820 }
1821 y_i[yi] = head_tmp3[0];
1822 y_i[yi + 1] = head_tmp3[1];
1823 }
1824 } else {
1825 /* uplo == blas_upper */
1826 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
1827 head_sum1[0] = head_sum1[1] = tail_sum1[0] = tail_sum1[1] = 0.0;
1828 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] = 0.0;
1829
1830 for (j = 0, aij = ai, xi = xi0; j < i;
1831 j++, aij += incaij, xi += incx) {
1832 a_elem[0] = a_i[aij];
1833 a_elem[1] = a_i[aij + 1];
1834 a_elem[1] = -a_elem[1];
1835 x_elem[0] = x_head_i[xi];
1836 x_elem[1] = x_head_i[xi + 1];
1837 {
1838 double head_e1, tail_e1;
1839 double d1;
1840 double d2;
1841 /* Real part */
1842 d1 = (double) a_elem[0] * x_elem[0];
1843 d2 = (double) -a_elem[1] * x_elem[1];
1844 {
1845 /* Compute double-double = double + double. */
1846 double e, t1, t2;
1847
1848 /* Knuth trick. */
1849 t1 = d1 + d2;
1850 e = t1 - d1;
1851 t2 = ((d2 - e) + (d1 - (t1 - e)));
1852
1853 /* The result is t1 + t2, after normalization. */
1854 head_e1 = t1 + t2;
1855 tail_e1 = t2 - (head_e1 - t1);
1856 }
1857 head_prod1[0] = head_e1;
1858 tail_prod1[0] = tail_e1;
1859 /* imaginary part */
1860 d1 = (double) a_elem[0] * x_elem[1];
1861 d2 = (double) a_elem[1] * x_elem[0];
1862 {
1863 /* Compute double-double = double + double. */
1864 double e, t1, t2;
1865
1866 /* Knuth trick. */
1867 t1 = d1 + d2;
1868 e = t1 - d1;
1869 t2 = ((d2 - e) + (d1 - (t1 - e)));
1870
1871 /* The result is t1 + t2, after normalization. */
1872 head_e1 = t1 + t2;
1873 tail_e1 = t2 - (head_e1 - t1);
1874 }
1875 head_prod1[1] = head_e1;
1876 tail_prod1[1] = tail_e1;
1877 }
1878 {
1879 double head_t, tail_t;
1880 double head_a, tail_a;
1881 double head_b, tail_b;
1882 /* Real part */
1883 head_a = head_sum1[0];
1884 tail_a = tail_sum1[0];
1885 head_b = head_prod1[0];
1886 tail_b = tail_prod1[0];
1887 {
1888 /* Compute double-double = double-double + double-double. */
1889 double bv;
1890 double s1, s2, t1, t2;
1891
1892 /* Add two hi words. */
1893 s1 = head_a + head_b;
1894 bv = s1 - head_a;
1895 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1896
1897 /* Add two lo words. */
1898 t1 = tail_a + tail_b;
1899 bv = t1 - tail_a;
1900 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1901
1902 s2 += t1;
1903
1904 /* Renormalize (s1, s2) to (t1, s2) */
1905 t1 = s1 + s2;
1906 s2 = s2 - (t1 - s1);
1907
1908 t2 += s2;
1909
1910 /* Renormalize (t1, t2) */
1911 head_t = t1 + t2;
1912 tail_t = t2 - (head_t - t1);
1913 }
1914 head_sum1[0] = head_t;
1915 tail_sum1[0] = tail_t;
1916 /* Imaginary part */
1917 head_a = head_sum1[1];
1918 tail_a = tail_sum1[1];
1919 head_b = head_prod1[1];
1920 tail_b = tail_prod1[1];
1921 {
1922 /* Compute double-double = double-double + double-double. */
1923 double bv;
1924 double s1, s2, t1, t2;
1925
1926 /* Add two hi words. */
1927 s1 = head_a + head_b;
1928 bv = s1 - head_a;
1929 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1930
1931 /* Add two lo words. */
1932 t1 = tail_a + tail_b;
1933 bv = t1 - tail_a;
1934 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1935
1936 s2 += t1;
1937
1938 /* Renormalize (s1, s2) to (t1, s2) */
1939 t1 = s1 + s2;
1940 s2 = s2 - (t1 - s1);
1941
1942 t2 += s2;
1943
1944 /* Renormalize (t1, t2) */
1945 head_t = t1 + t2;
1946 tail_t = t2 - (head_t - t1);
1947 }
1948 head_sum1[1] = head_t;
1949 tail_sum1[1] = tail_t;
1950 }
1951 x_elem[0] = x_tail_i[xi];
1952 x_elem[1] = x_tail_i[xi + 1];
1953 {
1954 double head_e1, tail_e1;
1955 double d1;
1956 double d2;
1957 /* Real part */
1958 d1 = (double) a_elem[0] * x_elem[0];
1959 d2 = (double) -a_elem[1] * x_elem[1];
1960 {
1961 /* Compute double-double = double + double. */
1962 double e, t1, t2;
1963
1964 /* Knuth trick. */
1965 t1 = d1 + d2;
1966 e = t1 - d1;
1967 t2 = ((d2 - e) + (d1 - (t1 - e)));
1968
1969 /* The result is t1 + t2, after normalization. */
1970 head_e1 = t1 + t2;
1971 tail_e1 = t2 - (head_e1 - t1);
1972 }
1973 head_prod2[0] = head_e1;
1974 tail_prod2[0] = tail_e1;
1975 /* imaginary part */
1976 d1 = (double) a_elem[0] * x_elem[1];
1977 d2 = (double) a_elem[1] * x_elem[0];
1978 {
1979 /* Compute double-double = double + double. */
1980 double e, t1, t2;
1981
1982 /* Knuth trick. */
1983 t1 = d1 + d2;
1984 e = t1 - d1;
1985 t2 = ((d2 - e) + (d1 - (t1 - e)));
1986
1987 /* The result is t1 + t2, after normalization. */
1988 head_e1 = t1 + t2;
1989 tail_e1 = t2 - (head_e1 - t1);
1990 }
1991 head_prod2[1] = head_e1;
1992 tail_prod2[1] = tail_e1;
1993 }
1994 {
1995 double head_t, tail_t;
1996 double head_a, tail_a;
1997 double head_b, tail_b;
1998 /* Real part */
1999 head_a = head_sum2[0];
2000 tail_a = tail_sum2[0];
2001 head_b = head_prod2[0];
2002 tail_b = tail_prod2[0];
2003 {
2004 /* Compute double-double = double-double + double-double. */
2005 double bv;
2006 double s1, s2, t1, t2;
2007
2008 /* Add two hi words. */
2009 s1 = head_a + head_b;
2010 bv = s1 - head_a;
2011 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2012
2013 /* Add two lo words. */
2014 t1 = tail_a + tail_b;
2015 bv = t1 - tail_a;
2016 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2017
2018 s2 += t1;
2019
2020 /* Renormalize (s1, s2) to (t1, s2) */
2021 t1 = s1 + s2;
2022 s2 = s2 - (t1 - s1);
2023
2024 t2 += s2;
2025
2026 /* Renormalize (t1, t2) */
2027 head_t = t1 + t2;
2028 tail_t = t2 - (head_t - t1);
2029 }
2030 head_sum2[0] = head_t;
2031 tail_sum2[0] = tail_t;
2032 /* Imaginary part */
2033 head_a = head_sum2[1];
2034 tail_a = tail_sum2[1];
2035 head_b = head_prod2[1];
2036 tail_b = tail_prod2[1];
2037 {
2038 /* Compute double-double = double-double + double-double. */
2039 double bv;
2040 double s1, s2, t1, t2;
2041
2042 /* Add two hi words. */
2043 s1 = head_a + head_b;
2044 bv = s1 - head_a;
2045 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2046
2047 /* Add two lo words. */
2048 t1 = tail_a + tail_b;
2049 bv = t1 - tail_a;
2050 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2051
2052 s2 += t1;
2053
2054 /* Renormalize (s1, s2) to (t1, s2) */
2055 t1 = s1 + s2;
2056 s2 = s2 - (t1 - s1);
2057
2058 t2 += s2;
2059
2060 /* Renormalize (t1, t2) */
2061 head_t = t1 + t2;
2062 tail_t = t2 - (head_t - t1);
2063 }
2064 head_sum2[1] = head_t;
2065 tail_sum2[1] = tail_t;
2066 }
2067 }
2068
2069 diag_elem = a_i[aij];
2070 x_elem[0] = x_head_i[xi];
2071 x_elem[1] = x_head_i[xi + 1];
2072 {
2073 head_prod1[0] = (double) x_elem[0] * diag_elem;
2074 tail_prod1[0] = 0.0;
2075 head_prod1[1] = (double) x_elem[1] * diag_elem;
2076 tail_prod1[1] = 0.0;
2077 }
2078 {
2079 double head_t, tail_t;
2080 double head_a, tail_a;
2081 double head_b, tail_b;
2082 /* Real part */
2083 head_a = head_sum1[0];
2084 tail_a = tail_sum1[0];
2085 head_b = head_prod1[0];
2086 tail_b = tail_prod1[0];
2087 {
2088 /* Compute double-double = double-double + double-double. */
2089 double bv;
2090 double s1, s2, t1, t2;
2091
2092 /* Add two hi words. */
2093 s1 = head_a + head_b;
2094 bv = s1 - head_a;
2095 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2096
2097 /* Add two lo words. */
2098 t1 = tail_a + tail_b;
2099 bv = t1 - tail_a;
2100 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2101
2102 s2 += t1;
2103
2104 /* Renormalize (s1, s2) to (t1, s2) */
2105 t1 = s1 + s2;
2106 s2 = s2 - (t1 - s1);
2107
2108 t2 += s2;
2109
2110 /* Renormalize (t1, t2) */
2111 head_t = t1 + t2;
2112 tail_t = t2 - (head_t - t1);
2113 }
2114 head_sum1[0] = head_t;
2115 tail_sum1[0] = tail_t;
2116 /* Imaginary part */
2117 head_a = head_sum1[1];
2118 tail_a = tail_sum1[1];
2119 head_b = head_prod1[1];
2120 tail_b = tail_prod1[1];
2121 {
2122 /* Compute double-double = double-double + double-double. */
2123 double bv;
2124 double s1, s2, t1, t2;
2125
2126 /* Add two hi words. */
2127 s1 = head_a + head_b;
2128 bv = s1 - head_a;
2129 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2130
2131 /* Add two lo words. */
2132 t1 = tail_a + tail_b;
2133 bv = t1 - tail_a;
2134 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2135
2136 s2 += t1;
2137
2138 /* Renormalize (s1, s2) to (t1, s2) */
2139 t1 = s1 + s2;
2140 s2 = s2 - (t1 - s1);
2141
2142 t2 += s2;
2143
2144 /* Renormalize (t1, t2) */
2145 head_t = t1 + t2;
2146 tail_t = t2 - (head_t - t1);
2147 }
2148 head_sum1[1] = head_t;
2149 tail_sum1[1] = tail_t;
2150 }
2151 x_elem[0] = x_tail_i[xi];
2152 x_elem[1] = x_tail_i[xi + 1];
2153 {
2154 head_prod2[0] = (double) x_elem[0] * diag_elem;
2155 tail_prod2[0] = 0.0;
2156 head_prod2[1] = (double) x_elem[1] * diag_elem;
2157 tail_prod2[1] = 0.0;
2158 }
2159 {
2160 double head_t, tail_t;
2161 double head_a, tail_a;
2162 double head_b, tail_b;
2163 /* Real part */
2164 head_a = head_sum2[0];
2165 tail_a = tail_sum2[0];
2166 head_b = head_prod2[0];
2167 tail_b = tail_prod2[0];
2168 {
2169 /* Compute double-double = double-double + double-double. */
2170 double bv;
2171 double s1, s2, t1, t2;
2172
2173 /* Add two hi words. */
2174 s1 = head_a + head_b;
2175 bv = s1 - head_a;
2176 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2177
2178 /* Add two lo words. */
2179 t1 = tail_a + tail_b;
2180 bv = t1 - tail_a;
2181 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2182
2183 s2 += t1;
2184
2185 /* Renormalize (s1, s2) to (t1, s2) */
2186 t1 = s1 + s2;
2187 s2 = s2 - (t1 - s1);
2188
2189 t2 += s2;
2190
2191 /* Renormalize (t1, t2) */
2192 head_t = t1 + t2;
2193 tail_t = t2 - (head_t - t1);
2194 }
2195 head_sum2[0] = head_t;
2196 tail_sum2[0] = tail_t;
2197 /* Imaginary part */
2198 head_a = head_sum2[1];
2199 tail_a = tail_sum2[1];
2200 head_b = head_prod2[1];
2201 tail_b = tail_prod2[1];
2202 {
2203 /* Compute double-double = double-double + double-double. */
2204 double bv;
2205 double s1, s2, t1, t2;
2206
2207 /* Add two hi words. */
2208 s1 = head_a + head_b;
2209 bv = s1 - head_a;
2210 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2211
2212 /* Add two lo words. */
2213 t1 = tail_a + tail_b;
2214 bv = t1 - tail_a;
2215 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2216
2217 s2 += t1;
2218
2219 /* Renormalize (s1, s2) to (t1, s2) */
2220 t1 = s1 + s2;
2221 s2 = s2 - (t1 - s1);
2222
2223 t2 += s2;
2224
2225 /* Renormalize (t1, t2) */
2226 head_t = t1 + t2;
2227 tail_t = t2 - (head_t - t1);
2228 }
2229 head_sum2[1] = head_t;
2230 tail_sum2[1] = tail_t;
2231 }
2232 j++;
2233 aij += incaij2;
2234 xi += incx;
2235
2236 for (; j < n; j++, aij += incaij2, xi += incx) {
2237 a_elem[0] = a_i[aij];
2238 a_elem[1] = a_i[aij + 1];
2239 x_elem[0] = x_head_i[xi];
2240 x_elem[1] = x_head_i[xi + 1];
2241 {
2242 double head_e1, tail_e1;
2243 double d1;
2244 double d2;
2245 /* Real part */
2246 d1 = (double) a_elem[0] * x_elem[0];
2247 d2 = (double) -a_elem[1] * x_elem[1];
2248 {
2249 /* Compute double-double = double + double. */
2250 double e, t1, t2;
2251
2252 /* Knuth trick. */
2253 t1 = d1 + d2;
2254 e = t1 - d1;
2255 t2 = ((d2 - e) + (d1 - (t1 - e)));
2256
2257 /* The result is t1 + t2, after normalization. */
2258 head_e1 = t1 + t2;
2259 tail_e1 = t2 - (head_e1 - t1);
2260 }
2261 head_prod1[0] = head_e1;
2262 tail_prod1[0] = tail_e1;
2263 /* imaginary part */
2264 d1 = (double) a_elem[0] * x_elem[1];
2265 d2 = (double) a_elem[1] * x_elem[0];
2266 {
2267 /* Compute double-double = double + double. */
2268 double e, t1, t2;
2269
2270 /* Knuth trick. */
2271 t1 = d1 + d2;
2272 e = t1 - d1;
2273 t2 = ((d2 - e) + (d1 - (t1 - e)));
2274
2275 /* The result is t1 + t2, after normalization. */
2276 head_e1 = t1 + t2;
2277 tail_e1 = t2 - (head_e1 - t1);
2278 }
2279 head_prod1[1] = head_e1;
2280 tail_prod1[1] = tail_e1;
2281 }
2282 {
2283 double head_t, tail_t;
2284 double head_a, tail_a;
2285 double head_b, tail_b;
2286 /* Real part */
2287 head_a = head_sum1[0];
2288 tail_a = tail_sum1[0];
2289 head_b = head_prod1[0];
2290 tail_b = tail_prod1[0];
2291 {
2292 /* Compute double-double = double-double + double-double. */
2293 double bv;
2294 double s1, s2, t1, t2;
2295
2296 /* Add two hi words. */
2297 s1 = head_a + head_b;
2298 bv = s1 - head_a;
2299 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2300
2301 /* Add two lo words. */
2302 t1 = tail_a + tail_b;
2303 bv = t1 - tail_a;
2304 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2305
2306 s2 += t1;
2307
2308 /* Renormalize (s1, s2) to (t1, s2) */
2309 t1 = s1 + s2;
2310 s2 = s2 - (t1 - s1);
2311
2312 t2 += s2;
2313
2314 /* Renormalize (t1, t2) */
2315 head_t = t1 + t2;
2316 tail_t = t2 - (head_t - t1);
2317 }
2318 head_sum1[0] = head_t;
2319 tail_sum1[0] = tail_t;
2320 /* Imaginary part */
2321 head_a = head_sum1[1];
2322 tail_a = tail_sum1[1];
2323 head_b = head_prod1[1];
2324 tail_b = tail_prod1[1];
2325 {
2326 /* Compute double-double = double-double + double-double. */
2327 double bv;
2328 double s1, s2, t1, t2;
2329
2330 /* Add two hi words. */
2331 s1 = head_a + head_b;
2332 bv = s1 - head_a;
2333 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2334
2335 /* Add two lo words. */
2336 t1 = tail_a + tail_b;
2337 bv = t1 - tail_a;
2338 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2339
2340 s2 += t1;
2341
2342 /* Renormalize (s1, s2) to (t1, s2) */
2343 t1 = s1 + s2;
2344 s2 = s2 - (t1 - s1);
2345
2346 t2 += s2;
2347
2348 /* Renormalize (t1, t2) */
2349 head_t = t1 + t2;
2350 tail_t = t2 - (head_t - t1);
2351 }
2352 head_sum1[1] = head_t;
2353 tail_sum1[1] = tail_t;
2354 }
2355 x_elem[0] = x_tail_i[xi];
2356 x_elem[1] = x_tail_i[xi + 1];
2357 {
2358 double head_e1, tail_e1;
2359 double d1;
2360 double d2;
2361 /* Real part */
2362 d1 = (double) a_elem[0] * x_elem[0];
2363 d2 = (double) -a_elem[1] * x_elem[1];
2364 {
2365 /* Compute double-double = double + double. */
2366 double e, t1, t2;
2367
2368 /* Knuth trick. */
2369 t1 = d1 + d2;
2370 e = t1 - d1;
2371 t2 = ((d2 - e) + (d1 - (t1 - e)));
2372
2373 /* The result is t1 + t2, after normalization. */
2374 head_e1 = t1 + t2;
2375 tail_e1 = t2 - (head_e1 - t1);
2376 }
2377 head_prod2[0] = head_e1;
2378 tail_prod2[0] = tail_e1;
2379 /* imaginary part */
2380 d1 = (double) a_elem[0] * x_elem[1];
2381 d2 = (double) a_elem[1] * x_elem[0];
2382 {
2383 /* Compute double-double = double + double. */
2384 double e, t1, t2;
2385
2386 /* Knuth trick. */
2387 t1 = d1 + d2;
2388 e = t1 - d1;
2389 t2 = ((d2 - e) + (d1 - (t1 - e)));
2390
2391 /* The result is t1 + t2, after normalization. */
2392 head_e1 = t1 + t2;
2393 tail_e1 = t2 - (head_e1 - t1);
2394 }
2395 head_prod2[1] = head_e1;
2396 tail_prod2[1] = tail_e1;
2397 }
2398 {
2399 double head_t, tail_t;
2400 double head_a, tail_a;
2401 double head_b, tail_b;
2402 /* Real part */
2403 head_a = head_sum2[0];
2404 tail_a = tail_sum2[0];
2405 head_b = head_prod2[0];
2406 tail_b = tail_prod2[0];
2407 {
2408 /* Compute double-double = double-double + double-double. */
2409 double bv;
2410 double s1, s2, t1, t2;
2411
2412 /* Add two hi words. */
2413 s1 = head_a + head_b;
2414 bv = s1 - head_a;
2415 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2416
2417 /* Add two lo words. */
2418 t1 = tail_a + tail_b;
2419 bv = t1 - tail_a;
2420 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2421
2422 s2 += t1;
2423
2424 /* Renormalize (s1, s2) to (t1, s2) */
2425 t1 = s1 + s2;
2426 s2 = s2 - (t1 - s1);
2427
2428 t2 += s2;
2429
2430 /* Renormalize (t1, t2) */
2431 head_t = t1 + t2;
2432 tail_t = t2 - (head_t - t1);
2433 }
2434 head_sum2[0] = head_t;
2435 tail_sum2[0] = tail_t;
2436 /* Imaginary part */
2437 head_a = head_sum2[1];
2438 tail_a = tail_sum2[1];
2439 head_b = head_prod2[1];
2440 tail_b = tail_prod2[1];
2441 {
2442 /* Compute double-double = double-double + double-double. */
2443 double bv;
2444 double s1, s2, t1, t2;
2445
2446 /* Add two hi words. */
2447 s1 = head_a + head_b;
2448 bv = s1 - head_a;
2449 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2450
2451 /* Add two lo words. */
2452 t1 = tail_a + tail_b;
2453 bv = t1 - tail_a;
2454 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2455
2456 s2 += t1;
2457
2458 /* Renormalize (s1, s2) to (t1, s2) */
2459 t1 = s1 + s2;
2460 s2 = s2 - (t1 - s1);
2461
2462 t2 += s2;
2463
2464 /* Renormalize (t1, t2) */
2465 head_t = t1 + t2;
2466 tail_t = t2 - (head_t - t1);
2467 }
2468 head_sum2[1] = head_t;
2469 tail_sum2[1] = tail_t;
2470 }
2471 }
2472 {
2473 double head_t, tail_t;
2474 double head_a, tail_a;
2475 double head_b, tail_b;
2476 /* Real part */
2477 head_a = head_sum1[0];
2478 tail_a = tail_sum1[0];
2479 head_b = head_sum2[0];
2480 tail_b = tail_sum2[0];
2481 {
2482 /* Compute double-double = double-double + double-double. */
2483 double bv;
2484 double s1, s2, t1, t2;
2485
2486 /* Add two hi words. */
2487 s1 = head_a + head_b;
2488 bv = s1 - head_a;
2489 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2490
2491 /* Add two lo words. */
2492 t1 = tail_a + tail_b;
2493 bv = t1 - tail_a;
2494 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2495
2496 s2 += t1;
2497
2498 /* Renormalize (s1, s2) to (t1, s2) */
2499 t1 = s1 + s2;
2500 s2 = s2 - (t1 - s1);
2501
2502 t2 += s2;
2503
2504 /* Renormalize (t1, t2) */
2505 head_t = t1 + t2;
2506 tail_t = t2 - (head_t - t1);
2507 }
2508 head_sum1[0] = head_t;
2509 tail_sum1[0] = tail_t;
2510 /* Imaginary part */
2511 head_a = head_sum1[1];
2512 tail_a = tail_sum1[1];
2513 head_b = head_sum2[1];
2514 tail_b = tail_sum2[1];
2515 {
2516 /* Compute double-double = double-double + double-double. */
2517 double bv;
2518 double s1, s2, t1, t2;
2519
2520 /* Add two hi words. */
2521 s1 = head_a + head_b;
2522 bv = s1 - head_a;
2523 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2524
2525 /* Add two lo words. */
2526 t1 = tail_a + tail_b;
2527 bv = t1 - tail_a;
2528 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2529
2530 s2 += t1;
2531
2532 /* Renormalize (s1, s2) to (t1, s2) */
2533 t1 = s1 + s2;
2534 s2 = s2 - (t1 - s1);
2535
2536 t2 += s2;
2537
2538 /* Renormalize (t1, t2) */
2539 head_t = t1 + t2;
2540 tail_t = t2 - (head_t - t1);
2541 }
2542 head_sum1[1] = head_t;
2543 tail_sum1[1] = tail_t;
2544 }
2545 {
2546 /* Compute complex-extra = complex-extra * complex-double. */
2547 double head_a0, tail_a0;
2548 double head_a1, tail_a1;
2549 double head_t1, tail_t1;
2550 double head_t2, tail_t2;
2551 head_a0 = head_sum1[0];
2552 tail_a0 = tail_sum1[0];
2553 head_a1 = head_sum1[1];
2554 tail_a1 = tail_sum1[1];
2555 /* real part */
2556 {
2557 /* Compute double-double = double-double * double. */
2558 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2559
2560 con = head_a0 * split;
2561 a11 = con - head_a0;
2562 a11 = con - a11;
2563 a21 = head_a0 - a11;
2564 con = alpha_i[0] * split;
2565 b1 = con - alpha_i[0];
2566 b1 = con - b1;
2567 b2 = alpha_i[0] - b1;
2568
2569 c11 = head_a0 * alpha_i[0];
2570 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2571
2572 c2 = tail_a0 * alpha_i[0];
2573 t1 = c11 + c2;
2574 t2 = (c2 - (t1 - c11)) + c21;
2575
2576 head_t1 = t1 + t2;
2577 tail_t1 = t2 - (head_t1 - t1);
2578 }
2579 {
2580 /* Compute double-double = double-double * double. */
2581 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2582
2583 con = head_a1 * split;
2584 a11 = con - head_a1;
2585 a11 = con - a11;
2586 a21 = head_a1 - a11;
2587 con = alpha_i[1] * split;
2588 b1 = con - alpha_i[1];
2589 b1 = con - b1;
2590 b2 = alpha_i[1] - b1;
2591
2592 c11 = head_a1 * alpha_i[1];
2593 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2594
2595 c2 = tail_a1 * alpha_i[1];
2596 t1 = c11 + c2;
2597 t2 = (c2 - (t1 - c11)) + c21;
2598
2599 head_t2 = t1 + t2;
2600 tail_t2 = t2 - (head_t2 - t1);
2601 }
2602 head_t2 = -head_t2;
2603 tail_t2 = -tail_t2;
2604 {
2605 /* Compute double-double = double-double + double-double. */
2606 double bv;
2607 double s1, s2, t1, t2;
2608
2609 /* Add two hi words. */
2610 s1 = head_t1 + head_t2;
2611 bv = s1 - head_t1;
2612 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2613
2614 /* Add two lo words. */
2615 t1 = tail_t1 + tail_t2;
2616 bv = t1 - tail_t1;
2617 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2618
2619 s2 += t1;
2620
2621 /* Renormalize (s1, s2) to (t1, s2) */
2622 t1 = s1 + s2;
2623 s2 = s2 - (t1 - s1);
2624
2625 t2 += s2;
2626
2627 /* Renormalize (t1, t2) */
2628 head_t1 = t1 + t2;
2629 tail_t1 = t2 - (head_t1 - t1);
2630 }
2631 head_tmp1[0] = head_t1;
2632 tail_tmp1[0] = tail_t1;
2633 /* imaginary part */
2634 {
2635 /* Compute double-double = double-double * double. */
2636 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2637
2638 con = head_a1 * split;
2639 a11 = con - head_a1;
2640 a11 = con - a11;
2641 a21 = head_a1 - a11;
2642 con = alpha_i[0] * split;
2643 b1 = con - alpha_i[0];
2644 b1 = con - b1;
2645 b2 = alpha_i[0] - b1;
2646
2647 c11 = head_a1 * alpha_i[0];
2648 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2649
2650 c2 = tail_a1 * alpha_i[0];
2651 t1 = c11 + c2;
2652 t2 = (c2 - (t1 - c11)) + c21;
2653
2654 head_t1 = t1 + t2;
2655 tail_t1 = t2 - (head_t1 - t1);
2656 }
2657 {
2658 /* Compute double-double = double-double * double. */
2659 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
2660
2661 con = head_a0 * split;
2662 a11 = con - head_a0;
2663 a11 = con - a11;
2664 a21 = head_a0 - a11;
2665 con = alpha_i[1] * split;
2666 b1 = con - alpha_i[1];
2667 b1 = con - b1;
2668 b2 = alpha_i[1] - b1;
2669
2670 c11 = head_a0 * alpha_i[1];
2671 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
2672
2673 c2 = tail_a0 * alpha_i[1];
2674 t1 = c11 + c2;
2675 t2 = (c2 - (t1 - c11)) + c21;
2676
2677 head_t2 = t1 + t2;
2678 tail_t2 = t2 - (head_t2 - t1);
2679 }
2680 {
2681 /* Compute double-double = double-double + double-double. */
2682 double bv;
2683 double s1, s2, t1, t2;
2684
2685 /* Add two hi words. */
2686 s1 = head_t1 + head_t2;
2687 bv = s1 - head_t1;
2688 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2689
2690 /* Add two lo words. */
2691 t1 = tail_t1 + tail_t2;
2692 bv = t1 - tail_t1;
2693 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2694
2695 s2 += t1;
2696
2697 /* Renormalize (s1, s2) to (t1, s2) */
2698 t1 = s1 + s2;
2699 s2 = s2 - (t1 - s1);
2700
2701 t2 += s2;
2702
2703 /* Renormalize (t1, t2) */
2704 head_t1 = t1 + t2;
2705 tail_t1 = t2 - (head_t1 - t1);
2706 }
2707 head_tmp1[1] = head_t1;
2708 tail_tmp1[1] = tail_t1;
2709 }
2710
2711 y_elem[0] = y_i[yi];
2712 y_elem[1] = y_i[yi + 1];
2713 {
2714 /* Compute complex-extra = complex-double * complex-double. */
2715 double head_t1, tail_t1;
2716 double head_t2, tail_t2;
2717 /* Real part */
2718 {
2719 /* Compute double_double = double * double. */
2720 double a1, a2, b1, b2, con;
2721
2722 con = y_elem[0] * split;
2723 a1 = con - y_elem[0];
2724 a1 = con - a1;
2725 a2 = y_elem[0] - a1;
2726 con = beta_i[0] * split;
2727 b1 = con - beta_i[0];
2728 b1 = con - b1;
2729 b2 = beta_i[0] - b1;
2730
2731 head_t1 = y_elem[0] * beta_i[0];
2732 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
2733 }
2734 {
2735 /* Compute double_double = double * double. */
2736 double a1, a2, b1, b2, con;
2737
2738 con = y_elem[1] * split;
2739 a1 = con - y_elem[1];
2740 a1 = con - a1;
2741 a2 = y_elem[1] - a1;
2742 con = beta_i[1] * split;
2743 b1 = con - beta_i[1];
2744 b1 = con - b1;
2745 b2 = beta_i[1] - b1;
2746
2747 head_t2 = y_elem[1] * beta_i[1];
2748 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
2749 }
2750 head_t2 = -head_t2;
2751 tail_t2 = -tail_t2;
2752 {
2753 /* Compute double-double = double-double + double-double. */
2754 double bv;
2755 double s1, s2, t1, t2;
2756
2757 /* Add two hi words. */
2758 s1 = head_t1 + head_t2;
2759 bv = s1 - head_t1;
2760 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2761
2762 /* Add two lo words. */
2763 t1 = tail_t1 + tail_t2;
2764 bv = t1 - tail_t1;
2765 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2766
2767 s2 += t1;
2768
2769 /* Renormalize (s1, s2) to (t1, s2) */
2770 t1 = s1 + s2;
2771 s2 = s2 - (t1 - s1);
2772
2773 t2 += s2;
2774
2775 /* Renormalize (t1, t2) */
2776 head_t1 = t1 + t2;
2777 tail_t1 = t2 - (head_t1 - t1);
2778 }
2779 head_tmp2[0] = head_t1;
2780 tail_tmp2[0] = tail_t1;
2781 /* Imaginary part */
2782 {
2783 /* Compute double_double = double * double. */
2784 double a1, a2, b1, b2, con;
2785
2786 con = y_elem[1] * split;
2787 a1 = con - y_elem[1];
2788 a1 = con - a1;
2789 a2 = y_elem[1] - a1;
2790 con = beta_i[0] * split;
2791 b1 = con - beta_i[0];
2792 b1 = con - b1;
2793 b2 = beta_i[0] - b1;
2794
2795 head_t1 = y_elem[1] * beta_i[0];
2796 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
2797 }
2798 {
2799 /* Compute double_double = double * double. */
2800 double a1, a2, b1, b2, con;
2801
2802 con = y_elem[0] * split;
2803 a1 = con - y_elem[0];
2804 a1 = con - a1;
2805 a2 = y_elem[0] - a1;
2806 con = beta_i[1] * split;
2807 b1 = con - beta_i[1];
2808 b1 = con - b1;
2809 b2 = beta_i[1] - b1;
2810
2811 head_t2 = y_elem[0] * beta_i[1];
2812 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
2813 }
2814 {
2815 /* Compute double-double = double-double + double-double. */
2816 double bv;
2817 double s1, s2, t1, t2;
2818
2819 /* Add two hi words. */
2820 s1 = head_t1 + head_t2;
2821 bv = s1 - head_t1;
2822 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
2823
2824 /* Add two lo words. */
2825 t1 = tail_t1 + tail_t2;
2826 bv = t1 - tail_t1;
2827 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
2828
2829 s2 += t1;
2830
2831 /* Renormalize (s1, s2) to (t1, s2) */
2832 t1 = s1 + s2;
2833 s2 = s2 - (t1 - s1);
2834
2835 t2 += s2;
2836
2837 /* Renormalize (t1, t2) */
2838 head_t1 = t1 + t2;
2839 tail_t1 = t2 - (head_t1 - t1);
2840 }
2841 head_tmp2[1] = head_t1;
2842 tail_tmp2[1] = tail_t1;
2843 }
2844 {
2845 double head_t, tail_t;
2846 double head_a, tail_a;
2847 double head_b, tail_b;
2848 /* Real part */
2849 head_a = head_tmp1[0];
2850 tail_a = tail_tmp1[0];
2851 head_b = head_tmp2[0];
2852 tail_b = tail_tmp2[0];
2853 {
2854 /* Compute double-double = double-double + double-double. */
2855 double bv;
2856 double s1, s2, t1, t2;
2857
2858 /* Add two hi words. */
2859 s1 = head_a + head_b;
2860 bv = s1 - head_a;
2861 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2862
2863 /* Add two lo words. */
2864 t1 = tail_a + tail_b;
2865 bv = t1 - tail_a;
2866 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2867
2868 s2 += t1;
2869
2870 /* Renormalize (s1, s2) to (t1, s2) */
2871 t1 = s1 + s2;
2872 s2 = s2 - (t1 - s1);
2873
2874 t2 += s2;
2875
2876 /* Renormalize (t1, t2) */
2877 head_t = t1 + t2;
2878 tail_t = t2 - (head_t - t1);
2879 }
2880 head_tmp3[0] = head_t;
2881 tail_tmp3[0] = tail_t;
2882 /* Imaginary part */
2883 head_a = head_tmp1[1];
2884 tail_a = tail_tmp1[1];
2885 head_b = head_tmp2[1];
2886 tail_b = tail_tmp2[1];
2887 {
2888 /* Compute double-double = double-double + double-double. */
2889 double bv;
2890 double s1, s2, t1, t2;
2891
2892 /* Add two hi words. */
2893 s1 = head_a + head_b;
2894 bv = s1 - head_a;
2895 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
2896
2897 /* Add two lo words. */
2898 t1 = tail_a + tail_b;
2899 bv = t1 - tail_a;
2900 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
2901
2902 s2 += t1;
2903
2904 /* Renormalize (s1, s2) to (t1, s2) */
2905 t1 = s1 + s2;
2906 s2 = s2 - (t1 - s1);
2907
2908 t2 += s2;
2909
2910 /* Renormalize (t1, t2) */
2911 head_t = t1 + t2;
2912 tail_t = t2 - (head_t - t1);
2913 }
2914 head_tmp3[1] = head_t;
2915 tail_tmp3[1] = tail_t;
2916 }
2917 y_i[yi] = head_tmp3[0];
2918 y_i[yi + 1] = head_tmp3[1];
2919 }
2920 }
2921
2922 FPU_FIX_STOP;
2923
2924 break;
2925 }
2926 }
2927 } /* end BLAS_zhemv2_c_c_x */
2928