1 #include <blas_extended.h>
2 #include <blas_extended_private.h>
3 #include <blas_fpu.h>
BLAS_csymv2_c_s_x(enum blas_order_type order,enum blas_uplo_type uplo,int n,const void * alpha,const void * a,int lda,const float * x_head,const float * x_tail,int incx,const void * beta,void * y,int incy,enum blas_prec_type prec)4 void BLAS_csymv2_c_s_x(enum blas_order_type order, enum blas_uplo_type uplo,
5 int n, const void *alpha, const void *a, int lda,
6 const float *x_head, const float *x_tail, int incx,
7 const void *beta, 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 symmetric 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) float*
42 * Vector x_head
43 *
44 * x_tail (input) float*
45 * Vector x_tail
46 *
47 * incx (input) int
48 * Stride for vector x.
49 *
50 * beta (input) const void*
51 *
52 * y (input) float*
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_csymv2_c_s_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 = x_head;
82 const float *x_tail_i = x_tail;
83 float *y_i = (float *) y;
84 float *alpha_i = (float *) alpha;
85 float *beta_i = (float *) beta;
86 float a_elem[2];
87 float x_elem;
88 float y_elem[2];
89 float prod1[2];
90 float prod2[2];
91 float sum1[2];
92 float sum2[2];
93 float tmp1[2];
94 float tmp2[2];
95 float tmp3[2];
96
97
98
99 /* Test for no-op */
100 if (n <= 0) {
101 return;
102 }
103 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0
104 && (beta_i[0] == 1.0 && beta_i[1] == 0.0)) {
105 return;
106 }
107
108 /* Check for error conditions. */
109 if (n < 0) {
110 BLAS_error(routine_name, -3, n, NULL);
111 }
112 if (lda < n) {
113 BLAS_error(routine_name, -6, n, NULL);
114 }
115 if (incx == 0) {
116 BLAS_error(routine_name, -9, incx, NULL);
117 }
118 if (incy == 0) {
119 BLAS_error(routine_name, -12, incy, NULL);
120 }
121
122 if ((order == blas_colmajor && uplo == blas_upper) ||
123 (order == blas_rowmajor && uplo == blas_lower)) {
124 incai = lda;
125 incaij = 1;
126 incaij2 = lda;
127 } else {
128 incai = 1;
129 incaij = lda;
130 incaij2 = 1;
131 }
132
133
134 incy *= 2;
135 incai *= 2;
136 incaij *= 2;
137 incaij2 *= 2;
138 xi0 = (incx > 0) ? 0 : ((-n + 1) * incx);
139 yi0 = (incy > 0) ? 0 : ((-n + 1) * incy);
140
141
142
143 /* The most general form, y <--- alpha * A * (x_head + x_tail) + beta * y */
144 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
145 sum1[0] = sum1[1] = 0.0;
146 sum2[0] = sum2[1] = 0.0;
147
148 for (j = 0, aij = ai, xi = xi0; j < i; j++, aij += incaij, xi += incx) {
149 a_elem[0] = a_i[aij];
150 a_elem[1] = a_i[aij + 1];
151 x_elem = x_head_i[xi];
152 {
153 prod1[0] = a_elem[0] * x_elem;
154 prod1[1] = a_elem[1] * x_elem;
155 }
156 sum1[0] = sum1[0] + prod1[0];
157 sum1[1] = sum1[1] + prod1[1];
158 x_elem = x_tail_i[xi];
159 {
160 prod2[0] = a_elem[0] * x_elem;
161 prod2[1] = a_elem[1] * x_elem;
162 }
163 sum2[0] = sum2[0] + prod2[0];
164 sum2[1] = sum2[1] + prod2[1];
165 }
166 for (; j < n; j++, aij += incaij2, xi += incx) {
167 a_elem[0] = a_i[aij];
168 a_elem[1] = a_i[aij + 1];
169 x_elem = x_head_i[xi];
170 {
171 prod1[0] = a_elem[0] * x_elem;
172 prod1[1] = a_elem[1] * x_elem;
173 }
174 sum1[0] = sum1[0] + prod1[0];
175 sum1[1] = sum1[1] + prod1[1];
176 x_elem = x_tail_i[xi];
177 {
178 prod2[0] = a_elem[0] * x_elem;
179 prod2[1] = a_elem[1] * x_elem;
180 }
181 sum2[0] = sum2[0] + prod2[0];
182 sum2[1] = sum2[1] + prod2[1];
183 }
184 sum1[0] = sum1[0] + sum2[0];
185 sum1[1] = sum1[1] + sum2[1];
186 {
187 tmp1[0] = sum1[0] * alpha_i[0] - sum1[1] * alpha_i[1];
188 tmp1[1] = sum1[0] * alpha_i[1] + sum1[1] * alpha_i[0];
189 }
190
191 y_elem[0] = y_i[yi];
192 y_elem[1] = y_i[yi + 1];
193 {
194 tmp2[0] = y_elem[0] * beta_i[0] - y_elem[1] * beta_i[1];
195 tmp2[1] = y_elem[0] * beta_i[1] + y_elem[1] * beta_i[0];
196 }
197
198 tmp3[0] = tmp1[0] + tmp2[0];
199 tmp3[1] = tmp1[1] + tmp2[1];
200 y_i[yi] = tmp3[0];
201 y_i[yi + 1] = tmp3[1];
202 }
203
204
205
206 break;
207 }
208
209 case blas_prec_double:
210 case blas_prec_indigenous:{
211
212 int i, j;
213 int xi, yi, xi0, yi0;
214 int aij, ai;
215 int incai;
216 int incaij, incaij2;
217
218 const float *a_i = (float *) a;
219 const float *x_head_i = x_head;
220 const float *x_tail_i = x_tail;
221 float *y_i = (float *) y;
222 float *alpha_i = (float *) alpha;
223 float *beta_i = (float *) beta;
224 float a_elem[2];
225 float x_elem;
226 float y_elem[2];
227 double prod1[2];
228 double prod2[2];
229 double sum1[2];
230 double sum2[2];
231 double tmp1[2];
232 double tmp2[2];
233 double tmp3[2];
234
235
236
237 /* Test for no-op */
238 if (n <= 0) {
239 return;
240 }
241 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0
242 && (beta_i[0] == 1.0 && beta_i[1] == 0.0)) {
243 return;
244 }
245
246 /* Check for error conditions. */
247 if (n < 0) {
248 BLAS_error(routine_name, -3, n, NULL);
249 }
250 if (lda < n) {
251 BLAS_error(routine_name, -6, n, NULL);
252 }
253 if (incx == 0) {
254 BLAS_error(routine_name, -9, incx, NULL);
255 }
256 if (incy == 0) {
257 BLAS_error(routine_name, -12, incy, NULL);
258 }
259
260 if ((order == blas_colmajor && uplo == blas_upper) ||
261 (order == blas_rowmajor && uplo == blas_lower)) {
262 incai = lda;
263 incaij = 1;
264 incaij2 = lda;
265 } else {
266 incai = 1;
267 incaij = lda;
268 incaij2 = 1;
269 }
270
271
272 incy *= 2;
273 incai *= 2;
274 incaij *= 2;
275 incaij2 *= 2;
276 xi0 = (incx > 0) ? 0 : ((-n + 1) * incx);
277 yi0 = (incy > 0) ? 0 : ((-n + 1) * incy);
278
279
280
281 /* The most general form, y <--- alpha * A * (x_head + x_tail) + beta * y */
282 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
283 sum1[0] = sum1[1] = 0.0;
284 sum2[0] = sum2[1] = 0.0;
285
286 for (j = 0, aij = ai, xi = xi0; j < i; j++, aij += incaij, xi += incx) {
287 a_elem[0] = a_i[aij];
288 a_elem[1] = a_i[aij + 1];
289 x_elem = x_head_i[xi];
290 {
291 prod1[0] = (double) a_elem[0] * x_elem;
292 prod1[1] = (double) a_elem[1] * x_elem;
293 }
294 sum1[0] = sum1[0] + prod1[0];
295 sum1[1] = sum1[1] + prod1[1];
296 x_elem = x_tail_i[xi];
297 {
298 prod2[0] = (double) a_elem[0] * x_elem;
299 prod2[1] = (double) a_elem[1] * x_elem;
300 }
301 sum2[0] = sum2[0] + prod2[0];
302 sum2[1] = sum2[1] + prod2[1];
303 }
304 for (; j < n; j++, aij += incaij2, xi += incx) {
305 a_elem[0] = a_i[aij];
306 a_elem[1] = a_i[aij + 1];
307 x_elem = x_head_i[xi];
308 {
309 prod1[0] = (double) a_elem[0] * x_elem;
310 prod1[1] = (double) a_elem[1] * x_elem;
311 }
312 sum1[0] = sum1[0] + prod1[0];
313 sum1[1] = sum1[1] + prod1[1];
314 x_elem = x_tail_i[xi];
315 {
316 prod2[0] = (double) a_elem[0] * x_elem;
317 prod2[1] = (double) a_elem[1] * x_elem;
318 }
319 sum2[0] = sum2[0] + prod2[0];
320 sum2[1] = sum2[1] + prod2[1];
321 }
322 sum1[0] = sum1[0] + sum2[0];
323 sum1[1] = sum1[1] + sum2[1];
324 {
325 tmp1[0] =
326 (double) sum1[0] * alpha_i[0] - (double) sum1[1] * alpha_i[1];
327 tmp1[1] =
328 (double) sum1[0] * alpha_i[1] + (double) sum1[1] * alpha_i[0];
329 }
330 y_elem[0] = y_i[yi];
331 y_elem[1] = y_i[yi + 1];
332 {
333 tmp2[0] =
334 (double) y_elem[0] * beta_i[0] - (double) y_elem[1] * beta_i[1];
335 tmp2[1] =
336 (double) y_elem[0] * beta_i[1] + (double) y_elem[1] * beta_i[0];
337 }
338 tmp3[0] = tmp1[0] + tmp2[0];
339 tmp3[1] = tmp1[1] + tmp2[1];
340 y_i[yi] = tmp3[0];
341 y_i[yi + 1] = tmp3[1];
342 }
343
344
345
346 break;
347 }
348
349 case blas_prec_extra:{
350
351 int i, j;
352 int xi, yi, xi0, yi0;
353 int aij, ai;
354 int incai;
355 int incaij, incaij2;
356
357 const float *a_i = (float *) a;
358 const float *x_head_i = x_head;
359 const float *x_tail_i = x_tail;
360 float *y_i = (float *) y;
361 float *alpha_i = (float *) alpha;
362 float *beta_i = (float *) beta;
363 float a_elem[2];
364 float x_elem;
365 float y_elem[2];
366 double head_prod1[2], tail_prod1[2];
367 double head_prod2[2], tail_prod2[2];
368 double head_sum1[2], tail_sum1[2];
369 double head_sum2[2], tail_sum2[2];
370 double head_tmp1[2], tail_tmp1[2];
371 double head_tmp2[2], tail_tmp2[2];
372 double head_tmp3[2], tail_tmp3[2];
373
374 FPU_FIX_DECL;
375
376 /* Test for no-op */
377 if (n <= 0) {
378 return;
379 }
380 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0
381 && (beta_i[0] == 1.0 && beta_i[1] == 0.0)) {
382 return;
383 }
384
385 /* Check for error conditions. */
386 if (n < 0) {
387 BLAS_error(routine_name, -3, n, NULL);
388 }
389 if (lda < n) {
390 BLAS_error(routine_name, -6, n, NULL);
391 }
392 if (incx == 0) {
393 BLAS_error(routine_name, -9, incx, NULL);
394 }
395 if (incy == 0) {
396 BLAS_error(routine_name, -12, incy, NULL);
397 }
398
399 if ((order == blas_colmajor && uplo == blas_upper) ||
400 (order == blas_rowmajor && uplo == blas_lower)) {
401 incai = lda;
402 incaij = 1;
403 incaij2 = lda;
404 } else {
405 incai = 1;
406 incaij = lda;
407 incaij2 = 1;
408 }
409
410
411 incy *= 2;
412 incai *= 2;
413 incaij *= 2;
414 incaij2 *= 2;
415 xi0 = (incx > 0) ? 0 : ((-n + 1) * incx);
416 yi0 = (incy > 0) ? 0 : ((-n + 1) * incy);
417
418 FPU_FIX_START;
419
420 /* The most general form, y <--- alpha * A * (x_head + x_tail) + beta * y */
421 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
422 head_sum1[0] = head_sum1[1] = tail_sum1[0] = tail_sum1[1] = 0.0;
423 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] = 0.0;
424
425 for (j = 0, aij = ai, xi = xi0; j < i; j++, aij += incaij, xi += incx) {
426 a_elem[0] = a_i[aij];
427 a_elem[1] = a_i[aij + 1];
428 x_elem = x_head_i[xi];
429 {
430 head_prod1[0] = (double) a_elem[0] * x_elem;
431 tail_prod1[0] = 0.0;
432 head_prod1[1] = (double) a_elem[1] * x_elem;
433 tail_prod1[1] = 0.0;
434 }
435 {
436 double head_t, tail_t;
437 double head_a, tail_a;
438 double head_b, tail_b;
439 /* Real part */
440 head_a = head_sum1[0];
441 tail_a = tail_sum1[0];
442 head_b = head_prod1[0];
443 tail_b = tail_prod1[0];
444 {
445 /* Compute double-double = double-double + double-double. */
446 double bv;
447 double s1, s2, t1, t2;
448
449 /* Add two hi words. */
450 s1 = head_a + head_b;
451 bv = s1 - head_a;
452 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
453
454 /* Add two lo words. */
455 t1 = tail_a + tail_b;
456 bv = t1 - tail_a;
457 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
458
459 s2 += t1;
460
461 /* Renormalize (s1, s2) to (t1, s2) */
462 t1 = s1 + s2;
463 s2 = s2 - (t1 - s1);
464
465 t2 += s2;
466
467 /* Renormalize (t1, t2) */
468 head_t = t1 + t2;
469 tail_t = t2 - (head_t - t1);
470 }
471 head_sum1[0] = head_t;
472 tail_sum1[0] = tail_t;
473 /* Imaginary part */
474 head_a = head_sum1[1];
475 tail_a = tail_sum1[1];
476 head_b = head_prod1[1];
477 tail_b = tail_prod1[1];
478 {
479 /* Compute double-double = double-double + double-double. */
480 double bv;
481 double s1, s2, t1, t2;
482
483 /* Add two hi words. */
484 s1 = head_a + head_b;
485 bv = s1 - head_a;
486 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
487
488 /* Add two lo words. */
489 t1 = tail_a + tail_b;
490 bv = t1 - tail_a;
491 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
492
493 s2 += t1;
494
495 /* Renormalize (s1, s2) to (t1, s2) */
496 t1 = s1 + s2;
497 s2 = s2 - (t1 - s1);
498
499 t2 += s2;
500
501 /* Renormalize (t1, t2) */
502 head_t = t1 + t2;
503 tail_t = t2 - (head_t - t1);
504 }
505 head_sum1[1] = head_t;
506 tail_sum1[1] = tail_t;
507 }
508 x_elem = x_tail_i[xi];
509 {
510 head_prod2[0] = (double) a_elem[0] * x_elem;
511 tail_prod2[0] = 0.0;
512 head_prod2[1] = (double) a_elem[1] * x_elem;
513 tail_prod2[1] = 0.0;
514 }
515 {
516 double head_t, tail_t;
517 double head_a, tail_a;
518 double head_b, tail_b;
519 /* Real part */
520 head_a = head_sum2[0];
521 tail_a = tail_sum2[0];
522 head_b = head_prod2[0];
523 tail_b = tail_prod2[0];
524 {
525 /* Compute double-double = double-double + double-double. */
526 double bv;
527 double s1, s2, t1, t2;
528
529 /* Add two hi words. */
530 s1 = head_a + head_b;
531 bv = s1 - head_a;
532 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
533
534 /* Add two lo words. */
535 t1 = tail_a + tail_b;
536 bv = t1 - tail_a;
537 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
538
539 s2 += t1;
540
541 /* Renormalize (s1, s2) to (t1, s2) */
542 t1 = s1 + s2;
543 s2 = s2 - (t1 - s1);
544
545 t2 += s2;
546
547 /* Renormalize (t1, t2) */
548 head_t = t1 + t2;
549 tail_t = t2 - (head_t - t1);
550 }
551 head_sum2[0] = head_t;
552 tail_sum2[0] = tail_t;
553 /* Imaginary part */
554 head_a = head_sum2[1];
555 tail_a = tail_sum2[1];
556 head_b = head_prod2[1];
557 tail_b = tail_prod2[1];
558 {
559 /* Compute double-double = double-double + double-double. */
560 double bv;
561 double s1, s2, t1, t2;
562
563 /* Add two hi words. */
564 s1 = head_a + head_b;
565 bv = s1 - head_a;
566 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
567
568 /* Add two lo words. */
569 t1 = tail_a + tail_b;
570 bv = t1 - tail_a;
571 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
572
573 s2 += t1;
574
575 /* Renormalize (s1, s2) to (t1, s2) */
576 t1 = s1 + s2;
577 s2 = s2 - (t1 - s1);
578
579 t2 += s2;
580
581 /* Renormalize (t1, t2) */
582 head_t = t1 + t2;
583 tail_t = t2 - (head_t - t1);
584 }
585 head_sum2[1] = head_t;
586 tail_sum2[1] = tail_t;
587 }
588 }
589 for (; j < n; j++, aij += incaij2, xi += incx) {
590 a_elem[0] = a_i[aij];
591 a_elem[1] = a_i[aij + 1];
592 x_elem = x_head_i[xi];
593 {
594 head_prod1[0] = (double) a_elem[0] * x_elem;
595 tail_prod1[0] = 0.0;
596 head_prod1[1] = (double) a_elem[1] * x_elem;
597 tail_prod1[1] = 0.0;
598 }
599 {
600 double head_t, tail_t;
601 double head_a, tail_a;
602 double head_b, tail_b;
603 /* Real part */
604 head_a = head_sum1[0];
605 tail_a = tail_sum1[0];
606 head_b = head_prod1[0];
607 tail_b = tail_prod1[0];
608 {
609 /* Compute double-double = double-double + double-double. */
610 double bv;
611 double s1, s2, t1, t2;
612
613 /* Add two hi words. */
614 s1 = head_a + head_b;
615 bv = s1 - head_a;
616 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
617
618 /* Add two lo words. */
619 t1 = tail_a + tail_b;
620 bv = t1 - tail_a;
621 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
622
623 s2 += t1;
624
625 /* Renormalize (s1, s2) to (t1, s2) */
626 t1 = s1 + s2;
627 s2 = s2 - (t1 - s1);
628
629 t2 += s2;
630
631 /* Renormalize (t1, t2) */
632 head_t = t1 + t2;
633 tail_t = t2 - (head_t - t1);
634 }
635 head_sum1[0] = head_t;
636 tail_sum1[0] = tail_t;
637 /* Imaginary part */
638 head_a = head_sum1[1];
639 tail_a = tail_sum1[1];
640 head_b = head_prod1[1];
641 tail_b = tail_prod1[1];
642 {
643 /* Compute double-double = double-double + double-double. */
644 double bv;
645 double s1, s2, t1, t2;
646
647 /* Add two hi words. */
648 s1 = head_a + head_b;
649 bv = s1 - head_a;
650 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
651
652 /* Add two lo words. */
653 t1 = tail_a + tail_b;
654 bv = t1 - tail_a;
655 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
656
657 s2 += t1;
658
659 /* Renormalize (s1, s2) to (t1, s2) */
660 t1 = s1 + s2;
661 s2 = s2 - (t1 - s1);
662
663 t2 += s2;
664
665 /* Renormalize (t1, t2) */
666 head_t = t1 + t2;
667 tail_t = t2 - (head_t - t1);
668 }
669 head_sum1[1] = head_t;
670 tail_sum1[1] = tail_t;
671 }
672 x_elem = x_tail_i[xi];
673 {
674 head_prod2[0] = (double) a_elem[0] * x_elem;
675 tail_prod2[0] = 0.0;
676 head_prod2[1] = (double) a_elem[1] * x_elem;
677 tail_prod2[1] = 0.0;
678 }
679 {
680 double head_t, tail_t;
681 double head_a, tail_a;
682 double head_b, tail_b;
683 /* Real part */
684 head_a = head_sum2[0];
685 tail_a = tail_sum2[0];
686 head_b = head_prod2[0];
687 tail_b = tail_prod2[0];
688 {
689 /* Compute double-double = double-double + double-double. */
690 double bv;
691 double s1, s2, t1, t2;
692
693 /* Add two hi words. */
694 s1 = head_a + head_b;
695 bv = s1 - head_a;
696 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
697
698 /* Add two lo words. */
699 t1 = tail_a + tail_b;
700 bv = t1 - tail_a;
701 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
702
703 s2 += t1;
704
705 /* Renormalize (s1, s2) to (t1, s2) */
706 t1 = s1 + s2;
707 s2 = s2 - (t1 - s1);
708
709 t2 += s2;
710
711 /* Renormalize (t1, t2) */
712 head_t = t1 + t2;
713 tail_t = t2 - (head_t - t1);
714 }
715 head_sum2[0] = head_t;
716 tail_sum2[0] = tail_t;
717 /* Imaginary part */
718 head_a = head_sum2[1];
719 tail_a = tail_sum2[1];
720 head_b = head_prod2[1];
721 tail_b = tail_prod2[1];
722 {
723 /* Compute double-double = double-double + double-double. */
724 double bv;
725 double s1, s2, t1, t2;
726
727 /* Add two hi words. */
728 s1 = head_a + head_b;
729 bv = s1 - head_a;
730 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
731
732 /* Add two lo words. */
733 t1 = tail_a + tail_b;
734 bv = t1 - tail_a;
735 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
736
737 s2 += t1;
738
739 /* Renormalize (s1, s2) to (t1, s2) */
740 t1 = s1 + s2;
741 s2 = s2 - (t1 - s1);
742
743 t2 += s2;
744
745 /* Renormalize (t1, t2) */
746 head_t = t1 + t2;
747 tail_t = t2 - (head_t - t1);
748 }
749 head_sum2[1] = head_t;
750 tail_sum2[1] = tail_t;
751 }
752 }
753 {
754 double head_t, tail_t;
755 double head_a, tail_a;
756 double head_b, tail_b;
757 /* Real part */
758 head_a = head_sum1[0];
759 tail_a = tail_sum1[0];
760 head_b = head_sum2[0];
761 tail_b = tail_sum2[0];
762 {
763 /* Compute double-double = double-double + double-double. */
764 double bv;
765 double s1, s2, t1, t2;
766
767 /* Add two hi words. */
768 s1 = head_a + head_b;
769 bv = s1 - head_a;
770 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
771
772 /* Add two lo words. */
773 t1 = tail_a + tail_b;
774 bv = t1 - tail_a;
775 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
776
777 s2 += t1;
778
779 /* Renormalize (s1, s2) to (t1, s2) */
780 t1 = s1 + s2;
781 s2 = s2 - (t1 - s1);
782
783 t2 += s2;
784
785 /* Renormalize (t1, t2) */
786 head_t = t1 + t2;
787 tail_t = t2 - (head_t - t1);
788 }
789 head_sum1[0] = head_t;
790 tail_sum1[0] = tail_t;
791 /* Imaginary part */
792 head_a = head_sum1[1];
793 tail_a = tail_sum1[1];
794 head_b = head_sum2[1];
795 tail_b = tail_sum2[1];
796 {
797 /* Compute double-double = double-double + double-double. */
798 double bv;
799 double s1, s2, t1, t2;
800
801 /* Add two hi words. */
802 s1 = head_a + head_b;
803 bv = s1 - head_a;
804 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
805
806 /* Add two lo words. */
807 t1 = tail_a + tail_b;
808 bv = t1 - tail_a;
809 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
810
811 s2 += t1;
812
813 /* Renormalize (s1, s2) to (t1, s2) */
814 t1 = s1 + s2;
815 s2 = s2 - (t1 - s1);
816
817 t2 += s2;
818
819 /* Renormalize (t1, t2) */
820 head_t = t1 + t2;
821 tail_t = t2 - (head_t - t1);
822 }
823 head_sum1[1] = head_t;
824 tail_sum1[1] = tail_t;
825 }
826 {
827 double cd[2];
828 cd[0] = (double) alpha_i[0];
829 cd[1] = (double) alpha_i[1];
830 {
831 /* Compute complex-extra = complex-extra * complex-double. */
832 double head_a0, tail_a0;
833 double head_a1, tail_a1;
834 double head_t1, tail_t1;
835 double head_t2, tail_t2;
836 head_a0 = head_sum1[0];
837 tail_a0 = tail_sum1[0];
838 head_a1 = head_sum1[1];
839 tail_a1 = tail_sum1[1];
840 /* real part */
841 {
842 /* Compute double-double = double-double * double. */
843 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
844
845 con = head_a0 * split;
846 a11 = con - head_a0;
847 a11 = con - a11;
848 a21 = head_a0 - a11;
849 con = cd[0] * split;
850 b1 = con - cd[0];
851 b1 = con - b1;
852 b2 = cd[0] - b1;
853
854 c11 = head_a0 * cd[0];
855 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
856
857 c2 = tail_a0 * cd[0];
858 t1 = c11 + c2;
859 t2 = (c2 - (t1 - c11)) + c21;
860
861 head_t1 = t1 + t2;
862 tail_t1 = t2 - (head_t1 - t1);
863 }
864 {
865 /* Compute double-double = double-double * double. */
866 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
867
868 con = head_a1 * split;
869 a11 = con - head_a1;
870 a11 = con - a11;
871 a21 = head_a1 - a11;
872 con = cd[1] * split;
873 b1 = con - cd[1];
874 b1 = con - b1;
875 b2 = cd[1] - b1;
876
877 c11 = head_a1 * cd[1];
878 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
879
880 c2 = tail_a1 * cd[1];
881 t1 = c11 + c2;
882 t2 = (c2 - (t1 - c11)) + c21;
883
884 head_t2 = t1 + t2;
885 tail_t2 = t2 - (head_t2 - t1);
886 }
887 head_t2 = -head_t2;
888 tail_t2 = -tail_t2;
889 {
890 /* Compute double-double = double-double + double-double. */
891 double bv;
892 double s1, s2, t1, t2;
893
894 /* Add two hi words. */
895 s1 = head_t1 + head_t2;
896 bv = s1 - head_t1;
897 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
898
899 /* Add two lo words. */
900 t1 = tail_t1 + tail_t2;
901 bv = t1 - tail_t1;
902 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
903
904 s2 += t1;
905
906 /* Renormalize (s1, s2) to (t1, s2) */
907 t1 = s1 + s2;
908 s2 = s2 - (t1 - s1);
909
910 t2 += s2;
911
912 /* Renormalize (t1, t2) */
913 head_t1 = t1 + t2;
914 tail_t1 = t2 - (head_t1 - t1);
915 }
916 head_tmp1[0] = head_t1;
917 tail_tmp1[0] = tail_t1;
918 /* imaginary part */
919 {
920 /* Compute double-double = double-double * double. */
921 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
922
923 con = head_a1 * split;
924 a11 = con - head_a1;
925 a11 = con - a11;
926 a21 = head_a1 - a11;
927 con = cd[0] * split;
928 b1 = con - cd[0];
929 b1 = con - b1;
930 b2 = cd[0] - b1;
931
932 c11 = head_a1 * cd[0];
933 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
934
935 c2 = tail_a1 * cd[0];
936 t1 = c11 + c2;
937 t2 = (c2 - (t1 - c11)) + c21;
938
939 head_t1 = t1 + t2;
940 tail_t1 = t2 - (head_t1 - t1);
941 }
942 {
943 /* Compute double-double = double-double * double. */
944 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
945
946 con = head_a0 * split;
947 a11 = con - head_a0;
948 a11 = con - a11;
949 a21 = head_a0 - a11;
950 con = cd[1] * split;
951 b1 = con - cd[1];
952 b1 = con - b1;
953 b2 = cd[1] - b1;
954
955 c11 = head_a0 * cd[1];
956 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
957
958 c2 = tail_a0 * cd[1];
959 t1 = c11 + c2;
960 t2 = (c2 - (t1 - c11)) + c21;
961
962 head_t2 = t1 + t2;
963 tail_t2 = t2 - (head_t2 - t1);
964 }
965 {
966 /* Compute double-double = double-double + double-double. */
967 double bv;
968 double s1, s2, t1, t2;
969
970 /* Add two hi words. */
971 s1 = head_t1 + head_t2;
972 bv = s1 - head_t1;
973 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
974
975 /* Add two lo words. */
976 t1 = tail_t1 + tail_t2;
977 bv = t1 - tail_t1;
978 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
979
980 s2 += t1;
981
982 /* Renormalize (s1, s2) to (t1, s2) */
983 t1 = s1 + s2;
984 s2 = s2 - (t1 - s1);
985
986 t2 += s2;
987
988 /* Renormalize (t1, t2) */
989 head_t1 = t1 + t2;
990 tail_t1 = t2 - (head_t1 - t1);
991 }
992 head_tmp1[1] = head_t1;
993 tail_tmp1[1] = tail_t1;
994 }
995
996 }
997 y_elem[0] = y_i[yi];
998 y_elem[1] = y_i[yi + 1];
999 {
1000 double head_e1, tail_e1;
1001 double d1;
1002 double d2;
1003 /* Real part */
1004 d1 = (double) y_elem[0] * beta_i[0];
1005 d2 = (double) -y_elem[1] * beta_i[1];
1006 {
1007 /* Compute double-double = double + double. */
1008 double e, t1, t2;
1009
1010 /* Knuth trick. */
1011 t1 = d1 + d2;
1012 e = t1 - d1;
1013 t2 = ((d2 - e) + (d1 - (t1 - e)));
1014
1015 /* The result is t1 + t2, after normalization. */
1016 head_e1 = t1 + t2;
1017 tail_e1 = t2 - (head_e1 - t1);
1018 }
1019 head_tmp2[0] = head_e1;
1020 tail_tmp2[0] = tail_e1;
1021 /* imaginary part */
1022 d1 = (double) y_elem[0] * beta_i[1];
1023 d2 = (double) y_elem[1] * beta_i[0];
1024 {
1025 /* Compute double-double = double + double. */
1026 double e, t1, t2;
1027
1028 /* Knuth trick. */
1029 t1 = d1 + d2;
1030 e = t1 - d1;
1031 t2 = ((d2 - e) + (d1 - (t1 - e)));
1032
1033 /* The result is t1 + t2, after normalization. */
1034 head_e1 = t1 + t2;
1035 tail_e1 = t2 - (head_e1 - t1);
1036 }
1037 head_tmp2[1] = head_e1;
1038 tail_tmp2[1] = tail_e1;
1039 }
1040 {
1041 double head_t, tail_t;
1042 double head_a, tail_a;
1043 double head_b, tail_b;
1044 /* Real part */
1045 head_a = head_tmp1[0];
1046 tail_a = tail_tmp1[0];
1047 head_b = head_tmp2[0];
1048 tail_b = tail_tmp2[0];
1049 {
1050 /* Compute double-double = double-double + double-double. */
1051 double bv;
1052 double s1, s2, t1, t2;
1053
1054 /* Add two hi words. */
1055 s1 = head_a + head_b;
1056 bv = s1 - head_a;
1057 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1058
1059 /* Add two lo words. */
1060 t1 = tail_a + tail_b;
1061 bv = t1 - tail_a;
1062 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1063
1064 s2 += t1;
1065
1066 /* Renormalize (s1, s2) to (t1, s2) */
1067 t1 = s1 + s2;
1068 s2 = s2 - (t1 - s1);
1069
1070 t2 += s2;
1071
1072 /* Renormalize (t1, t2) */
1073 head_t = t1 + t2;
1074 tail_t = t2 - (head_t - t1);
1075 }
1076 head_tmp3[0] = head_t;
1077 tail_tmp3[0] = tail_t;
1078 /* Imaginary part */
1079 head_a = head_tmp1[1];
1080 tail_a = tail_tmp1[1];
1081 head_b = head_tmp2[1];
1082 tail_b = tail_tmp2[1];
1083 {
1084 /* Compute double-double = double-double + double-double. */
1085 double bv;
1086 double s1, s2, t1, t2;
1087
1088 /* Add two hi words. */
1089 s1 = head_a + head_b;
1090 bv = s1 - head_a;
1091 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1092
1093 /* Add two lo words. */
1094 t1 = tail_a + tail_b;
1095 bv = t1 - tail_a;
1096 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1097
1098 s2 += t1;
1099
1100 /* Renormalize (s1, s2) to (t1, s2) */
1101 t1 = s1 + s2;
1102 s2 = s2 - (t1 - s1);
1103
1104 t2 += s2;
1105
1106 /* Renormalize (t1, t2) */
1107 head_t = t1 + t2;
1108 tail_t = t2 - (head_t - t1);
1109 }
1110 head_tmp3[1] = head_t;
1111 tail_tmp3[1] = tail_t;
1112 }
1113 y_i[yi] = head_tmp3[0];
1114 y_i[yi + 1] = head_tmp3[1];
1115 }
1116
1117 FPU_FIX_STOP;
1118
1119 break;
1120 }
1121 }
1122 } /* end BLAS_csymv2_c_s_x */
1123