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