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