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