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