1 #include <blas_extended.h>
2 #include <blas_extended_private.h>
3 #include <blas_fpu.h>
BLAS_zsymv2_z_c_x(enum blas_order_type order,enum blas_uplo_type uplo,int n,const void * alpha,const void * a,int lda,const void * x_head,const void * x_tail,int incx,const void * beta,void * y,int incy,enum blas_prec_type prec)4 void BLAS_zsymv2_z_c_x(enum blas_order_type order, enum blas_uplo_type uplo,
5 int n, const void *alpha, const void *a, int lda,
6 const void *x_head, const void *x_tail, int incx,
7 const void *beta, 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_zsymv2_z_c_x";
70 switch (prec) {
71
72 case blas_prec_single:{
73
74 int i, j;
75 int xi, yi, xi0, yi0;
76 int aij, ai;
77 int incai;
78 int incaij, incaij2;
79
80 const double *a_i = (double *) a;
81 const float *x_head_i = (float *) x_head;
82 const float *x_tail_i = (float *) x_tail;
83 double *y_i = (double *) y;
84 double *alpha_i = (double *) alpha;
85 double *beta_i = (double *) beta;
86 double a_elem[2];
87 float 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 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] =
155 (double) a_elem[0] * x_elem[0] - (double) a_elem[1] * x_elem[1];
156 prod1[1] =
157 (double) a_elem[0] * x_elem[1] + (double) a_elem[1] * x_elem[0];
158 }
159 sum1[0] = sum1[0] + prod1[0];
160 sum1[1] = sum1[1] + prod1[1];
161 x_elem[0] = x_tail_i[xi];
162 x_elem[1] = x_tail_i[xi + 1];
163 {
164 prod2[0] =
165 (double) a_elem[0] * x_elem[0] - (double) a_elem[1] * x_elem[1];
166 prod2[1] =
167 (double) a_elem[0] * x_elem[1] + (double) a_elem[1] * x_elem[0];
168 }
169 sum2[0] = sum2[0] + prod2[0];
170 sum2[1] = sum2[1] + prod2[1];
171 }
172 for (; j < n; j++, aij += incaij2, xi += incx) {
173 a_elem[0] = a_i[aij];
174 a_elem[1] = a_i[aij + 1];
175 x_elem[0] = x_head_i[xi];
176 x_elem[1] = x_head_i[xi + 1];
177 {
178 prod1[0] =
179 (double) a_elem[0] * x_elem[0] - (double) a_elem[1] * x_elem[1];
180 prod1[1] =
181 (double) a_elem[0] * x_elem[1] + (double) a_elem[1] * x_elem[0];
182 }
183 sum1[0] = sum1[0] + prod1[0];
184 sum1[1] = sum1[1] + prod1[1];
185 x_elem[0] = x_tail_i[xi];
186 x_elem[1] = x_tail_i[xi + 1];
187 {
188 prod2[0] =
189 (double) a_elem[0] * x_elem[0] - (double) a_elem[1] * x_elem[1];
190 prod2[1] =
191 (double) a_elem[0] * x_elem[1] + (double) a_elem[1] * x_elem[0];
192 }
193 sum2[0] = sum2[0] + prod2[0];
194 sum2[1] = sum2[1] + prod2[1];
195 }
196 sum1[0] = sum1[0] + sum2[0];
197 sum1[1] = sum1[1] + sum2[1];
198 {
199 tmp1[0] =
200 (double) sum1[0] * alpha_i[0] - (double) sum1[1] * alpha_i[1];
201 tmp1[1] =
202 (double) sum1[0] * alpha_i[1] + (double) sum1[1] * alpha_i[0];
203 }
204 y_elem[0] = y_i[yi];
205 y_elem[1] = y_i[yi + 1];
206 {
207 tmp2[0] =
208 (double) y_elem[0] * beta_i[0] - (double) y_elem[1] * beta_i[1];
209 tmp2[1] =
210 (double) y_elem[0] * beta_i[1] + (double) y_elem[1] * beta_i[0];
211 }
212 tmp3[0] = tmp1[0] + tmp2[0];
213 tmp3[1] = tmp1[1] + tmp2[1];
214 y_i[yi] = tmp3[0];
215 y_i[yi + 1] = tmp3[1];
216 }
217
218
219
220 break;
221 }
222
223 case blas_prec_double:
224 case blas_prec_indigenous:{
225
226 int i, j;
227 int xi, yi, xi0, yi0;
228 int aij, ai;
229 int incai;
230 int incaij, incaij2;
231
232 const double *a_i = (double *) a;
233 const float *x_head_i = (float *) x_head;
234 const float *x_tail_i = (float *) x_tail;
235 double *y_i = (double *) y;
236 double *alpha_i = (double *) alpha;
237 double *beta_i = (double *) beta;
238 double a_elem[2];
239 float x_elem[2];
240 double y_elem[2];
241 double prod1[2];
242 double prod2[2];
243 double sum1[2];
244 double sum2[2];
245 double tmp1[2];
246 double tmp2[2];
247 double tmp3[2];
248
249
250
251 /* Test for no-op */
252 if (n <= 0) {
253 return;
254 }
255 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0
256 && (beta_i[0] == 1.0 && beta_i[1] == 0.0)) {
257 return;
258 }
259
260 /* Check for error conditions. */
261 if (n < 0) {
262 BLAS_error(routine_name, -3, n, NULL);
263 }
264 if (lda < n) {
265 BLAS_error(routine_name, -6, n, NULL);
266 }
267 if (incx == 0) {
268 BLAS_error(routine_name, -9, incx, NULL);
269 }
270 if (incy == 0) {
271 BLAS_error(routine_name, -12, incy, NULL);
272 }
273
274 if ((order == blas_colmajor && uplo == blas_upper) ||
275 (order == blas_rowmajor && uplo == blas_lower)) {
276 incai = lda;
277 incaij = 1;
278 incaij2 = lda;
279 } else {
280 incai = 1;
281 incaij = lda;
282 incaij2 = 1;
283 }
284
285 incx *= 2;
286 incy *= 2;
287 incai *= 2;
288 incaij *= 2;
289 incaij2 *= 2;
290 xi0 = (incx > 0) ? 0 : ((-n + 1) * incx);
291 yi0 = (incy > 0) ? 0 : ((-n + 1) * incy);
292
293
294
295 /* The most general form, y <--- alpha * A * (x_head + x_tail) + beta * y */
296 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
297 sum1[0] = sum1[1] = 0.0;
298 sum2[0] = sum2[1] = 0.0;
299
300 for (j = 0, aij = ai, xi = xi0; j < i; j++, aij += incaij, xi += incx) {
301 a_elem[0] = a_i[aij];
302 a_elem[1] = a_i[aij + 1];
303 x_elem[0] = x_head_i[xi];
304 x_elem[1] = x_head_i[xi + 1];
305 {
306 prod1[0] =
307 (double) a_elem[0] * x_elem[0] - (double) a_elem[1] * x_elem[1];
308 prod1[1] =
309 (double) a_elem[0] * x_elem[1] + (double) a_elem[1] * x_elem[0];
310 }
311 sum1[0] = sum1[0] + prod1[0];
312 sum1[1] = sum1[1] + prod1[1];
313 x_elem[0] = x_tail_i[xi];
314 x_elem[1] = x_tail_i[xi + 1];
315 {
316 prod2[0] =
317 (double) a_elem[0] * x_elem[0] - (double) a_elem[1] * x_elem[1];
318 prod2[1] =
319 (double) a_elem[0] * x_elem[1] + (double) a_elem[1] * x_elem[0];
320 }
321 sum2[0] = sum2[0] + prod2[0];
322 sum2[1] = sum2[1] + prod2[1];
323 }
324 for (; j < n; j++, aij += incaij2, xi += incx) {
325 a_elem[0] = a_i[aij];
326 a_elem[1] = a_i[aij + 1];
327 x_elem[0] = x_head_i[xi];
328 x_elem[1] = x_head_i[xi + 1];
329 {
330 prod1[0] =
331 (double) a_elem[0] * x_elem[0] - (double) a_elem[1] * x_elem[1];
332 prod1[1] =
333 (double) a_elem[0] * x_elem[1] + (double) a_elem[1] * x_elem[0];
334 }
335 sum1[0] = sum1[0] + prod1[0];
336 sum1[1] = sum1[1] + prod1[1];
337 x_elem[0] = x_tail_i[xi];
338 x_elem[1] = x_tail_i[xi + 1];
339 {
340 prod2[0] =
341 (double) a_elem[0] * x_elem[0] - (double) a_elem[1] * x_elem[1];
342 prod2[1] =
343 (double) a_elem[0] * x_elem[1] + (double) a_elem[1] * x_elem[0];
344 }
345 sum2[0] = sum2[0] + prod2[0];
346 sum2[1] = sum2[1] + prod2[1];
347 }
348 sum1[0] = sum1[0] + sum2[0];
349 sum1[1] = sum1[1] + sum2[1];
350 {
351 tmp1[0] =
352 (double) sum1[0] * alpha_i[0] - (double) sum1[1] * alpha_i[1];
353 tmp1[1] =
354 (double) sum1[0] * alpha_i[1] + (double) sum1[1] * alpha_i[0];
355 }
356 y_elem[0] = y_i[yi];
357 y_elem[1] = y_i[yi + 1];
358 {
359 tmp2[0] =
360 (double) y_elem[0] * beta_i[0] - (double) y_elem[1] * beta_i[1];
361 tmp2[1] =
362 (double) y_elem[0] * beta_i[1] + (double) y_elem[1] * beta_i[0];
363 }
364 tmp3[0] = tmp1[0] + tmp2[0];
365 tmp3[1] = tmp1[1] + tmp2[1];
366 y_i[yi] = tmp3[0];
367 y_i[yi + 1] = tmp3[1];
368 }
369
370
371
372 break;
373 }
374
375 case blas_prec_extra:{
376
377 int i, j;
378 int xi, yi, xi0, yi0;
379 int aij, ai;
380 int incai;
381 int incaij, incaij2;
382
383 const double *a_i = (double *) a;
384 const float *x_head_i = (float *) x_head;
385 const float *x_tail_i = (float *) x_tail;
386 double *y_i = (double *) y;
387 double *alpha_i = (double *) alpha;
388 double *beta_i = (double *) beta;
389 double a_elem[2];
390 float x_elem[2];
391 double y_elem[2];
392 double head_prod1[2], tail_prod1[2];
393 double head_prod2[2], tail_prod2[2];
394 double head_sum1[2], tail_sum1[2];
395 double head_sum2[2], tail_sum2[2];
396 double head_tmp1[2], tail_tmp1[2];
397 double head_tmp2[2], tail_tmp2[2];
398 double head_tmp3[2], tail_tmp3[2];
399
400 FPU_FIX_DECL;
401
402 /* Test for no-op */
403 if (n <= 0) {
404 return;
405 }
406 if (alpha_i[0] == 0.0 && alpha_i[1] == 0.0
407 && (beta_i[0] == 1.0 && beta_i[1] == 0.0)) {
408 return;
409 }
410
411 /* Check for error conditions. */
412 if (n < 0) {
413 BLAS_error(routine_name, -3, n, NULL);
414 }
415 if (lda < n) {
416 BLAS_error(routine_name, -6, n, NULL);
417 }
418 if (incx == 0) {
419 BLAS_error(routine_name, -9, incx, NULL);
420 }
421 if (incy == 0) {
422 BLAS_error(routine_name, -12, incy, NULL);
423 }
424
425 if ((order == blas_colmajor && uplo == blas_upper) ||
426 (order == blas_rowmajor && uplo == blas_lower)) {
427 incai = lda;
428 incaij = 1;
429 incaij2 = lda;
430 } else {
431 incai = 1;
432 incaij = lda;
433 incaij2 = 1;
434 }
435
436 incx *= 2;
437 incy *= 2;
438 incai *= 2;
439 incaij *= 2;
440 incaij2 *= 2;
441 xi0 = (incx > 0) ? 0 : ((-n + 1) * incx);
442 yi0 = (incy > 0) ? 0 : ((-n + 1) * incy);
443
444 FPU_FIX_START;
445
446 /* The most general form, y <--- alpha * A * (x_head + x_tail) + beta * y */
447 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
448 head_sum1[0] = head_sum1[1] = tail_sum1[0] = tail_sum1[1] = 0.0;
449 head_sum2[0] = head_sum2[1] = tail_sum2[0] = tail_sum2[1] = 0.0;
450
451 for (j = 0, aij = ai, xi = xi0; j < i; j++, aij += incaij, xi += incx) {
452 a_elem[0] = a_i[aij];
453 a_elem[1] = a_i[aij + 1];
454 x_elem[0] = x_head_i[xi];
455 x_elem[1] = x_head_i[xi + 1];
456 {
457 double cd[2];
458 cd[0] = (double) x_elem[0];
459 cd[1] = (double) x_elem[1];
460 {
461 /* Compute complex-extra = complex-double * complex-double. */
462 double head_t1, tail_t1;
463 double head_t2, tail_t2;
464 /* Real part */
465 {
466 /* Compute double_double = double * double. */
467 double a1, a2, b1, b2, con;
468
469 con = a_elem[0] * split;
470 a1 = con - a_elem[0];
471 a1 = con - a1;
472 a2 = a_elem[0] - a1;
473 con = cd[0] * split;
474 b1 = con - cd[0];
475 b1 = con - b1;
476 b2 = cd[0] - b1;
477
478 head_t1 = a_elem[0] * cd[0];
479 tail_t1 =
480 (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
481 }
482 {
483 /* Compute double_double = double * double. */
484 double a1, a2, b1, b2, con;
485
486 con = a_elem[1] * split;
487 a1 = con - a_elem[1];
488 a1 = con - a1;
489 a2 = a_elem[1] - a1;
490 con = cd[1] * split;
491 b1 = con - cd[1];
492 b1 = con - b1;
493 b2 = cd[1] - b1;
494
495 head_t2 = a_elem[1] * cd[1];
496 tail_t2 =
497 (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
498 }
499 head_t2 = -head_t2;
500 tail_t2 = -tail_t2;
501 {
502 /* Compute double-double = double-double + double-double. */
503 double bv;
504 double s1, s2, t1, t2;
505
506 /* Add two hi words. */
507 s1 = head_t1 + head_t2;
508 bv = s1 - head_t1;
509 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
510
511 /* Add two lo words. */
512 t1 = tail_t1 + tail_t2;
513 bv = t1 - tail_t1;
514 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
515
516 s2 += t1;
517
518 /* Renormalize (s1, s2) to (t1, s2) */
519 t1 = s1 + s2;
520 s2 = s2 - (t1 - s1);
521
522 t2 += s2;
523
524 /* Renormalize (t1, t2) */
525 head_t1 = t1 + t2;
526 tail_t1 = t2 - (head_t1 - t1);
527 }
528 head_prod1[0] = head_t1;
529 tail_prod1[0] = tail_t1;
530 /* Imaginary part */
531 {
532 /* Compute double_double = double * double. */
533 double a1, a2, b1, b2, con;
534
535 con = a_elem[1] * split;
536 a1 = con - a_elem[1];
537 a1 = con - a1;
538 a2 = a_elem[1] - a1;
539 con = cd[0] * split;
540 b1 = con - cd[0];
541 b1 = con - b1;
542 b2 = cd[0] - b1;
543
544 head_t1 = a_elem[1] * cd[0];
545 tail_t1 =
546 (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
547 }
548 {
549 /* Compute double_double = double * double. */
550 double a1, a2, b1, b2, con;
551
552 con = a_elem[0] * split;
553 a1 = con - a_elem[0];
554 a1 = con - a1;
555 a2 = a_elem[0] - a1;
556 con = cd[1] * split;
557 b1 = con - cd[1];
558 b1 = con - b1;
559 b2 = cd[1] - b1;
560
561 head_t2 = a_elem[0] * cd[1];
562 tail_t2 =
563 (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
564 }
565 {
566 /* Compute double-double = double-double + double-double. */
567 double bv;
568 double s1, s2, t1, t2;
569
570 /* Add two hi words. */
571 s1 = head_t1 + head_t2;
572 bv = s1 - head_t1;
573 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
574
575 /* Add two lo words. */
576 t1 = tail_t1 + tail_t2;
577 bv = t1 - tail_t1;
578 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
579
580 s2 += t1;
581
582 /* Renormalize (s1, s2) to (t1, s2) */
583 t1 = s1 + s2;
584 s2 = s2 - (t1 - s1);
585
586 t2 += s2;
587
588 /* Renormalize (t1, t2) */
589 head_t1 = t1 + t2;
590 tail_t1 = t2 - (head_t1 - t1);
591 }
592 head_prod1[1] = head_t1;
593 tail_prod1[1] = tail_t1;
594 }
595 }
596 {
597 double head_t, tail_t;
598 double head_a, tail_a;
599 double head_b, tail_b;
600 /* Real part */
601 head_a = head_sum1[0];
602 tail_a = tail_sum1[0];
603 head_b = head_prod1[0];
604 tail_b = tail_prod1[0];
605 {
606 /* Compute double-double = double-double + double-double. */
607 double bv;
608 double s1, s2, t1, t2;
609
610 /* Add two hi words. */
611 s1 = head_a + head_b;
612 bv = s1 - head_a;
613 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
614
615 /* Add two lo words. */
616 t1 = tail_a + tail_b;
617 bv = t1 - tail_a;
618 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
619
620 s2 += t1;
621
622 /* Renormalize (s1, s2) to (t1, s2) */
623 t1 = s1 + s2;
624 s2 = s2 - (t1 - s1);
625
626 t2 += s2;
627
628 /* Renormalize (t1, t2) */
629 head_t = t1 + t2;
630 tail_t = t2 - (head_t - t1);
631 }
632 head_sum1[0] = head_t;
633 tail_sum1[0] = tail_t;
634 /* Imaginary part */
635 head_a = head_sum1[1];
636 tail_a = tail_sum1[1];
637 head_b = head_prod1[1];
638 tail_b = tail_prod1[1];
639 {
640 /* Compute double-double = double-double + double-double. */
641 double bv;
642 double s1, s2, t1, t2;
643
644 /* Add two hi words. */
645 s1 = head_a + head_b;
646 bv = s1 - head_a;
647 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
648
649 /* Add two lo words. */
650 t1 = tail_a + tail_b;
651 bv = t1 - tail_a;
652 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
653
654 s2 += t1;
655
656 /* Renormalize (s1, s2) to (t1, s2) */
657 t1 = s1 + s2;
658 s2 = s2 - (t1 - s1);
659
660 t2 += s2;
661
662 /* Renormalize (t1, t2) */
663 head_t = t1 + t2;
664 tail_t = t2 - (head_t - t1);
665 }
666 head_sum1[1] = head_t;
667 tail_sum1[1] = tail_t;
668 }
669 x_elem[0] = x_tail_i[xi];
670 x_elem[1] = x_tail_i[xi + 1];
671 {
672 double cd[2];
673 cd[0] = (double) x_elem[0];
674 cd[1] = (double) x_elem[1];
675 {
676 /* Compute complex-extra = complex-double * complex-double. */
677 double head_t1, tail_t1;
678 double head_t2, tail_t2;
679 /* Real part */
680 {
681 /* Compute double_double = double * double. */
682 double a1, a2, b1, b2, con;
683
684 con = a_elem[0] * split;
685 a1 = con - a_elem[0];
686 a1 = con - a1;
687 a2 = a_elem[0] - a1;
688 con = cd[0] * split;
689 b1 = con - cd[0];
690 b1 = con - b1;
691 b2 = cd[0] - b1;
692
693 head_t1 = a_elem[0] * cd[0];
694 tail_t1 =
695 (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
696 }
697 {
698 /* Compute double_double = double * double. */
699 double a1, a2, b1, b2, con;
700
701 con = a_elem[1] * split;
702 a1 = con - a_elem[1];
703 a1 = con - a1;
704 a2 = a_elem[1] - a1;
705 con = cd[1] * split;
706 b1 = con - cd[1];
707 b1 = con - b1;
708 b2 = cd[1] - b1;
709
710 head_t2 = a_elem[1] * cd[1];
711 tail_t2 =
712 (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
713 }
714 head_t2 = -head_t2;
715 tail_t2 = -tail_t2;
716 {
717 /* Compute double-double = double-double + double-double. */
718 double bv;
719 double s1, s2, t1, t2;
720
721 /* Add two hi words. */
722 s1 = head_t1 + head_t2;
723 bv = s1 - head_t1;
724 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
725
726 /* Add two lo words. */
727 t1 = tail_t1 + tail_t2;
728 bv = t1 - tail_t1;
729 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
730
731 s2 += t1;
732
733 /* Renormalize (s1, s2) to (t1, s2) */
734 t1 = s1 + s2;
735 s2 = s2 - (t1 - s1);
736
737 t2 += s2;
738
739 /* Renormalize (t1, t2) */
740 head_t1 = t1 + t2;
741 tail_t1 = t2 - (head_t1 - t1);
742 }
743 head_prod2[0] = head_t1;
744 tail_prod2[0] = tail_t1;
745 /* Imaginary part */
746 {
747 /* Compute double_double = double * double. */
748 double a1, a2, b1, b2, con;
749
750 con = a_elem[1] * split;
751 a1 = con - a_elem[1];
752 a1 = con - a1;
753 a2 = a_elem[1] - a1;
754 con = cd[0] * split;
755 b1 = con - cd[0];
756 b1 = con - b1;
757 b2 = cd[0] - b1;
758
759 head_t1 = a_elem[1] * cd[0];
760 tail_t1 =
761 (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
762 }
763 {
764 /* Compute double_double = double * double. */
765 double a1, a2, b1, b2, con;
766
767 con = a_elem[0] * split;
768 a1 = con - a_elem[0];
769 a1 = con - a1;
770 a2 = a_elem[0] - a1;
771 con = cd[1] * split;
772 b1 = con - cd[1];
773 b1 = con - b1;
774 b2 = cd[1] - b1;
775
776 head_t2 = a_elem[0] * cd[1];
777 tail_t2 =
778 (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
779 }
780 {
781 /* Compute double-double = double-double + double-double. */
782 double bv;
783 double s1, s2, t1, t2;
784
785 /* Add two hi words. */
786 s1 = head_t1 + head_t2;
787 bv = s1 - head_t1;
788 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
789
790 /* Add two lo words. */
791 t1 = tail_t1 + tail_t2;
792 bv = t1 - tail_t1;
793 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
794
795 s2 += t1;
796
797 /* Renormalize (s1, s2) to (t1, s2) */
798 t1 = s1 + s2;
799 s2 = s2 - (t1 - s1);
800
801 t2 += s2;
802
803 /* Renormalize (t1, t2) */
804 head_t1 = t1 + t2;
805 tail_t1 = t2 - (head_t1 - t1);
806 }
807 head_prod2[1] = head_t1;
808 tail_prod2[1] = tail_t1;
809 }
810 }
811 {
812 double head_t, tail_t;
813 double head_a, tail_a;
814 double head_b, tail_b;
815 /* Real part */
816 head_a = head_sum2[0];
817 tail_a = tail_sum2[0];
818 head_b = head_prod2[0];
819 tail_b = tail_prod2[0];
820 {
821 /* Compute double-double = double-double + double-double. */
822 double bv;
823 double s1, s2, t1, t2;
824
825 /* Add two hi words. */
826 s1 = head_a + head_b;
827 bv = s1 - head_a;
828 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
829
830 /* Add two lo words. */
831 t1 = tail_a + tail_b;
832 bv = t1 - tail_a;
833 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
834
835 s2 += t1;
836
837 /* Renormalize (s1, s2) to (t1, s2) */
838 t1 = s1 + s2;
839 s2 = s2 - (t1 - s1);
840
841 t2 += s2;
842
843 /* Renormalize (t1, t2) */
844 head_t = t1 + t2;
845 tail_t = t2 - (head_t - t1);
846 }
847 head_sum2[0] = head_t;
848 tail_sum2[0] = tail_t;
849 /* Imaginary part */
850 head_a = head_sum2[1];
851 tail_a = tail_sum2[1];
852 head_b = head_prod2[1];
853 tail_b = tail_prod2[1];
854 {
855 /* Compute double-double = double-double + double-double. */
856 double bv;
857 double s1, s2, t1, t2;
858
859 /* Add two hi words. */
860 s1 = head_a + head_b;
861 bv = s1 - head_a;
862 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
863
864 /* Add two lo words. */
865 t1 = tail_a + tail_b;
866 bv = t1 - tail_a;
867 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
868
869 s2 += t1;
870
871 /* Renormalize (s1, s2) to (t1, s2) */
872 t1 = s1 + s2;
873 s2 = s2 - (t1 - s1);
874
875 t2 += s2;
876
877 /* Renormalize (t1, t2) */
878 head_t = t1 + t2;
879 tail_t = t2 - (head_t - t1);
880 }
881 head_sum2[1] = head_t;
882 tail_sum2[1] = tail_t;
883 }
884 }
885 for (; j < n; j++, aij += incaij2, xi += incx) {
886 a_elem[0] = a_i[aij];
887 a_elem[1] = a_i[aij + 1];
888 x_elem[0] = x_head_i[xi];
889 x_elem[1] = x_head_i[xi + 1];
890 {
891 double cd[2];
892 cd[0] = (double) x_elem[0];
893 cd[1] = (double) x_elem[1];
894 {
895 /* Compute complex-extra = complex-double * complex-double. */
896 double head_t1, tail_t1;
897 double head_t2, tail_t2;
898 /* Real part */
899 {
900 /* Compute double_double = double * double. */
901 double a1, a2, b1, b2, con;
902
903 con = a_elem[0] * split;
904 a1 = con - a_elem[0];
905 a1 = con - a1;
906 a2 = a_elem[0] - a1;
907 con = cd[0] * split;
908 b1 = con - cd[0];
909 b1 = con - b1;
910 b2 = cd[0] - b1;
911
912 head_t1 = a_elem[0] * cd[0];
913 tail_t1 =
914 (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
915 }
916 {
917 /* Compute double_double = double * double. */
918 double a1, a2, b1, b2, con;
919
920 con = a_elem[1] * split;
921 a1 = con - a_elem[1];
922 a1 = con - a1;
923 a2 = a_elem[1] - a1;
924 con = cd[1] * split;
925 b1 = con - cd[1];
926 b1 = con - b1;
927 b2 = cd[1] - b1;
928
929 head_t2 = a_elem[1] * cd[1];
930 tail_t2 =
931 (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
932 }
933 head_t2 = -head_t2;
934 tail_t2 = -tail_t2;
935 {
936 /* Compute double-double = double-double + double-double. */
937 double bv;
938 double s1, s2, t1, t2;
939
940 /* Add two hi words. */
941 s1 = head_t1 + head_t2;
942 bv = s1 - head_t1;
943 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
944
945 /* Add two lo words. */
946 t1 = tail_t1 + tail_t2;
947 bv = t1 - tail_t1;
948 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
949
950 s2 += t1;
951
952 /* Renormalize (s1, s2) to (t1, s2) */
953 t1 = s1 + s2;
954 s2 = s2 - (t1 - s1);
955
956 t2 += s2;
957
958 /* Renormalize (t1, t2) */
959 head_t1 = t1 + t2;
960 tail_t1 = t2 - (head_t1 - t1);
961 }
962 head_prod1[0] = head_t1;
963 tail_prod1[0] = tail_t1;
964 /* Imaginary part */
965 {
966 /* Compute double_double = double * double. */
967 double a1, a2, b1, b2, con;
968
969 con = a_elem[1] * split;
970 a1 = con - a_elem[1];
971 a1 = con - a1;
972 a2 = a_elem[1] - a1;
973 con = cd[0] * split;
974 b1 = con - cd[0];
975 b1 = con - b1;
976 b2 = cd[0] - b1;
977
978 head_t1 = a_elem[1] * cd[0];
979 tail_t1 =
980 (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
981 }
982 {
983 /* Compute double_double = double * double. */
984 double a1, a2, b1, b2, con;
985
986 con = a_elem[0] * split;
987 a1 = con - a_elem[0];
988 a1 = con - a1;
989 a2 = a_elem[0] - a1;
990 con = cd[1] * split;
991 b1 = con - cd[1];
992 b1 = con - b1;
993 b2 = cd[1] - b1;
994
995 head_t2 = a_elem[0] * cd[1];
996 tail_t2 =
997 (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
998 }
999 {
1000 /* Compute double-double = double-double + double-double. */
1001 double bv;
1002 double s1, s2, t1, t2;
1003
1004 /* Add two hi words. */
1005 s1 = head_t1 + head_t2;
1006 bv = s1 - head_t1;
1007 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1008
1009 /* Add two lo words. */
1010 t1 = tail_t1 + tail_t2;
1011 bv = t1 - tail_t1;
1012 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1013
1014 s2 += t1;
1015
1016 /* Renormalize (s1, s2) to (t1, s2) */
1017 t1 = s1 + s2;
1018 s2 = s2 - (t1 - s1);
1019
1020 t2 += s2;
1021
1022 /* Renormalize (t1, t2) */
1023 head_t1 = t1 + t2;
1024 tail_t1 = t2 - (head_t1 - t1);
1025 }
1026 head_prod1[1] = head_t1;
1027 tail_prod1[1] = tail_t1;
1028 }
1029 }
1030 {
1031 double head_t, tail_t;
1032 double head_a, tail_a;
1033 double head_b, tail_b;
1034 /* Real part */
1035 head_a = head_sum1[0];
1036 tail_a = tail_sum1[0];
1037 head_b = head_prod1[0];
1038 tail_b = tail_prod1[0];
1039 {
1040 /* Compute double-double = double-double + double-double. */
1041 double bv;
1042 double s1, s2, t1, t2;
1043
1044 /* Add two hi words. */
1045 s1 = head_a + head_b;
1046 bv = s1 - head_a;
1047 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1048
1049 /* Add two lo words. */
1050 t1 = tail_a + tail_b;
1051 bv = t1 - tail_a;
1052 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1053
1054 s2 += t1;
1055
1056 /* Renormalize (s1, s2) to (t1, s2) */
1057 t1 = s1 + s2;
1058 s2 = s2 - (t1 - s1);
1059
1060 t2 += s2;
1061
1062 /* Renormalize (t1, t2) */
1063 head_t = t1 + t2;
1064 tail_t = t2 - (head_t - t1);
1065 }
1066 head_sum1[0] = head_t;
1067 tail_sum1[0] = tail_t;
1068 /* Imaginary part */
1069 head_a = head_sum1[1];
1070 tail_a = tail_sum1[1];
1071 head_b = head_prod1[1];
1072 tail_b = tail_prod1[1];
1073 {
1074 /* Compute double-double = double-double + double-double. */
1075 double bv;
1076 double s1, s2, t1, t2;
1077
1078 /* Add two hi words. */
1079 s1 = head_a + head_b;
1080 bv = s1 - head_a;
1081 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1082
1083 /* Add two lo words. */
1084 t1 = tail_a + tail_b;
1085 bv = t1 - tail_a;
1086 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1087
1088 s2 += t1;
1089
1090 /* Renormalize (s1, s2) to (t1, s2) */
1091 t1 = s1 + s2;
1092 s2 = s2 - (t1 - s1);
1093
1094 t2 += s2;
1095
1096 /* Renormalize (t1, t2) */
1097 head_t = t1 + t2;
1098 tail_t = t2 - (head_t - t1);
1099 }
1100 head_sum1[1] = head_t;
1101 tail_sum1[1] = tail_t;
1102 }
1103 x_elem[0] = x_tail_i[xi];
1104 x_elem[1] = x_tail_i[xi + 1];
1105 {
1106 double cd[2];
1107 cd[0] = (double) x_elem[0];
1108 cd[1] = (double) x_elem[1];
1109 {
1110 /* Compute complex-extra = complex-double * complex-double. */
1111 double head_t1, tail_t1;
1112 double head_t2, tail_t2;
1113 /* Real part */
1114 {
1115 /* Compute double_double = double * double. */
1116 double a1, a2, b1, b2, con;
1117
1118 con = a_elem[0] * split;
1119 a1 = con - a_elem[0];
1120 a1 = con - a1;
1121 a2 = a_elem[0] - a1;
1122 con = cd[0] * split;
1123 b1 = con - cd[0];
1124 b1 = con - b1;
1125 b2 = cd[0] - b1;
1126
1127 head_t1 = a_elem[0] * cd[0];
1128 tail_t1 =
1129 (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1130 }
1131 {
1132 /* Compute double_double = double * double. */
1133 double a1, a2, b1, b2, con;
1134
1135 con = a_elem[1] * split;
1136 a1 = con - a_elem[1];
1137 a1 = con - a1;
1138 a2 = a_elem[1] - a1;
1139 con = cd[1] * split;
1140 b1 = con - cd[1];
1141 b1 = con - b1;
1142 b2 = cd[1] - b1;
1143
1144 head_t2 = a_elem[1] * cd[1];
1145 tail_t2 =
1146 (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1147 }
1148 head_t2 = -head_t2;
1149 tail_t2 = -tail_t2;
1150 {
1151 /* Compute double-double = double-double + double-double. */
1152 double bv;
1153 double s1, s2, t1, t2;
1154
1155 /* Add two hi words. */
1156 s1 = head_t1 + head_t2;
1157 bv = s1 - head_t1;
1158 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1159
1160 /* Add two lo words. */
1161 t1 = tail_t1 + tail_t2;
1162 bv = t1 - tail_t1;
1163 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1164
1165 s2 += t1;
1166
1167 /* Renormalize (s1, s2) to (t1, s2) */
1168 t1 = s1 + s2;
1169 s2 = s2 - (t1 - s1);
1170
1171 t2 += s2;
1172
1173 /* Renormalize (t1, t2) */
1174 head_t1 = t1 + t2;
1175 tail_t1 = t2 - (head_t1 - t1);
1176 }
1177 head_prod2[0] = head_t1;
1178 tail_prod2[0] = tail_t1;
1179 /* Imaginary part */
1180 {
1181 /* Compute double_double = double * double. */
1182 double a1, a2, b1, b2, con;
1183
1184 con = a_elem[1] * split;
1185 a1 = con - a_elem[1];
1186 a1 = con - a1;
1187 a2 = a_elem[1] - a1;
1188 con = cd[0] * split;
1189 b1 = con - cd[0];
1190 b1 = con - b1;
1191 b2 = cd[0] - b1;
1192
1193 head_t1 = a_elem[1] * cd[0];
1194 tail_t1 =
1195 (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1196 }
1197 {
1198 /* Compute double_double = double * double. */
1199 double a1, a2, b1, b2, con;
1200
1201 con = a_elem[0] * split;
1202 a1 = con - a_elem[0];
1203 a1 = con - a1;
1204 a2 = a_elem[0] - a1;
1205 con = cd[1] * split;
1206 b1 = con - cd[1];
1207 b1 = con - b1;
1208 b2 = cd[1] - b1;
1209
1210 head_t2 = a_elem[0] * cd[1];
1211 tail_t2 =
1212 (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1213 }
1214 {
1215 /* Compute double-double = double-double + double-double. */
1216 double bv;
1217 double s1, s2, t1, t2;
1218
1219 /* Add two hi words. */
1220 s1 = head_t1 + head_t2;
1221 bv = s1 - head_t1;
1222 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1223
1224 /* Add two lo words. */
1225 t1 = tail_t1 + tail_t2;
1226 bv = t1 - tail_t1;
1227 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1228
1229 s2 += t1;
1230
1231 /* Renormalize (s1, s2) to (t1, s2) */
1232 t1 = s1 + s2;
1233 s2 = s2 - (t1 - s1);
1234
1235 t2 += s2;
1236
1237 /* Renormalize (t1, t2) */
1238 head_t1 = t1 + t2;
1239 tail_t1 = t2 - (head_t1 - t1);
1240 }
1241 head_prod2[1] = head_t1;
1242 tail_prod2[1] = tail_t1;
1243 }
1244 }
1245 {
1246 double head_t, tail_t;
1247 double head_a, tail_a;
1248 double head_b, tail_b;
1249 /* Real part */
1250 head_a = head_sum2[0];
1251 tail_a = tail_sum2[0];
1252 head_b = head_prod2[0];
1253 tail_b = tail_prod2[0];
1254 {
1255 /* Compute double-double = double-double + double-double. */
1256 double bv;
1257 double s1, s2, t1, t2;
1258
1259 /* Add two hi words. */
1260 s1 = head_a + head_b;
1261 bv = s1 - head_a;
1262 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1263
1264 /* Add two lo words. */
1265 t1 = tail_a + tail_b;
1266 bv = t1 - tail_a;
1267 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1268
1269 s2 += t1;
1270
1271 /* Renormalize (s1, s2) to (t1, s2) */
1272 t1 = s1 + s2;
1273 s2 = s2 - (t1 - s1);
1274
1275 t2 += s2;
1276
1277 /* Renormalize (t1, t2) */
1278 head_t = t1 + t2;
1279 tail_t = t2 - (head_t - t1);
1280 }
1281 head_sum2[0] = head_t;
1282 tail_sum2[0] = tail_t;
1283 /* Imaginary part */
1284 head_a = head_sum2[1];
1285 tail_a = tail_sum2[1];
1286 head_b = head_prod2[1];
1287 tail_b = tail_prod2[1];
1288 {
1289 /* Compute double-double = double-double + double-double. */
1290 double bv;
1291 double s1, s2, t1, t2;
1292
1293 /* Add two hi words. */
1294 s1 = head_a + head_b;
1295 bv = s1 - head_a;
1296 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1297
1298 /* Add two lo words. */
1299 t1 = tail_a + tail_b;
1300 bv = t1 - tail_a;
1301 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1302
1303 s2 += t1;
1304
1305 /* Renormalize (s1, s2) to (t1, s2) */
1306 t1 = s1 + s2;
1307 s2 = s2 - (t1 - s1);
1308
1309 t2 += s2;
1310
1311 /* Renormalize (t1, t2) */
1312 head_t = t1 + t2;
1313 tail_t = t2 - (head_t - t1);
1314 }
1315 head_sum2[1] = head_t;
1316 tail_sum2[1] = tail_t;
1317 }
1318 }
1319 {
1320 double head_t, tail_t;
1321 double head_a, tail_a;
1322 double head_b, tail_b;
1323 /* Real part */
1324 head_a = head_sum1[0];
1325 tail_a = tail_sum1[0];
1326 head_b = head_sum2[0];
1327 tail_b = tail_sum2[0];
1328 {
1329 /* Compute double-double = double-double + double-double. */
1330 double bv;
1331 double s1, s2, t1, t2;
1332
1333 /* Add two hi words. */
1334 s1 = head_a + head_b;
1335 bv = s1 - head_a;
1336 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1337
1338 /* Add two lo words. */
1339 t1 = tail_a + tail_b;
1340 bv = t1 - tail_a;
1341 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1342
1343 s2 += t1;
1344
1345 /* Renormalize (s1, s2) to (t1, s2) */
1346 t1 = s1 + s2;
1347 s2 = s2 - (t1 - s1);
1348
1349 t2 += s2;
1350
1351 /* Renormalize (t1, t2) */
1352 head_t = t1 + t2;
1353 tail_t = t2 - (head_t - t1);
1354 }
1355 head_sum1[0] = head_t;
1356 tail_sum1[0] = tail_t;
1357 /* Imaginary part */
1358 head_a = head_sum1[1];
1359 tail_a = tail_sum1[1];
1360 head_b = head_sum2[1];
1361 tail_b = tail_sum2[1];
1362 {
1363 /* Compute double-double = double-double + double-double. */
1364 double bv;
1365 double s1, s2, t1, t2;
1366
1367 /* Add two hi words. */
1368 s1 = head_a + head_b;
1369 bv = s1 - head_a;
1370 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1371
1372 /* Add two lo words. */
1373 t1 = tail_a + tail_b;
1374 bv = t1 - tail_a;
1375 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1376
1377 s2 += t1;
1378
1379 /* Renormalize (s1, s2) to (t1, s2) */
1380 t1 = s1 + s2;
1381 s2 = s2 - (t1 - s1);
1382
1383 t2 += s2;
1384
1385 /* Renormalize (t1, t2) */
1386 head_t = t1 + t2;
1387 tail_t = t2 - (head_t - t1);
1388 }
1389 head_sum1[1] = head_t;
1390 tail_sum1[1] = tail_t;
1391 }
1392 {
1393 /* Compute complex-extra = complex-extra * complex-double. */
1394 double head_a0, tail_a0;
1395 double head_a1, tail_a1;
1396 double head_t1, tail_t1;
1397 double head_t2, tail_t2;
1398 head_a0 = head_sum1[0];
1399 tail_a0 = tail_sum1[0];
1400 head_a1 = head_sum1[1];
1401 tail_a1 = tail_sum1[1];
1402 /* real part */
1403 {
1404 /* Compute double-double = double-double * double. */
1405 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1406
1407 con = head_a0 * split;
1408 a11 = con - head_a0;
1409 a11 = con - a11;
1410 a21 = head_a0 - a11;
1411 con = alpha_i[0] * split;
1412 b1 = con - alpha_i[0];
1413 b1 = con - b1;
1414 b2 = alpha_i[0] - b1;
1415
1416 c11 = head_a0 * alpha_i[0];
1417 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1418
1419 c2 = tail_a0 * alpha_i[0];
1420 t1 = c11 + c2;
1421 t2 = (c2 - (t1 - c11)) + c21;
1422
1423 head_t1 = t1 + t2;
1424 tail_t1 = t2 - (head_t1 - t1);
1425 }
1426 {
1427 /* Compute double-double = double-double * double. */
1428 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1429
1430 con = head_a1 * split;
1431 a11 = con - head_a1;
1432 a11 = con - a11;
1433 a21 = head_a1 - a11;
1434 con = alpha_i[1] * split;
1435 b1 = con - alpha_i[1];
1436 b1 = con - b1;
1437 b2 = alpha_i[1] - b1;
1438
1439 c11 = head_a1 * alpha_i[1];
1440 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1441
1442 c2 = tail_a1 * alpha_i[1];
1443 t1 = c11 + c2;
1444 t2 = (c2 - (t1 - c11)) + c21;
1445
1446 head_t2 = t1 + t2;
1447 tail_t2 = t2 - (head_t2 - t1);
1448 }
1449 head_t2 = -head_t2;
1450 tail_t2 = -tail_t2;
1451 {
1452 /* Compute double-double = double-double + double-double. */
1453 double bv;
1454 double s1, s2, t1, t2;
1455
1456 /* Add two hi words. */
1457 s1 = head_t1 + head_t2;
1458 bv = s1 - head_t1;
1459 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1460
1461 /* Add two lo words. */
1462 t1 = tail_t1 + tail_t2;
1463 bv = t1 - tail_t1;
1464 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1465
1466 s2 += t1;
1467
1468 /* Renormalize (s1, s2) to (t1, s2) */
1469 t1 = s1 + s2;
1470 s2 = s2 - (t1 - s1);
1471
1472 t2 += s2;
1473
1474 /* Renormalize (t1, t2) */
1475 head_t1 = t1 + t2;
1476 tail_t1 = t2 - (head_t1 - t1);
1477 }
1478 head_tmp1[0] = head_t1;
1479 tail_tmp1[0] = tail_t1;
1480 /* imaginary part */
1481 {
1482 /* Compute double-double = double-double * double. */
1483 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1484
1485 con = head_a1 * split;
1486 a11 = con - head_a1;
1487 a11 = con - a11;
1488 a21 = head_a1 - a11;
1489 con = alpha_i[0] * split;
1490 b1 = con - alpha_i[0];
1491 b1 = con - b1;
1492 b2 = alpha_i[0] - b1;
1493
1494 c11 = head_a1 * alpha_i[0];
1495 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1496
1497 c2 = tail_a1 * alpha_i[0];
1498 t1 = c11 + c2;
1499 t2 = (c2 - (t1 - c11)) + c21;
1500
1501 head_t1 = t1 + t2;
1502 tail_t1 = t2 - (head_t1 - t1);
1503 }
1504 {
1505 /* Compute double-double = double-double * double. */
1506 double a11, a21, b1, b2, c11, c21, c2, con, t1, t2;
1507
1508 con = head_a0 * split;
1509 a11 = con - head_a0;
1510 a11 = con - a11;
1511 a21 = head_a0 - a11;
1512 con = alpha_i[1] * split;
1513 b1 = con - alpha_i[1];
1514 b1 = con - b1;
1515 b2 = alpha_i[1] - b1;
1516
1517 c11 = head_a0 * alpha_i[1];
1518 c21 = (((a11 * b1 - c11) + a11 * b2) + a21 * b1) + a21 * b2;
1519
1520 c2 = tail_a0 * alpha_i[1];
1521 t1 = c11 + c2;
1522 t2 = (c2 - (t1 - c11)) + c21;
1523
1524 head_t2 = t1 + t2;
1525 tail_t2 = t2 - (head_t2 - t1);
1526 }
1527 {
1528 /* Compute double-double = double-double + double-double. */
1529 double bv;
1530 double s1, s2, t1, t2;
1531
1532 /* Add two hi words. */
1533 s1 = head_t1 + head_t2;
1534 bv = s1 - head_t1;
1535 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1536
1537 /* Add two lo words. */
1538 t1 = tail_t1 + tail_t2;
1539 bv = t1 - tail_t1;
1540 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1541
1542 s2 += t1;
1543
1544 /* Renormalize (s1, s2) to (t1, s2) */
1545 t1 = s1 + s2;
1546 s2 = s2 - (t1 - s1);
1547
1548 t2 += s2;
1549
1550 /* Renormalize (t1, t2) */
1551 head_t1 = t1 + t2;
1552 tail_t1 = t2 - (head_t1 - t1);
1553 }
1554 head_tmp1[1] = head_t1;
1555 tail_tmp1[1] = tail_t1;
1556 }
1557
1558 y_elem[0] = y_i[yi];
1559 y_elem[1] = y_i[yi + 1];
1560 {
1561 /* Compute complex-extra = complex-double * complex-double. */
1562 double head_t1, tail_t1;
1563 double head_t2, tail_t2;
1564 /* Real part */
1565 {
1566 /* Compute double_double = double * double. */
1567 double a1, a2, b1, b2, con;
1568
1569 con = y_elem[0] * split;
1570 a1 = con - y_elem[0];
1571 a1 = con - a1;
1572 a2 = y_elem[0] - a1;
1573 con = beta_i[0] * split;
1574 b1 = con - beta_i[0];
1575 b1 = con - b1;
1576 b2 = beta_i[0] - b1;
1577
1578 head_t1 = y_elem[0] * beta_i[0];
1579 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1580 }
1581 {
1582 /* Compute double_double = double * double. */
1583 double a1, a2, b1, b2, con;
1584
1585 con = y_elem[1] * split;
1586 a1 = con - y_elem[1];
1587 a1 = con - a1;
1588 a2 = y_elem[1] - a1;
1589 con = beta_i[1] * split;
1590 b1 = con - beta_i[1];
1591 b1 = con - b1;
1592 b2 = beta_i[1] - b1;
1593
1594 head_t2 = y_elem[1] * beta_i[1];
1595 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1596 }
1597 head_t2 = -head_t2;
1598 tail_t2 = -tail_t2;
1599 {
1600 /* Compute double-double = double-double + double-double. */
1601 double bv;
1602 double s1, s2, t1, t2;
1603
1604 /* Add two hi words. */
1605 s1 = head_t1 + head_t2;
1606 bv = s1 - head_t1;
1607 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1608
1609 /* Add two lo words. */
1610 t1 = tail_t1 + tail_t2;
1611 bv = t1 - tail_t1;
1612 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1613
1614 s2 += t1;
1615
1616 /* Renormalize (s1, s2) to (t1, s2) */
1617 t1 = s1 + s2;
1618 s2 = s2 - (t1 - s1);
1619
1620 t2 += s2;
1621
1622 /* Renormalize (t1, t2) */
1623 head_t1 = t1 + t2;
1624 tail_t1 = t2 - (head_t1 - t1);
1625 }
1626 head_tmp2[0] = head_t1;
1627 tail_tmp2[0] = tail_t1;
1628 /* Imaginary part */
1629 {
1630 /* Compute double_double = double * double. */
1631 double a1, a2, b1, b2, con;
1632
1633 con = y_elem[1] * split;
1634 a1 = con - y_elem[1];
1635 a1 = con - a1;
1636 a2 = y_elem[1] - a1;
1637 con = beta_i[0] * split;
1638 b1 = con - beta_i[0];
1639 b1 = con - b1;
1640 b2 = beta_i[0] - b1;
1641
1642 head_t1 = y_elem[1] * beta_i[0];
1643 tail_t1 = (((a1 * b1 - head_t1) + a1 * b2) + a2 * b1) + a2 * b2;
1644 }
1645 {
1646 /* Compute double_double = double * double. */
1647 double a1, a2, b1, b2, con;
1648
1649 con = y_elem[0] * split;
1650 a1 = con - y_elem[0];
1651 a1 = con - a1;
1652 a2 = y_elem[0] - a1;
1653 con = beta_i[1] * split;
1654 b1 = con - beta_i[1];
1655 b1 = con - b1;
1656 b2 = beta_i[1] - b1;
1657
1658 head_t2 = y_elem[0] * beta_i[1];
1659 tail_t2 = (((a1 * b1 - head_t2) + a1 * b2) + a2 * b1) + a2 * b2;
1660 }
1661 {
1662 /* Compute double-double = double-double + double-double. */
1663 double bv;
1664 double s1, s2, t1, t2;
1665
1666 /* Add two hi words. */
1667 s1 = head_t1 + head_t2;
1668 bv = s1 - head_t1;
1669 s2 = ((head_t2 - bv) + (head_t1 - (s1 - bv)));
1670
1671 /* Add two lo words. */
1672 t1 = tail_t1 + tail_t2;
1673 bv = t1 - tail_t1;
1674 t2 = ((tail_t2 - bv) + (tail_t1 - (t1 - bv)));
1675
1676 s2 += t1;
1677
1678 /* Renormalize (s1, s2) to (t1, s2) */
1679 t1 = s1 + s2;
1680 s2 = s2 - (t1 - s1);
1681
1682 t2 += s2;
1683
1684 /* Renormalize (t1, t2) */
1685 head_t1 = t1 + t2;
1686 tail_t1 = t2 - (head_t1 - t1);
1687 }
1688 head_tmp2[1] = head_t1;
1689 tail_tmp2[1] = tail_t1;
1690 }
1691 {
1692 double head_t, tail_t;
1693 double head_a, tail_a;
1694 double head_b, tail_b;
1695 /* Real part */
1696 head_a = head_tmp1[0];
1697 tail_a = tail_tmp1[0];
1698 head_b = head_tmp2[0];
1699 tail_b = tail_tmp2[0];
1700 {
1701 /* Compute double-double = double-double + double-double. */
1702 double bv;
1703 double s1, s2, t1, t2;
1704
1705 /* Add two hi words. */
1706 s1 = head_a + head_b;
1707 bv = s1 - head_a;
1708 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1709
1710 /* Add two lo words. */
1711 t1 = tail_a + tail_b;
1712 bv = t1 - tail_a;
1713 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1714
1715 s2 += t1;
1716
1717 /* Renormalize (s1, s2) to (t1, s2) */
1718 t1 = s1 + s2;
1719 s2 = s2 - (t1 - s1);
1720
1721 t2 += s2;
1722
1723 /* Renormalize (t1, t2) */
1724 head_t = t1 + t2;
1725 tail_t = t2 - (head_t - t1);
1726 }
1727 head_tmp3[0] = head_t;
1728 tail_tmp3[0] = tail_t;
1729 /* Imaginary part */
1730 head_a = head_tmp1[1];
1731 tail_a = tail_tmp1[1];
1732 head_b = head_tmp2[1];
1733 tail_b = tail_tmp2[1];
1734 {
1735 /* Compute double-double = double-double + double-double. */
1736 double bv;
1737 double s1, s2, t1, t2;
1738
1739 /* Add two hi words. */
1740 s1 = head_a + head_b;
1741 bv = s1 - head_a;
1742 s2 = ((head_b - bv) + (head_a - (s1 - bv)));
1743
1744 /* Add two lo words. */
1745 t1 = tail_a + tail_b;
1746 bv = t1 - tail_a;
1747 t2 = ((tail_b - bv) + (tail_a - (t1 - bv)));
1748
1749 s2 += t1;
1750
1751 /* Renormalize (s1, s2) to (t1, s2) */
1752 t1 = s1 + s2;
1753 s2 = s2 - (t1 - s1);
1754
1755 t2 += s2;
1756
1757 /* Renormalize (t1, t2) */
1758 head_t = t1 + t2;
1759 tail_t = t2 - (head_t - t1);
1760 }
1761 head_tmp3[1] = head_t;
1762 tail_tmp3[1] = tail_t;
1763 }
1764 y_i[yi] = head_tmp3[0];
1765 y_i[yi + 1] = head_tmp3[1];
1766 }
1767
1768 FPU_FIX_STOP;
1769
1770 break;
1771 }
1772 }
1773 } /* end BLAS_zsymv2_z_c_x */
1774